(* :Title: MetropolisMCMC *)

(* :Author: Mark Fisher *)

(* :Date: March 2005 *)

(* :Package Version: 1.1
	Added Gelfand-Dey functions March 2005
	Added LogLikelihood option April 2005
	Added BooleanFunction option August 2005
	Some modifications January 2006
*)

(* :Summary:
This package implements the (random-walk) Metropolis version of the Markov
Chain Monte Carlo (MCMC) algorithm. MCMC algorithms are used to draw from
probability distributions. GelfandDeyEstimate (which calls
GelfandDeyFunction) can be used to compute the normalizing constant.
*)

(* :Discussion:
The Markov Chain is implemented via FoldList. The "state" vector is the
pair {parameter vector, value} and the random "shock" is the pair {random
step vector, random test deviate}. The transition function is

Function[{state, shock},
	proposal = state[[1]] + shock[[1]]; (* random step *)
	proposalValue = fun[proposal];
	If[proposalValue >= state[[2]] * shock[[2]], (* random test *)
		(* then: accept proposal as new state *)
		{proposal, proposalValue},
		(* else: reject proposal and retain current state *)
		state]
	]

where fun is a probability density function (or a function proportional to
it). The transistion works by taking a random step and comparing the new
value with the old value scaled by a uniform(0,1) deviate.

Note the random-walk version of the Metropolis algorithm requires a symmetric
proposal distribution (such as the Gaussian distribution).
*)

(* :References:
	MetropolisMCMC implements Agorithm A.29 (p. 288) in
	Robert, Christian P. and George Casella (2004) Monte Carlo
		Statistical Methods, Second Edition, Springer.

	For the Gelfand-Dey method of marginal likelihood estimation, see
	Koop, Gary (2003) Bayesian Econometrics, John Wiley & Sons,
		pp. 104-106.
	For additional details, Koop refers the reader to
	Geweke, John (1999) Using Simulation Methods for Bayesian Econometric
		Models: Inference, Development, and Communication (with discussion
		and rejoinder), Econometric Reviews, 18, pp. 1-126.
*)

BeginPackage["MetropolisMCMC`",
	{"Statistics`NormalDistribution`"}]

MetropolisMCMC::usage = "MetropolisMCMC[fun, start, step, n] returns n
steps of a Markov Chain Monte Carlo (MCMC) simulation from the function fun
computed via the random-walk Metropolis algorithm, beginning at the vector
start of starting values, where step is a positive scalar or vector of
positive numbers (with the same dimension as start) that determines the
scale of the step size in each dimension (given the default multivariate
standard normal random steps). MetropolisMCMC takes the options
RandomStepFunction, RandomTestFunction, ReturnValues, LogLikelihood, and
BooleanFunction. (See their usage notes.)"

