Plot the "asymptotic" probability density generated by my Code
Bonjour,
Voilà je suis un cours de Monte Carlo suivant les chaînes de Markov et nous avons à répondre à ces questions :
a) Using the Metropolis algorithm, write a code that generates random numbers
according to the Gauss law : p(x) = A exp( -x^2/2). Is the constant A pertinent ?
b) Analyze the influence of , the maximum amplitude of the jumps, on the
acceptance rate and on the equilibration time. For this purpose, you may com-
pute the variance < x^2 >, that should converge to 1
VOICI MON CODE SOUS FORTRAN :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
| program Markov
real(kind=8) ksi, lambda, lambda_i, x0, x_trial,p0, p_trial,t, beta, sum_sqr_x
lambda = 0.001
print *,' Nbre de Pas ?'
read *, n_steps
x0 = 0.
p0 = exp(-x0**2/2)
sum_sqr_x = x0**2
do k = 1, n_steps
!--- Definition de la fonction p(x) ----------------------
call random_number(ksi)
lambda_i = lambda*(2*ksi-1)
x_trial = x0 + lambda_i
p_trial = exp(-x_trial**2/2)
t = min (1.0, p_trial/p0)
call random_number(beta)
if (beta < t ) then
x0 = x_trial
p0 = p_trial
end if
sum_sqr_x = sum_sqr_x + x0**2
end do
sum_sqr_x = sum_sqr_x / n_steps
print *, '<x**2 > = ', sum_sqr_x
end program MARKOV |
À partir de là, tout se passe bien ; tout marche. Mais quand il s’agit de tracer (grâce à la fonction "plot") ...the ”asymptotic” probability density... je n'y arrive pas :/ Auriez-vous des solutions ?
Merci d'avances,
Mes amitiés à mes confrères les programmeurs :)
