#include <stdlib.h> #include <math.h> #define pi 3.1415926 double rand_mean(int n, double *tmp) { int i; double ans=0; for (i=0;i<n;i++) ans+=tmp[i]; ans = ans/(n*1.0); return ans; } double rand_var(int n,double *tmp) { int i; double ans=0; double mean; mean = rand_mean(n,tmp); for (i=0;i<n;i++) ans+=(tmp[i]-mean)*(tmp[i]-mean); ans = ans/(n*1.0); return ans; } double *rand_uniform(int n, double a,double b, int seed) { int i; double *ans; ans = (double *) malloc (n*sizeof(double)); for (i=0;i<n;i++) { ans[i] = (rand()*seed+1)%32768/32768.0; ans[i] = (b-a)*ans[i]+a; } return ans; } double *rand_normal(int n,double mean, double var,int seed)//輸入var為variance { int i; double *unif; double *ans; ans = (double *) malloc (n*sizeof(double)); unif = rand_uniform(2*n,0,1,seed); for (i=0;i<n;i++) { ans[i] = cos(2*pi*unif[i])*sqrt(-2*log(unif[n+i])); ans[i] = sqrt(var)*ans[i]+mean; } free(unif); return ans; } double *rand_exponential(int n,double lambda, int seed)//輸入lambda為mean的倒數 { int i; double *unif; double *ans; ans = (double *) malloc (n*sizeof(*ans)); unif = rand_uniform(n,0,1,seed); for (i=0;i<n;i++) ans[i] = -log(unif[i])/lambda; return ans; } double *rand_normal2(int n,double mean,double var,int seed) { int i; double *rexp; double *unif; double *ans; int count=0; ans = (double *) malloc (n*sizeof(*ans)); do { rexp = rand_exponential(1,1,seed); unif = rand_uniform(1,0,1,seed); if (unif[0]<=exp(-rexp[0]*rexp[0]/2+rexp[0]-1/2) ) { ans[count] = rexp[0]; count ++; } free(rexp); free(unif); }while (count <n); unif = rand_uniform(n,0,1,seed); for (i=0;i<n;i++) { if (unif[i]<=0.5) ans[i] = -ans[i]; ans[i] = sqrt(var)*ans[i]+mean; } free(unif); return ans; } double *rand_normal3(int n,double mean,double var,int seed) { int i; double *rexp; double *unif; double *ans; int count=0; ans = (double *) malloc (n*sizeof(*ans)); do { rexp = rand_exponential(2,1,seed); if (rexp[1]>=(rexp[0]-1)*(rexp[0]-1)/2 ) { ans[count] = rexp[0]; count ++; } free(rexp); }while (count <n); unif = rand_uniform(n,0,1,seed); for (i=0;i<n;i++) { if (unif[i]<=0.5) ans[i] = -ans[i]; ans[i] = sqrt(var)*ans[i]+mean; } free(unif); return ans; } double *rand_chisquare(int n, int df, int seed)//df:degree of freedom { int i,j; double *rnorm; double *ans; ans = (double *) malloc (n*sizeof(double)); rnorm = rand_normal(n*df,0,1,seed); for (i=0;i<n;i++) { ans[i]=0; for (j=0;j<df;j++) { ans[i]+=rnorm[j+i*df]*rnorm[j+i*df]; } } return ans; } double *rand_t(int n,int df,int seed) { int i; double *ans; double *rchi; double *rnorm; ans = (double *) malloc (n*sizeof(*ans)); rchi = rand_chisquare(n,df,seed); rnorm = rand_normal(n,0,1,seed); for (i=0;i<n;i++) ans[i] = rnorm[i]/sqrt(rchi[i]/df); return ans; } ※ 編輯: markoo 來自: 140.113.7.248 (07/29 16:43)