※ 引述《roussas (...真無言)》之銘言：
我也有看第一題的說法
他好像是先設一個thata初值為1
然後取值一次取一個U U從UNIFORM(0,N)中取
做出alpha值 再抽Z從U(0,1)中取跟alpha比較
好像不太像一次抽許多點出來
而且我發現了一個小問題
就是他要求q(x,y)對稱 但是老師給的卡方自由度一分布不是對稱方程式(上面有提)
所以他又用了另外一種方法求alpha值
大家再去看看吧 英文爛說的不是很正確
: 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.185.163