Problem with "class bernoulli_distribution"

[gravit123]

hi,

im trying to use "class bernoulli_distribution" to simulate the number of detected particle considering the detector´s deadtime. Trigger probability is 30%. Say from N particles only 0.3N are expected to be initially recorded. And from these 30% entries each ninth particle can be further processed by detector (because of the deadtime after a data processing). So I want to know the survived particle number after data processing at the end. I used the following code which I modified from http://www.cplusplus.com/reference/random/bernoulli_distribution/ (parameter names were not modified).

// bernoulli_distribution
#include <iostream>
#include <random>

int main()
{
	const int nrolls=526;
	
	std::default_random_engine generator;
	std::bernoulli_distribution distribution(0.3);
	
	int count=0;  // count number of trues
	
	for (int i=1; i<nrolls; i++) 
        if (count%9==0 && distribution(generator)) //count%9 counts each ninth
        count++; 
	
	std::cout << "bernoulli_distribution (0.3) x 526:" << std::endl;
	std::cout << "true:  " << count << std::endl;
	std::cout << "false: " << nrolls-count << std::endl;
	
	return 0;
}

the result "ture" is about 20 here (I dont have work PC to tell exact number at the moment, could have it later, sry). But both experiment and formula of deadtime show a number of about 50! Anyway if I correct the trigger probability of 30% up to 70% in std::bernoulli_distribution distribution(0.3);, then the simulated result will agree with other two results. Certainly I´ve tested this with other number than int nrolls=526, but the output seems to agree with a trigger probability of 70%.

I cant find the documentation about "std::bernoulli_distribution distribution()" But is it here really correct to input 0.7 instead of 0.3? Did I overlook or misunderstand something?

thank you in advance

[MiiNiPaa]

if (count%9==0 && distribution(generator)) count++
if count%9 != 0 then no furher action done: no increase in count, no actual generation (because of lazy argument evaluation of &&). So count should be equal 1 at the end of generation

[JLBorges]

// std::default_random_engine generator;
std::default_random_engine generator( std::time(nullptr) ); // *** seed it 

> I cant find the documentation about "std::bernoulli_distribution distribution()"

std::bernoulli_distribution:
http://en.cppreference.com/w/cpp/numeric/random/bernoulli_distribution

C++11 random number generation:
http://en.cppreference.com/w/cpp/numeric/random

[gravit123]

Thanks for your relies

if count%9 != 0 then no furher action done: no increase in count, no actual generation (because of lazy argument evaluation of &&). So count should be equal 1 at the end of generation

My apologies it should be i%9 instead of count%9, then you should have

true: 16

as output.

// bernoulli_distribution
#include <iostream>
#include <random>

int main()
{
	const int nrolls=526;
	
	std::default_random_engine generator;
	std::bernoulli_distribution distribution(0.3);
	
	int count=0;  // count number of trues
	
	for (int i=1; i<nrolls; i++) 
        if (i%9==0 && distribution(generator)) //i%9 counts each ninth
        count++; 
	
	std::cout << "bernoulli_distribution (0.3) x 526:" << std::endl;
	std::cout << "true:  " << count << std::endl;
	std::cout << "false: " << nrolls-count << std::endl;
	
	return 0;
}

0.25 from std::bernoulli_distribution d(0.25); is the probability of truth as shown in example from http://en.cppreference.com/w/cpp/numeric/random/bernoulli_distribution

...well somehow not in my case. As said if I correct 0.3 to 1-0.3=0.7, then the result is tolerably. Its really weird for me.

[MiiNiPaa]

It is expected: && operator works if both operands is true, so it will have a chance of checking result of generation every 9 iterations. It is 526/9 = 58 checks. And from these 58 checks only 30% will evaluate to true: 58 * 0.3 ≈ 17.5. Kinda close to what you got. And then you deduct those 17 or similar from total iterations. Even those where you do not actually check your generator state.

[cire]

Count is incremented only when i is evenly divisible by 9 and distribution(generator) evaluates to true.

As already indicated, the only time that distribution(generator) is evaluated (ie, a random number is generated) is when i is evenly divisible by 9.

So your number of rolls (526) is wrong with regard to the actual number of random numbers generated (which would be about 526/9.)

Try:

    for ( unsigned i=0; i<nrolls; ++i )
        if ( distribution(generator) || i%9 == 0 )
            ++count ;

[gravit123]

Even those where you do not actually check your generator state.

yeah, i I see your point now. If the first condition i%9 does not fit, the second will not be evaluated. ...So we have to evaluate the entries of truth 526*0.3=158 instead of generating 526/9=58 during initialisation?

Using

   for ( unsigned i=0; i<nrolls; ++i )
        if ( distribution(generator) || i%9 == 0 )
            ++count ;

returns

true: 197

Its even more than over all recorded number regardless of deadtime 526*0.3=158 (simulation returns 160). But experiment (and also the well-known formula) shows a value around 46. In that case we have gained extra entries from second expression i%9==0 additionally to first one.

[JLBorges]

#include <iostream>
#include <random>
#include <ctime>
#include <iomanip>
 
int main()
{
    const int nrolls = 526 ;
    const int divisor = 9 ;
    const int ntrials = 4 ;
    std::mt19937 generator( std::time(nullptr) ); // mersenne twister
 
    std::cout << "bernoulli distribution - deadtime distortion (divided sample)\n\n" ;
    std::cout << std::fixed << std::setprecision(2) ;
 
    for( double p = 0.1 ; p < 0.95 ; p += 0.1 )
    {
        for( int n = 0 ; n < ntrials ; ++n )
        {
            std::cout << "bernoulli (" << p << ") x " << nrolls << '/' << divisor << ": " ;
            std::bernoulli_distribution distribution(p);
 
            int count = {0} ;  // count number of trues
            for( int i=1 ; i <= nrolls ; ++i )
                if( distribution(generator) && i%divisor == 0 ) ++count ;
 
            std::cout << std::setw(3) << count<< " (" << std::setw(5)
                      << count * 100.0 / ( nrolls / divisor ) << "%)\n" ;
        }
        std::cout << '\n' ;
    }
}

http://ideone.com/jMXa24