#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)