作者lalalarc (初来乍到)
看板NCTU-STAT94G
标题Re: [闲聊] MCMC
时间Thu Aug 11 10:45:56 2005
※ 引述《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