[閒聊] MCMC

[roussas (...真無言)]

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
推 littlehana:又又真的是太強了~~噓!我沒有叫他又又唷!!