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 :
À 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 ?
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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
Merci d'avances,
Mes amitiés à mes confrères les programmeurs
Partager