[Help-gsl] Re: Discrete PRNG

[Rodney Sparapani]

On 01/10/11 03:36 PM, Rodney Sparapani wrote:
> Hi Ralph:
>
> With my naive coding of Marsaglia's method, I find it is about twice
> as fast as Walker for K=3. Here's a naive snippet:
>
> double
> w_r[] ={ 0.2522, 0.58523, 0.16257},
> sd_r[]={ 1.10140819, 1.730751282, 2.746961958},
> t_r[] ={ 0.82433435, 0.3338340845, 0.1325240531};
>
> double prob[3], tot=0.;
> int j, k, k1, k2, k3, d1k=0, d1[3], d2k=0, d2[3], d3k=0, d3[3];
>
> for(j=0; j<3; j++){
> prob[j]=exp(-0.5*t_r[j]*pow(y_temp-log_lambda_i, 2.))*w_r[j]/sd_r[j];
> tot+=prob[j];
> }
>
> for(j=0; j<3; j++){
> prob[j]/=tot;
>
> k1=floor(10*prob[j]);
> d1k+=k1;
> d1[j]=k1;
>
> k2=floor(100*prob[j])-10*k1;
> d2k+=k2;
> d2[j]=k2;
>
> k3=floor(1000*prob[j])-100*k1-10*k2;
> d3k+=k3;
> d3[j]=k3;
> }
>
> unif=gsl_ran_flat(r0, 0., 1.);
>
> if (unif< (d1k/10.)) {
> k=gsl_ran_flat(r0, 1., d1k+1. );
>
> if(k<=d1[0]) k=0;
> else if(k<=(d1[0]+d1[1])) k=1;
> else k=2;
> }
> else if (unif< ((d1k/10.)+(d2k/100.)) ) {
> k=gsl_ran_flat(r0, 1., d2k+1. );
>
> if(k<=d2[0]) k=0;
> else if(k<=(d2[0]+d2[1])) k=1;
> else k=2;
> }
> else {
> k=gsl_ran_flat(r0, 1., d3k+1. );
>
> if(k<=d3[0]) k=0;
> else if(k<=(d3[0]+d3[1])) k=1;
> else k=2;
> }
>
> vs.
> double prob[3];
> int j, k;
>
> for(j=0; j<3; j++){
> prob[j]=exp(-0.5*t_r[j]*pow(y_temp-log_lambda_i, 2.))*w_r[j]/sd_r[j];
> }
>
> gsl_ran_discrete_t *G=gsl_ran_discrete_preproc(3, prob);
>
> k=gsl_ran_discrete(r0, G);
>
> gsl_ran_discrete_free(G);

Actually, I found just generating from the Multinomial distribution
with an N=1 is faster than either of the above.  DUH.

Rodney