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