作者markoo (上行下笑阿..)
看板NCTU-STAT98G
標題0728 統算 隨機分配生成 (離散部分)
時間Sat Aug 8 10:47:45 2009
同學，補上之前離散型分配 變數生成
那其他的就同學自行揣摩補上了:D
Bernoulli(p) :
P{X=1} = p,
P{X=0} = 1-p,
X=0 X=1
|_______|_____________|
0 1-p p 1
所以經由uniform(0,1)的落點可以得到一個樣本值X1=0或1, 直到生成Xn結束
double *rand_bernoulli(int n, double p, int seed) //Bernoulli(p)
{
int i;
double *ans;
double *unif;
ans = (double *) malloc (n*sizeof(double));
unif = rand_uniform(n,0,1,seed);
for (i=0;i<n;i++)
{
if (unif[i]<=p) ans[i]=1;
else ans[i]=0;
}
free(unif);
return ans;
}
Binomial(m,p) 可以直接由m個Bernoulli(p)加起來
double *rand_binomial(int n,int m, double p, int seed) //Binomial(m,p)
{
int i,j;
double *ans;
double *bernoulli;
ans = (double *) malloc (n*sizeof(double));
bernoulli = rand_bernoulli(n*m,p,seed);
for (i=0;i<n;i++)
{
ans[i]=0;
for (j=0;j<m;j++) ans[i]+=bernoulli[j+i*m];
}
free(bernoulli);
return ans;
}
這邊Geometric值由1,2,...,infinity.(X代表失敗次數)
X=1 X=2 X=3 .... X=infinity
|________|_________|_________|_____|____|__|_||||||||||
0 p (1-p)p (1-p)^2*p 1
一樣由uniform(0,1)去判斷他落在哪個區域就可以得到一個樣本
double *rand_geometric(int n, double p, int seed) // Geometric(p)
{
int i,k;
int check=0;
double *unif;
double *ans;
double pleft = 0;
double pright= 0;
double temp;
ans = (double *) malloc (n*sizeof(double));
unif = rand_uniform(n,0,1,seed);
for (i=0;i<n;i++)
{
pright =0;
temp = p;
k=0;
do
{
pleft =pright;
pright = pleft+ temp;
k++;
temp = (1-p)*temp;
//這邊temp是P{X=k+1}與P{X=k}之間的比例
} while (unif[i]<= pleft || unif[i]>pright);
ans[i]= k;
}
return ans;
}
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.166.54.115