棄却サンプリングのプログラム

#include <oxstd.h>
#include <oxfloat.h>
#include <oxprob.h>
main(){
decl x,y,c,n,i;
decl n_accept,n_count;
decl file;
 n=1000;
 y=zeros(n,1);
 n_accept=n_count=0;
 for(i=0;i<n;++i){
    do{
        n_count+=1;
        x=ranchi(1,1,1);
    c=sqrt(M_2PI)*3^(1.5)*exp(-1.5);
        }
        while(ranu(1,1) > x*exp(-x)/(c*1.0/sqrt(M_2PI*x)*exp(-0.5*x)));
        n_accept+=1;y[i]=x;
 }
// Output
 format(600);
 println("Acceptance Probability= ",n_accept/n_count);
 file = fopen("gamma.dat", "w");
    fprint(file,"%10.5f",y);fclose(file);   
}