作者roussas (...真无言)
看板NCTU-STAT94G
标题[闲聊] MCMC
时间Thu Aug 11 05:08:58 2005
1.选取的fv(x)和收敛次数有关,选到好的函数收敛速度会加快,选到不好的跑1000000次还
是不见得有收敛
2.我想不到更省记忆体的方法,所以大家就动动脑吧....我只能想到省下Z的前I_BURN个!
3.如果用imsl去生成乱数,建议用自定矩阵去存值,这样才能在不用时把记忆体free掉
4.参考资料:casella berger p.255...比那篇papper容易懂....
5.最重要的...大家有发现更好的方法或是需要注意的地方记得讲一下喔,还是我有错的话
也要讲~~XD
//--------------------------------------------
#define Min(x,y) (x<y)?x:y
Input sample size of Chi-suqare ,N
Input the degree of freedom of Chi-square(data type=int),degree
exp_N=N+I_BURN;//I_BURN=前面未收敛,必须要被舍弃的值
double *Z//乱数值所存放的阵列,lenth=N
double *U//随机U ~uniform存放的阵列 lenth=exp_N
double *V//随机V ~fv(x)存放的阵列 lenth=exp_N
//生成V,U二组乱数
imsls_d_random_uniform(exp_N,IMSLS_RETURN_USER,U,0);
生成V ~fv(x),共exp_N个;
//用MCMC生成数列Z
tempZ= V[0];
count=0;
do//空转I_BURN次,不写入Z值
{
count++;
temp= fy(V[count])*fv(V[count-1])/(fv(V[count])*fy(V[count-1]);
alpha= Min(temp,1);
if(U[count]<alpha)
{
tempZ=V[count];
}//if else tempZ 不变
}while(count!=I_BURN);
//接着真正取值,下面do的中的 count=I_BURN开始
Z[0]=tempZ;
i=1;
do
{
count++;//这是算Y,U的位置
temp= fy(V[count])*fv(V[count-1])/(fv(V[count])*fy(V[count-1 ]));
alpha= Min(temp,1);
if(U[count]<alpha)
{
Z[i]=V[count];
}else
{
Z[i]=Z[i-1];
}
//printf("Z[%d]=%lf\n",i,Z[i]);
i++;//这是算Z的位置
}while(count!=exp_N);
free(U);
free(V);
free(Z);
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.113.129.30
1F：推 littlehana:又又真的是太强了~~嘘!我没有叫他又又唷!! 140.113.68.221 08/11