metropolis_hastings Derived Type

type, public :: metropolis_hastings

An implementation of the Metropolis-Hastings algorithm for the generation of a Markov chain. This is a default implementation that allows sampling of normally distributed posterior distributions centered on zero with unit standard deviations. Proposals are generated from a multivariate normal distribution with an identity covariance matrix and centered on zero. To alter these sampling and target distributions simply create a new class inheriting from this class and override the appropriate routines.


Contents

Type-Bound Procedures


Type-Bound Procedures

procedure, public :: compute_hastings_ratio => mh_hastings_ratio

  • private 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.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:) :: xc

    The current state vector.
    real(kind=real64), intent(in), dimension(size(xc)) :: xp

    The proposed state vector.

    Return Value real(kind=real64)

    The ratio.

procedure, public :: evaluate_proposal_pdf => mh_eval_proposal

  • private pure function mh_eval_proposal(this, xc) result(rst)

    Evaluates the proposal distribution PDF at the specified set of variables.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:) :: xc

    The array of variables to evaluate.

    Return Value real(kind=real64)

    The value of the PDF at xc.

procedure, public :: generate_proposal => mh_proposal

  • private function mh_proposal(this, xc) result(rst)

    Proposes a new sample set of variables. The sample is generated by sampling a multivariate normal distribution.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:) :: xc

    The current set of variables.

    Return Value real(kind=real64), allocatable, dimension(:)

    The proposed set of variables.

procedure, public :: get_accepted_count => mh_get_num_accepted

  • private pure function mh_get_num_accepted(this) result(rst)

    Gets the number of accepted steps.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.

    Return Value integer(kind=int32)

    The number of accepted steps.

procedure, public :: get_chain => mh_get_chain

  • private function mh_get_chain(this, bin, err) result(rst)

    Gets a copy of the stored Markov chain.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), optional :: bin

    An optional input allowing for a burn-in region. The parameter represents the amount (percentage-based) of the overall chain to disregard as "burn-in" values. The value shoud exist on [0, 1). The default value is 0 such that no values are disregarded.
    class(errors), intent(inout), optional, target :: err

    The error handling object.

    Return Value real(kind=real64), allocatable, dimension(:,:)

    The resulting chain with each parameter represented by a column.

procedure, public :: get_chain_length => mh_get_chain_length

  • private pure function mh_get_chain_length(this) result(rst)

    Gets the length of the chain (number of stored state variables).

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.

    Return Value integer(kind=int32)

    The chain length.

procedure, public :: get_proposal_cholesky => mh_get_prop_chol_cov

  • private pure function mh_get_prop_chol_cov(this) result(rst)

    Gets the Cholesky-factored (lower-triangular) form of the proposal covariance matrix.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.

    Return Value real(kind=real64), allocatable, dimension(:,:)

    The Cholesky-factored form of the proposal covariance matrix store in lower-triangular form.

procedure, public :: get_proposal_covariance => mh_get_prop_cov

  • private pure function mh_get_prop_cov(this) result(rst)

    Gets the covariance matrix of the proposal distribution.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.

    Return Value real(kind=real64), allocatable, dimension(:,:)

    The covariance matrix.

procedure, public :: get_proposal_initialized => mh_get_is_prop_init

  • private pure function mh_get_is_prop_init(this) result(rst)

    Gets a value determining if the proposal distribution object has been initialized.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.

    Return Value logical

    Returns true if the object has been initialized; else, false.

procedure, public :: get_proposal_means => mh_get_prop_mean

  • private pure function mh_get_prop_mean(this) result(rst)

    Gets the mean values of the proposal distribution.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.

    Return Value real(kind=real64), allocatable, dimension(:)

    An array containing the mean values.

procedure, public :: get_state_variable_count => mh_get_nvars

  • private pure function mh_get_nvars(this) result(rst)

    Gets the number of state variables.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(in) :: this

    The metropolis_hastings object.

    Return Value integer(kind=int32)

    The number of state variables.

