ガンマ乱数発生のプログラム
#include <oxstd.h> main(){ decl alpha,i,u1,u2,x1,x2,a,b,c,x,y; // set the value of alpha alpha=0.5; // ranseed(100); // change seed for(i=0;i<50;++i) { if(alpha > 1.0){ a=((alpha-1)/exp(1))^((alpha-1)/2); b=((alpha+1)/exp(1))^((alpha+1)/2); do{ u1=ranu(1,1);u2=ranu(1,1); x1=u1*a;x2=u2*b; y=x2/x1; }while(2*log(x1)-(alpha-1)*log(y)+y >1); } else if(alpha==1.0){ u1=ranu(1,1);y=-log(u1); } else{ do{ u1=ranu(1,1);u2=ranu(1,1); c=1+alpha/exp(1); if(u1 <=1/c){ x=(c*u1)^(1/alpha);y=exp(-x); } else{ x=-log(c*(1-u1)/alpha);y=x^(alpha-1); } } while(u2 >=y); } println(y); } }