RandomStepFunction::usage = "RandomStepFunction is an option for
MetropolisMCMC which specifies the function to use to generate the random
steps. The specified function takes two arguments, n and d, where n is the
number of steps and d is the dimension of the parameter vector. The default
setting is RandomStepFunction -> (RandomArray[NormalDistribution[0, 1],
{##}]&)."

RandomTestFunction::usage = "RandomTestFunction is an option for
MetropolisMCMC which specifies the function to use to generate the uniform
(0,1) random deviates to use in the test for whether to accept the next
random step. The specified function takes the number of steps as its sole
argument. The default setting is RandomTestFunction -> (Table[Random[],
{#1}] &)."

ReturnValues::usage = "ReturnValues is an option for MetropolisMCMC which
specifies whether to return the values of the sampled points along with the
sampled points. The default setting is ReturnValues -> False."

LogLikelihood::usage = "LogLikelihood is an option for MetropolisMCMC that
specifies whether the function returns values proportional to the
likelihood or to the loglikelihood. The default setting is LogLikelihood ->
False. (When this option set to True, the random test deviates are
transformed appropriately.)"

BooleanFunction::usage = "BooleanFunction is an option for MetropolisMCMC
that specifies a function that returns True of False to be applied to the
proposal prior to evaluating the likelihood function. For example,
BooleanFunction -> (And@@Thread[#>=0]&). The default setting is
BooleanFunction -> (True&). The option BooleanFunction is provided for use
with LogLikelihood -> True, in order to avoid taking the log of zero."

AcceptanceRate::usage = "AcceptanceRate[mc] computes, the fraction of times
the proposal is accepted; namely, N[Length[Split[mc]]/Length[mc]]."

GelfandDeyFunction::usage = "GelfandDeyFunction[mean, cov, p] returns a
function that can be used to compute the normalization constant where mean
is the vector of means of the MCMC parameter matrix from a Markov Chain
Monte Carlo simulation and cov is the covariance matrix. The third argument
p specifies the truncation for the Gaussian in terms of a p-value.
Alternatively, one may specify GelfandDeyFunction[mcmc, p], where mcmc is
the MCMC parameter matrix of parameters."

GelfandDeyEstimate::usage = "GelfandDeyEstimate[mcmcdata, p] returns the
Gelfand-Dey estimate of the normalization constant given the Markov Chain
Monte Carlo simulation, where mcmcdata is produced by MetropolisMCMC with
the option ReturnValues -> True and where the truncation point is
controlled by the \"p-value\" p."

Begin["Private`"]

AcceptanceRate[mc_List] := N[Length[Split[mc]]/Length[mc]]


Options[MetropolisMCMC] = {
	RandomStepFunction ->
	(Statistics`Common`DistributionsCommon`RandomArray[
		Statistics`NormalDistribution`NormalDistribution[0, 1], {##}] &),
	RandomTestFunction -> (Table[Random[], {#}]&),
	ReturnValues -> False,
	LogLikelihood -> False,
	BooleanFunction -> (True&)
	}

MetropolisMCMC[fun_, start_?(VectorQ[#, NumericQ]&),
		stepSize_?(VectorQ[#, NumericQ]&), n_Integer?Positive,
		opts___?OptionQ] :=
	Module[{ranStepFun, ranTestFun, vals, loglike, boole, foldFun,
		proposal, proposalValue, shocks},
	{ranStepFun, ranTestFun, vals, loglike, boole} =
		{RandomStepFunction, RandomTestFunction, ReturnValues,
			LogLikelihood, BooleanFunction} /. {opts} /.
				Options[MetropolisMCMC];
	{loglike, vals} = TrueQ /@ {loglike, vals};
	foldFun =
		With[{op = If[loglike, Plus, Times]},
		Function[{state, shock},
			proposal = state[[1]] + stepSize * shock[[1]];
			If[boole[proposal],
				(* then *)
				proposalValue = fun[proposal];
				If[proposalValue >= state[[2]] ~ op ~ shock[[2]],
					{proposal, proposalValue},
					state],
				(* else *)
				state]
			]
		];
	shocks =
		Transpose[{
			ranStepFun[n, Length[start]],
			If[loglike, Log, Identity][ranTestFun[n]]
			}];
	FoldList[
		foldFun,
		{start, fun[start]},
		shocks
		][[
			Sequence @@ If[vals, All, {All, 1}]
			]] (* then return values with parameters, else just parameters *)
	] /; Length[start] == Length[stepSize]

MetropolisMCMC[fun_, start_?(VectorQ[#, NumericQ]&),
		stepSize_?NumericQ, n_Integer?Positive,
		opts___?OptionQ] :=
	MetropolisMCMC[fun, start, Table[stepSize, {Length[start]}], n, opts]

MetropolisMCMC[__] := $Failed

(* Gelfand-Dey *)

GelfandDeyFunction[
		mean_?(VectorQ[#, NumericQ]&),
		cov_?(MatrixQ[#, NumericQ]&),
		p_:(.01)] :=
	Module[{n, vec, halfquad},
	n = Length[mean];
	vec = Table[Unique[], {n}];
	halfquad = (mean - vec).Inverse[cov].(mean - vec)/2;
	Compile @@ {vec,
		If @@ {halfquad <= InverseGammaRegularized[n/2, 0, 1. - p],
			Exp[-halfquad]/Sqrt[(2*Pi)^n * Det[cov]]/(1. - p), 0}
		}
	]

GelfandDeyFunction[mc_?(MatrixQ[#, NumericQ]&), p_:(.01)] :=
	Module[{data, mean, centered, cov},
	data = Transpose[mc];
	mean = Mean /@ data;
	centered = data - mean;
	cov = (centered.Transpose[centered])/Length[mc];
	GelfandDeyFunction[mean, cov, p]
	]

GelfandDeyEstimate[mc_, p_:(.01)] /;
			(Dimensions[mc][[-1]] == 2 &&
			MatrixQ[mc[[All, 1]], NumericQ] &&		(* parameters *)
			VectorQ[mc[[All, 2]], NumericQ]) :=		(* values *)
	1/Mean[
		(GelfandDeyFunction[mc[[All, 1]], p] @@@ mc[[All, 1]])/mc[[All, 2]]
		]

GelfandDeyFunction[__] := $Failed

GelfandDeyEstimate[__] := $Failed

End[]
EndPackage[]