Re: [Help-gsl] Discrete PRNG

[Rodney Sparapani]

On 01/10/11 03:00 PM, X Statistics wrote:
> Hi,
>
> I do not have it but I am interested if you code it up. Please post 
> here if you do.
>
> In my applications (sampling importance resampling) the number of 
> discrete values and the draws are large, say greater than 500. 
> Walker's method (available in GSL) works pretty well for that case. 
> Nowadays, however, I am using a stratified sampling in which just one 
> uniform random variate must be draw to obtain a sample of any size M 
> >= 1. That is fast and induces less variability in the estimation of 
> functions or parameters.
>
> Ralph.
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);


--
Rodney Sparapani       Center for Patient Care and Outcomes Research
Sr. Biostatistician               http://www.mcw.edu/pcor
4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA