Skip to content

Commit

Permalink
Update
Browse files Browse the repository at this point in the history
  • Loading branch information
@jchristopherson
jchristopherson committed Feb 14, 2025
1 parent e48403f commit 90cef53
Show file tree
Hide file tree
Showing 4 changed files with 223 additions and 146 deletions.
170 changes: 26 additions & 144 deletions examples/mcmc_example.f90
View file
Original file line number Diff line number Diff line change
@@ -1,152 +1,34 @@
module nl_example
use iso_fortran_env
contains
subroutine exfun(x, p, f, stop)
! Arguments
real(real64), intent(in) :: x(:), p(:)
real(real64), intent(out) :: f(:)
logical, intent(out) :: stop

! Function
f = p(4) * x**3 + p(3) * x**2 + p(2) * x + p(1)

! No need to stop
stop = .false.
end subroutine
end module

program example
use iso_fortran_env
use fstats
use fplot_core
use nl_example
implicit none

! Local Variables
! logical :: check
! integer(int32) :: i, n
! procedure(regression_function), pointer :: fcn
! real(real64) :: xp(21), yp(21), mdl(4), ym(21)
! real(real64), allocatable, dimension(:,:) :: chain
! type(metropolis_hastings) :: mcmc

! ! Plot Variables
! type(multiplot) :: plt, pplt
! class(terminal), pointer :: term
! type(plot_2d) :: plt1, plt2, plt3, plt4, xyplt
! type(plot_data_histogram) :: pdh
! type(plot_data_2d) :: pd

! ! Initialization
! fcn => exfun

! ! Data to fit
! xp = [0.0d0, 0.1d0, 0.2d0, 0.3d0, 0.4d0, 0.5d0, 0.6d0, 0.7d0, 0.8d0, &
! 0.9d0, 1.0d0, 1.1d0, 1.2d0, 1.3d0, 1.4d0, 1.5d0, 1.6d0, 1.7d0, &
! 1.8d0, 1.9d0, 2.0d0]
! yp = [1.216737514d0, 1.250032542d0, 1.305579195d0, 1.040182335d0, &
! 1.751867738d0, 1.109716707d0, 2.018141531d0, 1.992418729d0, &
! 1.807916923d0, 2.078806005d0, 2.698801324d0, 2.644662712d0, &
! 3.412756702d0, 4.406137221d0, 4.567156645d0, 4.999550779d0, &
! 5.652854194d0, 6.784320119d0, 8.307936836d0, 8.395126494d0, &
! 10.30252404d0]

! ! Generate an initial estimate - based upon prior knowledge of the problem
! mdl = [1.186d0, 0.447d0, -0.12d0, 1.06d0]
! call random_number(mdl)

! ! Initialize the MH object
! call mcmc%initialize(fcn, xp, yp, mdl)

! ! Form the Markov chain
! call mcmc%evaluate(fcn, xp, yp)

! ! Extract the chain and plot
! chain = mcmc%get_chain()
! n = mcmc%get_chain_length()
! print "(AI0)", "Chain Length: ", n

! ! Update the model using the means of each parameter
! mdl = chain(n,:)
! call fcn(xp, mdl, ym, check)

! ! Report out the results
! do i = 1, size(mdl)
! print "(AI0)", "Parameter ", i
! print "(AAF10.7)", achar(9), "Value: ", chain(n, i)
! print "(AAF10.7)", achar(9), "Mean: ", mean(chain(:,i))
! print "(AAF10.7)", achar(9), "Std. Dev.: ", standard_deviation(chain(:,i))
! end do

! ! Plot the histograms for each parameter
! call plt%initialize(2, 2)
! term => plt%get_terminal()
! call term%set_window_height(800)
! call term%set_window_width(1200)
! call plt1%initialize()
! call plt2%initialize()
! call plt3%initialize()
! call plt4%initialize()

! call plt1%set_title("p_1")
! call pdh%define_data(chain(:,1))
! call pdh%set_transparency(0.5)
! call plt1%push(pdh)

! call plt2%set_title("p_2")
! call pdh%define_data(chain(:,2))
! call plt2%push(pdh)

! call plt3%set_title("p_3")
! call pdh%define_data(chain(:,3))
! call plt3%push(pdh)

! call plt4%set_title("p_4")
! call pdh%define_data(chain(:,4))
! call plt4%push(pdh)

! call plt%set(1, 1, plt1)
! call plt%set(2, 1, plt2)
! call plt%set(1, 2, plt3)
! call plt%set(2, 2, plt4)
! call plt%draw()

! ! Generate a parameter plot
! call pplt%initialize(2, 2)
! term => pplt%get_terminal()
! call term%set_window_height(800)
! call term%set_window_width(1200)
! call plt1%clear_all()
! call plt2%clear_all()
! call plt3%clear_all()
! call plt4%clear_all()

! call pd%define_data(chain(:,1))
! call plt1%push(pd)

! call pd%define_data(chain(:,2))
! call plt2%push(pd)

! call pd%define_data(chain(:,3))
! call plt3%push(pd)

! call pd%define_data(chain(:,4))
! call plt4%push(pd)

! call pplt%set(1, 1, plt1)
! call pplt%set(2, 1, plt2)
! call pplt%set(1, 2, plt3)
! call pplt%set(2, 2, plt4)
! call pplt%draw()

! ! Generate an X-Y plot based upon the means of each parameter set
! call xyplt%initialize()
! call pd%define_data(xp, ym)
! call pd%set_line_width(2.0)
! call xyplt%push(pd)
! call pd%define_data(xp, yp)
! call pd%set_draw_line(.false.)
! call pd%set_draw_markers(.true.)
! call xyplt%push(pd)
! call xyplt%draw()
real(real64) :: xi(2)
real(real64), allocatable, dimension(:,:) :: chains
type(metropolis_hastings) :: mcmc

! Plot Variables
type(plot_2d) :: plt
type(plot_data_2d) :: pd

! Create an initial estimate
call random_number(xi)

! Sample a multivariate normal distribution using MH
call mcmc%sample(xi)

! Get the chains
chains = mcmc%get_chain()

! Plot the results
call plt%initialize()
call pd%define_data(chains(:,1), chains(:,2))
call pd%set_draw_line(.false.)
call pd%set_draw_markers(.true.)
call pd%set_marker_style(MARKER_FILLED_CIRCLE)
call pd%set_marker_scaling(0.5)
call plt%push(pd)
call plt%draw()
end program
8 changes: 6 additions & 2 deletions src/fstats.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,8 +76,7 @@ module fstats
public :: scaled_random_resample
public :: bootstrap_statistics
public :: bootstrap
public :: box_muller_sample
public :: rejection_sample

public :: lowess
public :: pooled_variance
public :: bartletts_test
Expand All @@ -90,6 +89,11 @@ module fstats
public :: doe_evaluate_model
public :: doe_model

! FSTATS_SAMPLING.F90
public :: box_muller_sample
public :: rejection_sample
public :: sample_normal_multivariate

! FSTATS_MCMC.F90
public :: metropolis_hastings

Expand Down
162 changes: 162 additions & 0 deletions src/fstats_mcmc.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ module fstats_mcmc
use fstats_descriptive_statistics
use fstats_types
use fstats_regression
use fstats_sampling
implicit none
private
public :: metropolis_hastings
Expand All @@ -22,11 +23,17 @@ module fstats_mcmc
!! The buffer where each new state is stored as a row in the matrix.
integer(int32), private :: m_numVars = 0
!! The number of state variables.
type(multivariate_normal_distribution) :: m_propDist
!! The proposal distribution from which to draw samples.
contains
procedure, public :: get_state_variable_count => mh_get_nvars
procedure, public :: get_chain_length => mh_get_chain_length
procedure, public :: push_new_state => mh_push
procedure, public :: get_chain => mh_get_chain
procedure, public :: generate_proposal => mh_proposal
procedure, public :: compute_hastings_ratio => mh_hastings_ratio
procedure, public :: target_distribution => mh_target
procedure, public :: sample => mh_sample
procedure, private :: resize_buffer => mh_resize_buffer
procedure, private :: get_buffer_length => mh_get_buffer_length
end type
Expand Down Expand Up @@ -203,5 +210,160 @@ function mh_get_chain(this, err) result(rst)
end if
end function

! ------------------------------------------------------------------------------
function mh_proposal(this, xc) result(rst)
!! Proposes a new sample set of variables. The sample is generated by
!! sampling a multivariate normal distribution.
class(metropolis_hastings), intent(inout) :: this
!! The metropolis_hastings object.
real(real64), intent(in), dimension(:) :: xc
!! The current set of variables.
real(real64), allocatable, dimension(:) :: rst
!! The proposed set of variables.

! Sample from the distribution
call this%m_propDist%set_means(xc) ! center the distribution at the current state
rst = sample_normal_multivariate(this%m_propDist)
end function

! ------------------------------------------------------------------------------
function mh_hastings_ratio(this, xc, xp) result(rst)
!! Evaluates the Hasting's ratio. If the proposal distribution is
!! symmetric, this ratio is unity; however, in the case of an asymmetric
!! distribution this ratio is not ensured to be unity.
class(metropolis_hastings), intent(inout) :: this
!! The metropolis_hastings object.
real(real64), intent(in), dimension(:) :: xc
!! The current state vector.
real(real64), intent(in), dimension(size(xc)) :: xp
!! The proposed state vector.
real(real64) :: rst
!! The ratio.

! Local Variables
real(real64), allocatable, dimension(:) :: means
real(real64) :: q1, q2

! Initialization
means = this%m_propDist%get_means()

! Process
call this%m_propDist%set_means(xc)
q1 = this%m_propDist%pdf(xp)

call this%m_propDist%set_means(xp)
q2 = this%m_propDist%pdf(xc)

! Restore the state of the object
call this%m_propDist%set_means(means)

! Compute the ratio
rst = q2 / q1
end function

! ------------------------------------------------------------------------------
function mh_target(this, x) result(rst)
!! Returns the probability value from the target distribution at the
!! specified state. The user is expected to overload this routine to
!! define the desired distribution. The default behavior of this
!! routine is to sammple a multivariate normal distribution with a mean
!! of zero and a variance of one (identity covariance matrix).
class(metropolis_hastings), intent(in) :: this
!! The metropolis_hastings object.
real(real64), intent(in), dimension(:) :: x
!! The state vector.
real(real64) :: rst
!! The value of the probability density function of the distribution
!! being sampled.

! Local Variables
integer(int32) :: i, n
real(real64), allocatable, dimension(:) :: mu
real(real64), allocatable, dimension(:,:) :: sigma
type(multivariate_normal_distribution) :: dist

! The default will be a multivariate normal distribution with an identity
! matrix for the covariance matrix
n = size(x)
allocate(mu(n), sigma(n, n), source = 0.0d0)
do i = 1, n
sigma(i,i) = 1.0d0
end do
call dist%initialize(mu, sigma)
rst = dist%pdf(x)
end function

! ------------------------------------------------------------------------------
subroutine mh_sample(this, xi, niter, err)
!! Samples the distribution using the Metropolis-Hastings algorithm.
class(metropolis_hastings), intent(inout) :: this
!! The metropolis_hastings object.
real(real64), intent(in), dimension(:) :: xi
!! An initial estimate of the state variables.
integer(int32), intent(in), optional :: niter
!! An optional input defining the number of iterations to take. The
!! default is 10,000.
class(errors), intent(inout), optional, target :: err
!! The error handling object.

! Local Variables
integer(int32) :: i, npts
real(real64) :: r, pp, pc, a, a1, a2, alpha
real(real64), allocatable, dimension(:) :: xc, xp
class(errors), pointer :: errmgr
type(errors), target :: deferr

! Initialization
if (present(err)) then
errmgr => err
else
errmgr => deferr
end if
if (present(niter)) then
npts = niter
else
npts = this%initial_iteration_estimate
end if

! TO DO: Reset the buffer state

! Initialize the stored distribution
! TO DO: Figure out a means for generating a meaningful covariance matrix.

! Store the initial value
call this%push_new_state(xi, err = errmgr)
if (errmgr%has_error_occurred()) return

! Iteration Process
xc = xi
pc = this%target_distribution(xc)
do i = 2, npts
! Create a proposal & evaluate it's PDF
xp = this%generate_proposal(xc)
pp = this%target_distribution(xp)

! Evaluate the probabilities
a1 = this%compute_hastings_ratio(xc, xp)
a2 = pp / pc
alpha = min(1.0d0, a1 * a2)

! Do we keep our current state or move to the new state?
call random_number(r)
if (r <= alpha) then
! Take the new value
call this%push_new_state(xp, err = errmgr)
if (errmgr%has_error_occurred()) return

! Update the values
xc = xp
pc = pp
else
! Keep our current estimate
call this%push_new_state(xc, err = errmgr)
if (errmgr%has_error_occurred()) return
end if
end do
end subroutine

! ------------------------------------------------------------------------------
end module
Loading

0 comments on commit 90cef53

Please sign in to comment.