added code and organise folders

code/mcmc.m

+(* ::Package:: *)
+
+(* Mathematica Package *)
+
+BeginPackage["MCMC`", {"Units`","DataFits`","BBHReduce`","NRTimeSeries`"}];
+
+MCMC::usage = "MCMC[plogexpr, paramspec, numsteps]
+
+Perform MCMC sampling of the supplied probability distribution.
+
+1. plogexpr should be an expression that gives the unnormalized log
+probability for a particular choice of parameter values.
+
+2. paramspec either gives the results of a previous MCMC run (w/ same
+plogexpr--just to add on more iterations), or lists the model parameters
+like so:
+{{param1, ival1, spread1, domain1}, ...}
+a) Each param should be symbolic.
+b) ival is the initial parameter value.
+c) spread is roughly how far to try to change the parameter each step in
+the Markov chain. In this routine we select new parameters values based
+on an exponential distribution of the form Exp[\[CapitalDelta]param/spread]. My
+Numerical Recipes book advises setting these spreads so that the average
+candidate acceptance is 10-40%.
+d) Each domain is either Reals or a list of all possible values the
+parameter can take on (needs to be a uniform grid).
+
+3. numsteps is the number of Markov chain steps to perform.";
+
+MCMCModelFit::usage = "MCMCModelFit[data, errors, model, paramspec, ivars, numsteps]
+
+Perform MCMC samping of the probability distribution resulting from
+modeling data with model, assuming Gaussian errors. Straightforward
+wrapper around MCMC and GetChisqExpr.
+
+1. data must be given as:
+{{ivar1, dvar1}, {ivar2, dvar2}, ..., {ivarN, dvarN}}
+where ivar is the independent variable, dvar is the dependent variable,
+and N is the number of data points. Either the ivars or dvars can be
+vector valued; if the independent variable is a vector, then we're just
+dealing with a function of multiple variables, and if the dependent variable
+is, then we have a vector field.
+
+2. errors must have the same length as data (= N), with each entry giving
+the errors in the corresponding dependent variable supplied in data. If
+each dvar is just a number, then so too should be each element of errors;
+if instead each dvar is a vector, then each element of errors should also
+be a vector of the same length.
+
+3. model should evaluate to either a number or a numerical vector
+(depending on dvar) when all parameters and independent variables are set.
+
+4. paramspec either gives the results of a previous MCMC run (w/ same
+plogexpr--just to add on more iterations), or lists the model parameters
+like so:
+{{param1, ival1, spread1, domain1}, ...}
+a) Each param should be symbolic.
+b) ival is the initial parameter value.
+c) spread is roughly how far to try to change the parameter each step in
+the Markov chain. In this routine we select new parameters values based
+on an exponential distribution of the form Exp[\[CapitalDelta]param/spread]. My
+Numerical Recipes book advises setting these spreads so that the average
+candidate acceptance is 10-40%.
+d) Each domain is either Reals or a list of all possible values the
+parameter can take on (needs to be a uniform grid).
+
+5. ivars gives a list of symbolic independent variables, in the same
+order as in data, on which model depends. If there's only one, then it
+need not be a list.
+
+6. numsteps is the number of Markov chain steps to perform.";
+
+GetChisqExpr::usage = "GetChisqExpr[data_List, errors_List, model_, ivars_List]
+
+Compute the chi^2 statistic for the comparison between data and model.
+(Not the reduced chi^2.)
+
+1. data must be given as:
+{{ivar1, dvar1}, {ivar2, dvar2}, ..., {ivarN, dvarN}}
+where ivar is the independent variable, dvar is the dependent variable,
+and N is the number of data points. Either the ivars or dvars can be
+vector valued; if the independent variable is a vector, then we're just
+dealing with a function of multiple variables, and if the dependent variable
+is, then we have a vector field.
+
+2. errors must have the same length as data (= N), with each entry giving
+the errors in the corresponding dependent variable supplied in data. If
+each dvar is just a number, then so too should be each element of errors;
+if instead each dvar is a vector, then each element of errors should also
+be a vector of the same length.
+
+3. ivars gives a list of symbolic independent variables, in the same
+order as in data, on which model depends. If there's only one, then it
+need not be a list.";
+
+MCMCResult::usage = "If your result object is named mcmcobj, try: mcmcobj[\"Properties\"].";
+
+Clear["MCMC`*"];
+
+Begin["`Private`"];
+
+Options[MCMC] = {
+	"BurnFraction" -> 0.1,
+	"Debug" -> False,
+	"ProgressMonitor" -> Column[{Row[{"Step", "/", "MaxSteps", "  ", TimeProgress["TimeElapsed", "DoneFraction"]}],
+		"CurrentParameters"}],
+	"ProgressInterval" -> 10,
+	"SaveTo" -> None,
+	"SaveInterval" -> 1000,
+	"FitResiduals"-> False
+};
+
+MCMC::nonnumer = "Log probability given supplied model does not evaluate to a number for initial parameters; instead evaluated to: `1`\nAbort!";
+MCMC::badinp = "Bad input: `1`.";
+MCMC[plogexpr_, paramspec_, num_Integer, opts : OptionsPattern[]] /;
+	TestMCMCInput[paramspec, num] :=
+	Module[{params, spreads, state, stateval, sets, Ns, discrete,
+		continuous, stateplog, candplog, cand, candval, hist, prevhist, prevnum,
+		prevtime, n, i, t1, t2, burn, alpha, transplog, status = "Initializing...", resume,
+		bestfitparams, z, corr, vars, plogfunc, resultlist, bestfitplot,fitresiduals},
+		Monitor[
+		If[Head[paramspec] === List, (*is user attempting to resume previous mcmc run?*)
+			resume = False; (*no*)
+
+			params = paramspec[[All, 1]];
+			stateval = paramspec[[All, 2]];
+			spreads = paramspec[[All, 3]];
+			sets = paramspec[[All, 4]];
+			prevhist = {};
+			prevnum = 0;
+			prevtime = 0;
+		,
+			resume = True; (*yes*)
+
+			params = "Parameters" /. paramspec[[1]];
+			stateval = Last["ParameterRun" /. paramspec[[1]]];
+			spreads = "ProposalSpreads" /. paramspec[[1]];
+			sets = "ParameterDomains" /. paramspec[[1]];
+
+			prevhist := Drop[Sp["ParameterRun", "ParametersLogPRun", "TransitionLogPRun"] /. paramspec[[1]], -1];
+			prevnum = Length[prevhist];
+			prevtime = "TimeSpent" /. paramspec[[1]]
+		];
+
+		n = Length[spreads]; (* # of parameters *)
+		Ns = Length /@ sets; (* size of each parameter's domain; 0 for real-valued parameters *)
+
+		status = "Evaluating chisq...";
+		plogfunc = Function[Evaluate[params], Evaluate[plogexpr]];
+   
+		discrete = Flatten[Position[sets, _List]]; (* list of parameters that are discrete valued *)
+		continuous = Complement[Range[n], discrete]; (* same, but instead continuous valued *)
+
+		cand = state = stateval; (* set initial condition *)
+
+		If[! discrete === {},
+			(* ensure discrete parameters' ICs are within their domains *)
+			state[[discrete]] = Nearest[sets[[#]] -> Range[Length[sets[[#]]]], stateval[[#]]][[1]] & /@ discrete;
+			(* convert discrete parameters' proposal dist spreads to index spreads. doesn't really work for nonuniform domains... *)
+			spreads[[discrete]] /= (Mean[GetDifferences[sets[[#]]]] & /@ discrete);
+		];
+
+		t1 = t2 = AbsoluteTime[];
+
+		status = "Initial step...";
+		If[resume,
+			stateplog = Last["ParametersLogPRun" /. paramspec[[1]]];
+		,
+			stateplog = plogfunc @@ stateval;
+			(*stateplog = Abs[plogfunc @@ stateval] (* Modified by Xisco *)*)
+		];
+		candplog = 0;
+
+		If[!NumericQ[stateplog],
+			Message[MCMC::nonnumer, stateplog];
+			Return[$Failed];
+		];
+
+		(* hist is complete run history. set initial point. *)
+		hist = Table[{0., 0., 0.}, {num}];
+
+		For[i = 2, i <= num, i++,
+			hist[[i-1]] = {stateval, stateplog, transplog(*Min[0, candplog - stateplog]*)};
+
+			(* save all information about run at intervals, if desired *)
+			If[Head[OptionValue["SaveTo"]] === String && Mod[i, OptionValue["SaveInterval"]] == 0,
+				Put[{
+						"BestFitParameters" -> Rule @@@ Sp[params, Mean[Join[prevhist, hist][[1 ;; i - 1 + prevnum, 1]]] // N],
+						"ParameterErrors" -> Rule @@@ Sp[params, StandardDeviation[Join[prevhist, hist][[1 ;; i - 1 + prevnum, 1]]] // N],
+						"AverageAcceptance" -> N[Mean[Exp[Join[prevhist, hist][[1 ;; i - 1 + prevnum, 3]]]]],
+						"TimeSpent" -> (t2 - t1) Second + prevtime,
+						"Parameters" -> params,
+						"ProposalSpreads" -> spreads,
+						"ParameterDomains" -> sets,
+						"BurnFraction" -> OptionValue["BurnFraction"],
+						"BurnEnd" -> burn,
+						"ParameterRun" -> Join[prevhist, hist][[1 ;; i - 1 + prevnum, 1]],
+						"ParametersLogPRun" -> Join[prevhist, hist][[1 ;; i - 1 + prevnum, 2]],
+						"TransitionLogPRun" -> Join[prevhist, hist][[1 ;; i - 1 + prevnum, 3]],
+						"BurnFraction" -> OptionValue["BurnFraction"]
+					}
+				,
+					OptionValue["SaveTo"]
+				]
+			];
+
+			(* update status indicator *)
+			If[Mod[i, OptionValue["ProgressInterval"]] == 0,
+				status = OptionValue["ProgressMonitor"] /. \
+					{
+						"CurrentParameters" -> Rule @@@ Sp[params, stateval],
+						"TimeElapsed" -> t2 - t1,
+						"DoneFraction" -> (i - 1)/(num),
+						"Step" -> i,
+						"MaxSteps" -> num,
+						"AverageAcceptance" -> If[And[! FreeQ[OptionValue["ProgressMonitor"], "AverageAcceptance"], i > 2],
+							Chop[N[Mean[Exp[hist[[2 ;; i-1, 3]]]]]], Null
+						]
+					}
+			];
+
+			alpha = RandomReal[{0, 1}, n]; (* random variables with which to generate candidate point *)
+			If[! continuous === {},
+				cand[[continuous]] = ExpSampleList[state[[continuous]], spreads[[continuous]], alpha[[continuous]]];
+			];
+			If[! discrete === {},
+				cand[[discrete]] = DiscExpSampleList[state[[discrete]], Ns[[discrete]], spreads[[discrete]], alpha[[discrete]]]
+			];
+
+			candval = cand;
+			If[! discrete === {},
+				candval[[discrete]] = sets[[#, cand[[#]]]] & /@ discrete;
+			];
+
+			If[OptionValue["Debug"], Print["cand ",cand," ","candval",candval]];
+
+			candplog = plogfunc @@ candval;
+			(*candplog = Abs[plogfunc @@ candval];(* Modified by Xisco *)*)
+
+
+			transplog = Min[0., candplog - stateplog +
+				If[! discrete === {}, (* discrete proposal dist is not symmetric (continuous is). take this into account. *)
+					Total[DiscExpPlog @@@ Sp[state[[discrete]], cand[[discrete]], Ns[[discrete]], spreads[[discrete]]]] -
+						Total[DiscExpPlog @@@ Sp[cand[[discrete]], state[[discrete]], Ns[[discrete]], spreads[[discrete]]]]
+				,
+					0
+				]
+			];
+
+			alpha = RandomReal[{0, 1}]; (* random variable with which to determine whether to accept candidate *)
+
+			If[OptionValue["Debug"], Print["transplog ",transplog," vs. random plog ",Log[alpha]]];
+
+			If[Log[alpha] < transplog,
+				 If[OptionValue["Debug"], Print["cand accepted; had plog ",transplog]];
+
+				 state = cand;
+				 stateval = candval;
+				 stateplog = candplog
+				 ,
+				 If[OptionValue["Debug"], Print["cand rejected; had plog ",transplog]];
+			];
+			t2 = AbsoluteTime[];
+		];
+
+		hist[[num]] = {stateval, stateplog, Min[0, candplog - stateplog]};
+
+		burn = Ceiling[Min[(num + prevnum)/2, Max[1000, (num + prevnum)*OptionValue["BurnFraction"]]]];
+
+		bestfitparams = Rule @@@ Sp[params, Mean[Join[prevhist, hist][[burn ;; num + prevnum, 1]]] // N];
+
+		status = "Computing correlation matrix...";
+		corr = Correlation[Join[prevhist, hist][[All, 1]]];
+
+		status = "Done!";
+		, status];
+
+		resultlist = {
+			"BestFitParameters" -> bestfitparams,
+			"ParameterErrors" -> Rule @@@ Sp[params, StandardDeviation[Join[prevhist, hist][[burn ;; num + prevnum, 1]]] // N],
+			"AverageAcceptance" -> N[Mean[Exp[Join[prevhist, hist][[burn ;; num + prevnum, 3]]]]],
+			"TimeSpent" -> (t2 - t1) Second + prevtime,
+			"NumSteps" -> num + prevnum,
+			"Parameters" -> params,
+			"ProposalSpreads" -> spreads,
+			"ParameterDomains" -> sets,
+			"BurnFraction" -> OptionValue["BurnFraction"],
+			"BurnEnd" -> burn,
+			"CorrelationMatrix" -> MatrixForm[corr],
+			"ParameterRun" -> Join[prevhist, hist][[All, 1]],
+			"ParametersLogPRun" -> Join[prevhist, hist][[All, 2]],
+			"TransitionLogPRun" -> Join[prevhist, hist][[All, 3]]
+		};
+		MCMCResult[resultlist]
+	];
+
+Options[MCMCModelFit] = Join[Options[MCMC], {"MakeBestFitPlot" -> False}];
+MCMCModelFit[data_List, errors_List, model_, paramspec_, vars_, num_Integer, opts : OptionsPattern[]] /;
+	TestMCMCMFInput[data, errors, model, vars] :=
+	Module[{plogexpr, result, bestfitplot},
+		(*1/2 comes from converting chi^2 to Gaussian*)
+		plogexpr = -GetChisqExpr[data, errors, model, vars] / 2;
+		
+		result = MCMC[plogexpr, paramspec, num, Sequence@@FilterRules[{opts}, Options[MCMC][[All, 1]]]];
+		
+		If[OptionValue["MakeBestFitPlot"],
+			bestfitplot = If[Length[vars] == 1,
+				Show[
+					Plot[Evaluate[model /. result["BestFitParameters"] /. vars[[1]] -> z], {z, Min[data[[All, 1]]], Max[data[[All, 1]]]}, PlotRange -> All],
+					If[Length[data[[1, 2]]] == 0,
+						ListPlot[data, Joined->False, PlotRange -> All],
+						ListPlot[Table[Sp[data[[All, 1]], data[[All, 2, i]]], {i, 1, Length[model]}], Joined->False, PlotRange -> All]
+					],
+					Frame -> True,
+					Axes -> False
+				]
+			,
+				"Number of ind. variables > 1."
+			];
+			AppendTo[result[[1]], "BestFitPlot" -> bestfitplot];
+		];
+		result
+	];
+
+DiscExpNorm = Compile[{{i, _Integer}, {NN, _Integer}, {t, _Real}},
+	(1 + Exp[1/t] - Exp[(i - NN)/t] - Exp[(1 - i)/t])/(Exp[1/t] - 1)];
+
+DiscExpNormList = Compile[{{i, _Integer, 1}, {NN, _Integer, 1}, {t, _Real, 1}},
+	(1 + Exp[1/t] - Exp[(i - NN)/t] - Exp[(1 - i)/t])/(Exp[1/t] - 1)];
+
+DiscExpPlog = Compile[{{i, _Integer}, {j, _Integer}, {NN, _Integer}, {t, _Real}},
+	-Abs[j - i]/t - Log[DiscExpNorm[i, NN, t]]];
+
+DiscExpSample = Compile[{{i, _Integer}, {NN, _Integer}, {t, _Real}, {alpha, _Real}},
+		Max[If[DiscExpNorm[i, NN, t] alpha <= Exp[1/t] (1 - Exp[-i/t])/(Exp[1/t] - 1),
+			Ceiling[t Log[1 + DiscExpNorm[i, NN, t] alpha (Exp[1/t] - 1) Exp[(i - 1)/t]]]
+		,
+			Ceiling[i - t Log[DiscExpNorm[i, NN, t] alpha (1 - Exp[1/t]) + Exp[1/t] + 1 - Exp[-(i - 1)/t]]]
+		], 1]
+	];
+
+DiscExpSampleList[i_List, NN_List, t_List, alpha_List] := DiscExpSample @@@ Transpose[{i, NN, t, alpha}];
+
+ExpSample = Compile[{state, spreads, alpha},
+	state - Sign[1/2 - alpha] spreads *	(Log[2.] + Log[Min[alpha, 1 - alpha]])
+];
+
+ExpSampleList[state_List, spreads_List, alpha_List] := ExpSample @@@ Transpose[{state, spreads, alpha}];
+
+Chisq[dpoints_List, modpoints_List, errors_List] /; Length[Dimensions[modpoints]] == 2 :=
+	Total[Flatten[(modpoints - dpoints)^2 / errors^2]];
+GetChisqExpr[data_List, errors_List, model_, invars_] :=
+	Module[{ipoints, dpoints, modpoints, modfunc, vars},
+		If[Length[invars] == 0,
+			vars = {invars},
+			vars = invars
+		];
+
+		If[NumericQ[data[[1, 1]]],
+			ipoints = List /@ data[[All, 1]],
+			ipoints = data[[All, 1]]
+		];
+
+		If[NumericQ[data[[1, 2]]],
+			dpoints = List /@ data[[All, 2]];
+			modfunc = Function[Evaluate[vars], {Evaluate[model]}];
+		,
+			dpoints = data[[All, 2]];
+			modfunc = Function[Evaluate[vars], Evaluate[model]];
+		];
+
+		modpoints = modfunc @@@ ipoints;
+
+		Chisq[dpoints, modpoints, errors]
+	];
+
+TimeLeft[timesofar_, fractiondone_] := If[fractiondone == 0., 60 * 60 * 24. - 1., timesofar * (1. / fractiondone - 1.)];
+
+Clear[TimeProgress];
+TimeProgress[timesofar_?NumericQ, fractiondone_?NumericQ] :=
+	Row[{ProgressIndicator[fractiondone],
+		", Time elapsed: " <> DateString[timesofar, {"Hour24", ":", "Minute", ":", "Second"}],
+		", Time left: " <>	DateString[TimeLeft[timesofar, fractiondone], {"Hour24", ":", "Minute", ":", "Second"}]}];
+
+Sp[x__List] /; (Equal @@ Length /@ {x}) && Length[{x}] > 1 :=
+		Transpose[{x}];
+
+(*Gets y[i+1] - y[i]*)
+GetDifferences[list_List] :=
+		Drop[(RotateLeft[list] - list), -1];
+
+TestMCMCInput[paramspec_, num_Integer] :=
+(
+	If[!MatchQ[paramspec, _MCMCResult],
+		If[!(Length[Dimensions[paramspec]] == 2 && Dimensions[paramspec][[2]] == 4 &&
+			MatchQ[paramspec[[All, 1]], {__Symbol}]),
+			Message[MCMC::badinp, "bad parameter specification"];
+			Return[False]
+		]
+	];
+
+	If[num < 2,
+		Message[MCMC::badinp, "need at least 2 steps"];
+		Return[False]
+	];
+
+	Return[True];
+);
+
+TestMCMCMFInput[data_List, errors_List, model_, vars_] :=
+(
+	If[!(MatchQ[data, {{{__?NumericQ}, {__?NumericQ}}..}] ||
+		MatchQ[data, {{_?NumericQ, {__?NumericQ}}..}] ||
+		MatchQ[data, {{_?NumericQ, _?NumericQ}..}] ||
+		MatchQ[data, {{{__?NumericQ}, _?NumericQ}..}]),
+
+		Message[MCMC::badinp, "data shaped inconsistently/incorrectly"];
+		Return[False];
+	];
+
+	If[Length[data /. _?NumericQ -> 1 // Union] > 1,
+		Message[MCMC::badinp, "data shaped inconsistently"];
+		Return[False];
+	];
+
+	If[!(data[[All, 2]] /. _?NumericQ -> 1) === (errors /. _?NumericQ -> 1),
+		Message[MCMC::badinp, "data shaped differently than errors"];
+		Return[False];
+	];
+
+	If[!If[# == 0, 1, #]&[Length[data[[1,1]]]] == If[# == 0, 1, #]&[Length[vars]],
+		Message[MCMC::badinp, "# of independent vars in data different from specified"];
+		Return[False];
+	];
+
+	Return[True];
+);
+
+Clear[MCMCResult];
+Format[MCMCResult[list_List]] := "MCMCResult"["BestFitParameters" /. list, "\[LeftSkeleton]" <> ToString[Length["ParameterRun" /. list]] <> "\[RightSkeleton]"];
+
+(this:MCMCResult[list_List])["ParameterRunPlots", opts___] :=
+	Table[
+		ListPlot[Transpose[("ParameterRun" /. list)][[i]],
+			AxesLabel -> {"Step", ToString[this["Parameters"][[i]]]},
+			FrameLabel -> {"Step", ToString[this["Parameters"][[i]]]},
+			opts
+		] // Rasterize,
+	{i, Length[this["Parameters"]]}];
+
+(this:MCMCResult[list_List])["ParameterHistograms", opts___] :=
+	Table[If[this["NumSteps"] > 1*^6,
+		SmoothHistogram[#,
+			Filling -> Axis,
+			Axes -> {True, False},
+			Ticks -> {Automatic, None},
+			AxesLabel -> {ToString[this["Parameters"][[i]]], None},
+			opts]&
+	,
+		Histogram[#,
+			Ticks -> {Automatic, None},
+			Axes -> {True, False},
+			AxesLabel -> {ToString[this["Parameters"][[i]]], None},
+			opts]&][
+		Transpose[("ParameterRun" /. list)[[this["BurnEnd"] ;; this["NumSteps"]]]][[i]]
+	], {i, Length[this["Parameters"]]}];
+
+(this:MCMCResult[list_List])["Properties"] := Join[list[[All, 1]], {"ParameterRunPlots", "ParameterHistograms"}];
+
+(this:MCMCResult[list_List])[str_String] := str /. list;
+(this:MCMCResult[list_List])[{str__String}] := Rule @@@ Sp[{str}, this /@ {str}];
+
+End[];
+
+EndPackage[];