generic, public :: initialize_proposal => mh_init_proposal_1, mh_init_proposal_2

  • private subroutine mh_init_proposal_1(this, mu, sigma, err)

    Initializes the multivariate normal distribution used to generate proposals.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:) :: mu

    An N-element array containing the mean values for the distribution.
    real(kind=real64), intent(in), dimension(:,:) :: sigma

    An N-by-N covariance matrix for the distribution. This matrix must be positive-definite.
    class(errors), intent(inout), optional, target :: err

    An error handling object.

  • private subroutine mh_init_proposal_2(this, n, err)

    Initializes the multivariate normal distribution to a mean of zero and a variance of one.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    integer(kind=int32), intent(in) :: n

    The number of state variables.
    class(errors), intent(inout), optional, target :: err

    An error handling object.

procedure, public :: on_acceptance => mh_on_success

  • private subroutine mh_on_success(this, iter, alpha, xc, xp, err)

    Currently, this routine does nothing and is a placeholder for the user that inherits this class to provide functionallity upon acceptance of a proposed value.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    integer(kind=int32), intent(in) :: iter

    The current iteration number.
    real(kind=real64), intent(in) :: alpha

    The proposal probabilty term used for acceptance criteria.
    real(kind=real64), intent(in), dimension(:) :: xc

    An N-element array containing the current state variables.
    real(kind=real64), intent(in), dimension(size(xc)) :: xp

    An N-element array containing the proposed state variables that were just accepted.
    class(errors), intent(inout), optional, target :: err

    An error handling object.

procedure, public :: on_rejection => mh_on_rejection

  • private subroutine mh_on_rejection(this, iter, alpha, xc, xp, err)

    Currently, this routine does nothing and is a placeholder for the user that inherits this class to provide functionallity upon rejection of a proposed value.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    integer(kind=int32), intent(in) :: iter

    The current iteration number.
    real(kind=real64), intent(in) :: alpha

    The proposal probabilty term used for acceptance criteria.
    real(kind=real64), intent(in), dimension(:) :: xc

    An N-element array containing the current state variables.
    real(kind=real64), intent(in), dimension(size(xc)) :: xp

    An N-element array containing the proposed state variables that were just rejected.
    class(errors), intent(inout), optional, target :: err

    An error handling object.

procedure, public :: push_new_state => mh_push

  • private subroutine mh_push(this, x, err)

    Pushes a new set of state variables onto the buffer.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:) :: x

    The new N-element state array.
    class(errors), intent(inout), optional, target :: err

    The error handling object.

procedure, public :: reset => mh_clear_chain

  • private subroutine mh_clear_chain(this)

    Resets the object and clears out the buffer storing the chain values.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.

procedure, public :: sample => mh_sample

  • private subroutine mh_sample(this, xi, niter, err)

    Samples the distribution using the Metropolis-Hastings algorithm.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:) :: xi

    An N-element array containing initial starting values of the state variables.
    integer(kind=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.

procedure, public :: set_proposal_covariance => mh_set_prop_cov

  • private subroutine mh_set_prop_cov(this, x, err)

    Sets the covariance matrix of the proposal distribution.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:,:) :: x

    The covariance matrix. This matrix must be positive-definite.
    class(errors), intent(inout), optional, target :: err

    An error handling object.

procedure, public :: set_proposal_means => mh_set_prop_mean

  • private subroutine mh_set_prop_mean(this, x, err)

    Sets the mean values of the proposal distribution.

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:) :: x

    The updated mean values.
    class(errors), intent(inout), optional, target :: err

    An error handling object.

procedure, public :: target_distribution => mh_target

  • private 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 sample a multivariate normal distribution with a mean of zero and a variance of one (identity covariance matrix).

    Arguments

    Type IntentOptional Attributes Name
    class(metropolis_hastings), intent(inout) :: this

    The metropolis_hastings object.
    real(kind=real64), intent(in), dimension(:) :: x

    The state vector.

    Return Value real(kind=real64)

    The value of the probability density function of the distribution being sampled.