Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6839 Discussions

## Question for vdRngGaussian and vdRngUniform.

Beginner
210 Views

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.

1 Solution
Employee
210 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

Employee
211 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