作者Masaki2005 (我是建威XD)
看板NCTU-STAT94G
标题[心得] 第一题~~大家参考XD
时间Thu Aug 11 14:28:42 2005
/*取一个Uniform(0,100)作为参考的函数(jumping distribution)
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "myfunction.h"
#define PI 3.14159265358979323846
#define autoseed(x) ((rand()/32767.)*x+x)
#define Chisq(x) exp(-x/2.)/sqrt(2*PI*x)
#define Min(x,y) (x<y)?x:y
FILE *out;
main()
{
int i=0,j=0;
double Seed=1;
double ThetaStar,ThetaT=1,Alpha;//起始位置=1
double U;
out=fopen("Result.txt","w");
do
{
ThetaStar=unif(autoseed(Seed))*100;
//printf("ThetaStar = %lf ",ThetaStar);
Seed++;
//printf("ThetaT = %lf ",ThetaT);
//printf("*= %lf ",Chisq(ThetaStar));
//printf("T= %lf ",Chisq(ThetaT));
Alpha=min((Chisq(ThetaStar)/Chisq(ThetaT)),1);
//printf("Alpha = %lf ",Alpha);
U=unif(autoseed(Seed));
//printf("U = %lf\n",U);
Seed++;
if(Alpha>=U)
{
ThetaT=ThetaStar;
if(j>=500000)
{
printf("Theta[%d] = %lf\t",i+1,ThetaT);
fprintf(out,"%lf\n",ThetaT);
i++;
}
}
j++;
}while(i!=1000);
fclose(out);
}
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.113.7.249