From 90cef537e19b9a9adc6a35f4ecb42d3f9d883158 Mon Sep 17 00:00:00 2001 From: Jason Christopherson Date: Fri, 14 Feb 2025 13:56:22 -0600 Subject: [PATCH] Update --- examples/mcmc_example.f90 | 170 ++++++-------------------------------- src/fstats.f90 | 8 +- src/fstats_mcmc.f90 | 162 ++++++++++++++++++++++++++++++++++++ src/fstats_sampling.f90 | 29 +++++++ 4 files changed, 223 insertions(+), 146 deletions(-) diff --git a/examples/mcmc_example.f90 b/examples/mcmc_example.f90 index 95c2493..59fd86e 100644 --- a/examples/mcmc_example.f90 +++ b/examples/mcmc_example.f90 @@ -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 \ No newline at end of file diff --git a/src/fstats.f90 b/src/fstats.f90 index 3e45a94..c3d9d4f 100644 --- a/src/fstats.f90 +++ b/src/fstats.f90 @@ -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 @@ -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 diff --git a/src/fstats_mcmc.f90 b/src/fstats_mcmc.f90 index f1915b1..18754ad9 100644 --- a/src/fstats_mcmc.f90 +++ b/src/fstats_mcmc.f90 @@ -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 @@ -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 @@ -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 \ No newline at end of file diff --git a/src/fstats_sampling.f90 b/src/fstats_sampling.f90 index 4d47bbc..42bcdeb 100644 --- a/src/fstats_sampling.f90 +++ b/src/fstats_sampling.f90 @@ -6,6 +6,7 @@ module fstats_sampling private public :: box_muller_sample public :: rejection_sample + public :: sample_normal_multivariate real(real64), parameter :: pi = 2.0d0 * acos(0.0d0) real(real64), parameter :: twopi = 2.0d0 * pi @@ -131,5 +132,33 @@ function rejection_sample(tdist, n, xmin, xmax) result(rst) end do end function +! ****************************************************************************** +! MULTIVARIATE SAMPLING +! ------------------------------------------------------------------------------ +function sample_normal_multivariate(dist) result(rst) + !! Samples a multivariate normal distribution such that \(\vec{x} = + !! \vec{mu} + L \vec{u}\), where \(L\) is the lower form of the Cholesky + !! factorization of the covariance matrix, and \(\vec{u}\) is a randomly + !! generated vector that exists on the set \([0 1]\) + class(multivariate_normal_distribution), intent(in) :: dist + !! The multivariate normal distribution to sample. + real(real64), allocatable, dimension(:) :: rst + !! The resulting vector. + + ! Local Variables + integer(int32) :: n + real(real64), allocatable, dimension(:) :: u + real(real64), allocatable, dimension(:,:) :: L + + ! Initialization + L = dist%get_cholesky_factored_matrix() + n = size(L, 1) + allocate(u(n)) + call random_number(u) ! populating u from [0, 1]. + + ! Process + rst = dist%get_means() + matmul(L, u) +end function + ! ------------------------------------------------------------------------------ end module \ No newline at end of file