- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page