Skip to content

Commit e48403f

Browse files
jchristophersonjchristopherson
committed
Clean up
1 parent 4b0973e commit e48403f

File tree

4 files changed

+178
-537
lines changed

4 files changed

+178
-537
lines changed

‎examples/mcmc_example.f90‎

Lines changed: 126 additions & 126 deletions
Original file line numberDiff line numberDiff line change
@@ -23,130 +23,130 @@ program example
2323
implicit none
2424

2525
! Local Variables
26-
logical :: check
27-
integer(int32) :: i, n
28-
procedure(regression_function), pointer :: fcn
29-
real(real64) :: xp(21), yp(21), mdl(4), ym(21)
30-
real(real64), allocatable, dimension(:,:) :: chain
31-
type(metropolis_hastings) :: mcmc
32-
33-
! Plot Variables
34-
type(multiplot) :: plt, pplt
35-
class(terminal), pointer :: term
36-
type(plot_2d) :: plt1, plt2, plt3, plt4, xyplt
37-
type(plot_data_histogram) :: pdh
38-
type(plot_data_2d) :: pd
39-
40-
! Initialization
41-
fcn => exfun
42-
43-
! Data to fit
44-
xp = [0.0d0, 0.1d0, 0.2d0, 0.3d0, 0.4d0, 0.5d0, 0.6d0, 0.7d0, 0.8d0, &
45-
0.9d0, 1.0d0, 1.1d0, 1.2d0, 1.3d0, 1.4d0, 1.5d0, 1.6d0, 1.7d0, &
46-
1.8d0, 1.9d0, 2.0d0]
47-
yp = [1.216737514d0, 1.250032542d0, 1.305579195d0, 1.040182335d0, &
48-
1.751867738d0, 1.109716707d0, 2.018141531d0, 1.992418729d0, &
49-
1.807916923d0, 2.078806005d0, 2.698801324d0, 2.644662712d0, &
50-
3.412756702d0, 4.406137221d0, 4.567156645d0, 4.999550779d0, &
51-
5.652854194d0, 6.784320119d0, 8.307936836d0, 8.395126494d0, &
52-
10.30252404d0]
53-
54-
! Generate an initial estimate - based upon prior knowledge of the problem
55-
mdl = [1.186d0, 0.447d0, -0.12d0, 1.06d0]
56-
call random_number(mdl)
57-
58-
! Initialize the MH object
59-
call mcmc%initialize(fcn, xp, yp, mdl)
60-
61-
! Form the Markov chain
62-
call mcmc%evaluate(fcn, xp, yp)
63-
64-
! Extract the chain and plot
65-
chain = mcmc%get_chain()
66-
n = mcmc%get_chain_length()
67-
print "(AI0)", "Chain Length: ", n
68-
69-
! Update the model using the means of each parameter
70-
mdl = chain(n,:)
71-
call fcn(xp, mdl, ym, check)
72-
73-
! Report out the results
74-
do i = 1, size(mdl)
75-
print "(AI0)", "Parameter ", i
76-
print "(AAF10.7)", achar(9), "Value: ", chain(n, i)
77-
print "(AAF10.7)", achar(9), "Mean: ", mean(chain(:,i))
78-
print "(AAF10.7)", achar(9), "Std. Dev.: ", standard_deviation(chain(:,i))
79-
end do
80-
81-
! Plot the histograms for each parameter
82-
call plt%initialize(2, 2)
83-
term => plt%get_terminal()
84-
call term%set_window_height(800)
85-
call term%set_window_width(1200)
86-
call plt1%initialize()
87-
call plt2%initialize()
88-
call plt3%initialize()
89-
call plt4%initialize()
90-
91-
call plt1%set_title("p_1")
92-
call pdh%define_data(chain(:,1))
93-
call pdh%set_transparency(0.5)
94-
call plt1%push(pdh)
95-
96-
call plt2%set_title("p_2")
97-
call pdh%define_data(chain(:,2))
98-
call plt2%push(pdh)
99-
100-
call plt3%set_title("p_3")
101-
call pdh%define_data(chain(:,3))
102-
call plt3%push(pdh)
103-
104-
call plt4%set_title("p_4")
105-
call pdh%define_data(chain(:,4))
106-
call plt4%push(pdh)
107-
108-
call plt%set(1, 1, plt1)
109-
call plt%set(2, 1, plt2)
110-
call plt%set(1, 2, plt3)
111-
call plt%set(2, 2, plt4)
112-
call plt%draw()
113-
114-
! Generate a parameter plot
115-
call pplt%initialize(2, 2)
116-
term => pplt%get_terminal()
117-
call term%set_window_height(800)
118-
call term%set_window_width(1200)
119-
call plt1%clear_all()
120-
call plt2%clear_all()
121-
call plt3%clear_all()
122-
call plt4%clear_all()
123-
124-
call pd%define_data(chain(:,1))
125-
call plt1%push(pd)
126-
127-
call pd%define_data(chain(:,2))
128-
call plt2%push(pd)
129-
130-
call pd%define_data(chain(:,3))
131-
call plt3%push(pd)
132-
133-
call pd%define_data(chain(:,4))
134-
call plt4%push(pd)
135-
136-
call pplt%set(1, 1, plt1)
137-
call pplt%set(2, 1, plt2)
138-
call pplt%set(1, 2, plt3)
139-
call pplt%set(2, 2, plt4)
140-
call pplt%draw()
141-
142-
! Generate an X-Y plot based upon the means of each parameter set
143-
call xyplt%initialize()
144-
call pd%define_data(xp, ym)
145-
call pd%set_line_width(2.0)
146-
call xyplt%push(pd)
147-
call pd%define_data(xp, yp)
148-
call pd%set_draw_line(.false.)
149-
call pd%set_draw_markers(.true.)
150-
call xyplt%push(pd)
151-
call xyplt%draw()
26+
! logical :: check
27+
! integer(int32) :: i, n
28+
! procedure(regression_function), pointer :: fcn
29+
! real(real64) :: xp(21), yp(21), mdl(4), ym(21)
30+
! real(real64), allocatable, dimension(:,:) :: chain
31+
! type(metropolis_hastings) :: mcmc
32+
33+
! ! Plot Variables
34+
! type(multiplot) :: plt, pplt
35+
! class(terminal), pointer :: term
36+
! type(plot_2d) :: plt1, plt2, plt3, plt4, xyplt
37+
! type(plot_data_histogram) :: pdh
38+
! type(plot_data_2d) :: pd
39+
40+
! ! Initialization
41+
! fcn => exfun
42+
43+
! ! Data to fit
44+
! xp = [0.0d0, 0.1d0, 0.2d0, 0.3d0, 0.4d0, 0.5d0, 0.6d0, 0.7d0, 0.8d0, &
45+
! 0.9d0, 1.0d0, 1.1d0, 1.2d0, 1.3d0, 1.4d0, 1.5d0, 1.6d0, 1.7d0, &
46+
! 1.8d0, 1.9d0, 2.0d0]
47+
! yp = [1.216737514d0, 1.250032542d0, 1.305579195d0, 1.040182335d0, &
48+
! 1.751867738d0, 1.109716707d0, 2.018141531d0, 1.992418729d0, &
49+
! 1.807916923d0, 2.078806005d0, 2.698801324d0, 2.644662712d0, &
50+
! 3.412756702d0, 4.406137221d0, 4.567156645d0, 4.999550779d0, &
51+
! 5.652854194d0, 6.784320119d0, 8.307936836d0, 8.395126494d0, &
52+
! 10.30252404d0]
53+
54+
! ! Generate an initial estimate - based upon prior knowledge of the problem
55+
! mdl = [1.186d0, 0.447d0, -0.12d0, 1.06d0]
56+
! call random_number(mdl)
57+
58+
! ! Initialize the MH object
59+
! call mcmc%initialize(fcn, xp, yp, mdl)
60+
61+
! ! Form the Markov chain
62+
! call mcmc%evaluate(fcn, xp, yp)
63+
64+
! ! Extract the chain and plot
65+
! chain = mcmc%get_chain()
66+
! n = mcmc%get_chain_length()
67+
! print "(AI0)", "Chain Length: ", n
68+
69+
! ! Update the model using the means of each parameter
70+
! mdl = chain(n,:)
71+
! call fcn(xp, mdl, ym, check)
72+
73+
! ! Report out the results
74+
! do i = 1, size(mdl)
75+
! print "(AI0)", "Parameter ", i
76+
! print "(AAF10.7)", achar(9), "Value: ", chain(n, i)
77+
! print "(AAF10.7)", achar(9), "Mean: ", mean(chain(:,i))
78+
! print "(AAF10.7)", achar(9), "Std. Dev.: ", standard_deviation(chain(:,i))
79+
! end do
80+
81+
! ! Plot the histograms for each parameter
82+
! call plt%initialize(2, 2)
83+
! term => plt%get_terminal()
84+
! call term%set_window_height(800)
85+
! call term%set_window_width(1200)
86+
! call plt1%initialize()
87+
! call plt2%initialize()
88+
! call plt3%initialize()
89+
! call plt4%initialize()
90+
91+
! call plt1%set_title("p_1")
92+
! call pdh%define_data(chain(:,1))
93+
! call pdh%set_transparency(0.5)
94+
! call plt1%push(pdh)
95+
96+
! call plt2%set_title("p_2")
97+
! call pdh%define_data(chain(:,2))
98+
! call plt2%push(pdh)
99+
100+
! call plt3%set_title("p_3")
101+
! call pdh%define_data(chain(:,3))
102+
! call plt3%push(pdh)
103+
104+
! call plt4%set_title("p_4")
105+
! call pdh%define_data(chain(:,4))
106+
! call plt4%push(pdh)
107+
108+
! call plt%set(1, 1, plt1)
109+
! call plt%set(2, 1, plt2)
110+
! call plt%set(1, 2, plt3)
111+
! call plt%set(2, 2, plt4)
112+
! call plt%draw()
113+
114+
! ! Generate a parameter plot
115+
! call pplt%initialize(2, 2)
116+
! term => pplt%get_terminal()
117+
! call term%set_window_height(800)
118+
! call term%set_window_width(1200)
119+
! call plt1%clear_all()
120+
! call plt2%clear_all()
121+
! call plt3%clear_all()
122+
! call plt4%clear_all()
123+
124+
! call pd%define_data(chain(:,1))
125+
! call plt1%push(pd)
126+
127+
! call pd%define_data(chain(:,2))
128+
! call plt2%push(pd)
129+
130+
! call pd%define_data(chain(:,3))
131+
! call plt3%push(pd)
132+
133+
! call pd%define_data(chain(:,4))
134+
! call plt4%push(pd)
135+
136+
! call pplt%set(1, 1, plt1)
137+
! call pplt%set(2, 1, plt2)
138+
! call pplt%set(1, 2, plt3)
139+
! call pplt%set(2, 2, plt4)
140+
! call pplt%draw()
141+
142+
! ! Generate an X-Y plot based upon the means of each parameter set
143+
! call xyplt%initialize()
144+
! call pd%define_data(xp, ym)
145+
! call pd%set_line_width(2.0)
146+
! call xyplt%push(pd)
147+
! call pd%define_data(xp, yp)
148+
! call pd%set_draw_line(.false.)
149+
! call pd%set_draw_markers(.true.)
150+
! call xyplt%push(pd)
151+
! call xyplt%draw()
152152
end program

0 commit comments

Comments
 (0)