/*取一個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