Question for vdRngGaussian and vdRngUniform.

[Fiori]

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.

[Andrey_N_Intel · Accepted solution]

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