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 _______________________________________________ Help-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/help-gsl
