Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Fiori
Beginner
70 Views

Question for vdRngGaussian and vdRngUniform.

Jump to solution

Hi!
I want to write the following algorithm using the functions vdRngGaussian and vdRngUniform.

Let X be a vector of n Bernoulli(phi) random variables.
sumx = sum(X);
accept<--aceept


for I=1:nN


1. Propose thetacan ~ Normal(theta, sigma2), wher theta = log(1 + phi) - log(1 -phi);


2.  phican = (exp(thetacan)-1)/(exp(thetacan)+1);


3.  Compute
           logcan = sumx*log(phican) + (n -sumx)*log(1-phican) + thetacan -2*log(1 +exp(thetacan));
           logold =  sumx*log(phi) + (n -sumx)*log(1-phi)+ theta-2*log(1 +exp(theta));
           logf = logcan - logold;

4. Propose u ~ Uniform(0, 1)

5. if log(u)<logf
    phi <--- phican;
    accept <-- accept  + 1
end

end of iterations

 

The criterion in order to choose sigma2 is the acceptance rate (accept/N) to range between 30%-50%.
If acceptance rate<30%, increase sigma2.
If acceptance/rata>30%, decrease sigma2.

I simulated n=5000 observations from Bernoulli(0.9). Using Matlab, I chose sigma2 to be  0.01. I run for
this dataset my matlab-code many times and the acceptance rate was around to 45%. The problem I face is that I can't specify sigma2
for the C-code. For example, I chose sigma2 to be 0.00001 and for 5 different iterations of the algorithm the rate was
 {1, 0.2290, 0.0206, 0.3550, 0.3550}.
This is a weird result because from theory it is known than for a value of sigma2, the rate should be equal to a number (for example 0.2290)
and not to range from 0.0206 to 1.

I attached my C-code, Matlab-code and the dataset.

Could you please tell me if I haven't use the functions properly.

Thank you very much.

 

0 Kudos
1 Solution
Andrey_N_Intel
Employee
70 Views

Hello,

In the loop below you call a codec() routine multiple times. Inside of the routine you use the same MT19937 BRNG initialized with time() output. Can you check how the value of seed varies between calls of the codec()? Alternatively, for experimental purposes can you pass the index of the iteration (or whatever you want but different for each iteration) into codec() and use it as the seed? Please, let me know, the results of the experiment.

% % C
sigma2c = 0.0001;
phi = phitrue;
drawsc = nan(1, numits);
for i=1:numits
   
    out = codec(n, x, phi, sigma2c, acceptc);
    acceptc = out(1, 3); % update accept
    phi = out(1, 2); % update phi
   
    drawsc(1, i) = phi;
end
acceptc/numits
 

Aside observation - use of Intel MKL RNGs on vector length 1 is not effective from performance perspective.

Also, just curious: you generate the array of Bernoulli distributed random numbers using Binomial RNG with ntrials=1 and then calculate the sum of the array entries, thus, get Binomial distribution again. Why not to generate the value using one call to Binomial RNG with the same success probability and number of trials 5000? If this case you will not need summation operation. Please, let me know if that makes sense to you.

Andrey

 

View solution in original post

1 Reply
Andrey_N_Intel
Employee
71 Views

Hello,

In the loop below you call a codec() routine multiple times. Inside of the routine you use the same MT19937 BRNG initialized with time() output. Can you check how the value of seed varies between calls of the codec()? Alternatively, for experimental purposes can you pass the index of the iteration (or whatever you want but different for each iteration) into codec() and use it as the seed? Please, let me know, the results of the experiment.

% % C
sigma2c = 0.0001;
phi = phitrue;
drawsc = nan(1, numits);
for i=1:numits
   
    out = codec(n, x, phi, sigma2c, acceptc);
    acceptc = out(1, 3); % update accept
    phi = out(1, 2); % update phi
   
    drawsc(1, i) = phi;
end
acceptc/numits
 

Aside observation - use of Intel MKL RNGs on vector length 1 is not effective from performance perspective.

Also, just curious: you generate the array of Bernoulli distributed random numbers using Binomial RNG with ntrials=1 and then calculate the sum of the array entries, thus, get Binomial distribution again. Why not to generate the value using one call to Binomial RNG with the same success probability and number of trials 5000? If this case you will not need summation operation. Please, let me know if that makes sense to you.

Andrey

 

View solution in original post

Reply