|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <time.h>
|
|
#include <math.h>
|
|
|
|
#define NSAMPLE 10000
|
|
#define N 10
|
|
|
|
void init_prob(double *p);
|
|
void metropolis(double *p, int *samples);
|
|
void frequency(int *samples, double *freq);
|
|
|
|
int main() {
|
|
double p[N];
|
|
double freq[N] = {0};
|
|
int samples[NSAMPLE];
|
|
|
|
srand48(time(NULL));
|
|
|
|
init_prob(p);
|
|
metropolis(p, samples);
|
|
frequency(samples, freq);
|
|
|
|
FILE *out = fopen("frequency.dat", "w");
|
|
int i;
|
|
for(i = 0; i < N; i++) {
|
|
fprintf(out, "%d %.4lf %.4lf\n", i, p[i], freq[i]);
|
|
}
|
|
fclose(out);
|
|
|
|
return 0;
|
|
}
|
|
|
|
void init_prob(double *p) {
|
|
int i;
|
|
double norm = 0.;
|
|
for(i = 0; i < N; i++) {
|
|
p[i] = drand48();
|
|
norm += p[i];
|
|
}
|
|
for(i = 0; i < N; i++) {
|
|
p[i] /= norm;
|
|
}
|
|
}
|
|
|
|
void metropolis(double *p, int *samples) {
|
|
int n_init = drand48() * N;
|
|
int i, j;
|
|
int n_prev = n_init;
|
|
for(i = 0; i < NSAMPLE; i++) {
|
|
int n_tmp = drand48() * N;
|
|
double p_tmp = p[n_tmp];
|
|
double p_prev = p[n_prev];
|
|
double A = p_tmp / p_prev;
|
|
if(A > 1 || drand48() <= A) {
|
|
n_prev = n_tmp;
|
|
}
|
|
samples[i] = n_prev;
|
|
}
|
|
}
|
|
|
|
void frequency(int *samples, double *freq) {
|
|
int i;
|
|
for(i = 0; i < NSAMPLE; i++) {
|
|
freq[samples[i]]++;
|
|
}
|
|
for(i = 0; i < N; i++) {
|
|
freq[i] /= NSAMPLE;
|
|
}
|
|
}
|