作者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