UAH LTG IDL Library

IDL Routines from Phillip Bitzer and UAH Lightning Group

summary     class     fields     routine details     file attributes

ats_mcmc__define.pro

ATS606, Statistics, Monte Carlo, Markov Chain, MCMC

includes main-level program

This object can be used to estimate the PDFs of parameters that best fit a model using the Markov Chain Monte Carlo (MCMC) method.

There are several advantages to use a MCMC to do this:

 No local minimization problems 
This particular MCMC calculates the likelihood assuming Gaussian Probabilites (so the likelihood is a function of $\chi^2$). Further, this uses uniform priors and uniform proposal distributions.

You can fit to any model you wish. The model should be a function that accepts two parameters; the first parameter is the array of independent data, and the second parameter is the array of parameters of the model. It is called using CALL_EXTERNAL. This file contains an example of this ats_mcmc_model_linear.

One limitation is this object uses a predefined number of steps in the chain.

Of course, MCMCs are discussed in detail in ATS606, so see the notes for more information.

Examples

This object contains "testing" data to get a feel for how it works. To use it:

 m = ats_mcmc() m->StepAll m->GetVars, param=p 
Running these commands will a fit to a linear model ($y=A + Bx$) using simiulated data. The data is simulated with $A=5$ and $B=2$ with a random Gaussian error thrown in.

To use your own data, you should provide the independent and dependent data as parameters when intializing the function. Further, you'll need to set up the priors and proposals.

Author information

Author

Phillip M. Bitzer, University of Alabama in Huntsville, pm.bitzer "AT" uah.edu

History

Modification History:

 First written: Apr 16, 2014 See Git log for incremental changes. 

Uses:

Guess init
Priors init
Proposals init
modelFct init
nSteps init

Routines

Routines from ats_mcmc__define.pro

result = ats_mcmc_model_quadratic(x, param)

This provides an(other) (internal) example of how the model function should be defined.

result = ats_mcmc_model_linear(x, param)

This provides an (internal) example of how the model function should be defined.

ats_mcmc::PlotChains [, PNG= boolean] [, XRANGE=array] [, NOERASE=NOERASE]

This method will plot a "time" series of the chain for each parameter and the likelihood (well, $\chi^2$).

ats_mcmc::GetStats [, MEANS=array] [, STDEVS=array] [, /Print]

This method will get some rundimentary statistics about the chain.

ats_mcmc::PlotHistos [, PNG= boolean]

This method will produce a histogram of the PDF for each parameter, using the thinned chain.

ats_mcmc::PlotFit [, PNG= boolean]

This method will plot the data provided and the "best fit" parameters to the model.

result = ats_mcmc::thinChainInd( [/NoBurn])

This method will go get the indices of a chain that should effectively "thin" the chain.

ats_mcmc::GetVars [, LL=array] [, Params=array]

This method "gets" properties of the object.

result = ats_mcmc::AcceptStep(llhood)

This method decides whether or not to accept the cnadidate link in the chain, based on the passed likelihood and the likelihood at the previous step.

result = ats_mcmc::GetLikelihood(modelData)

Given some model data, calcluate the likelihood.

result = ats_mcmc::CalcModelData(param)

Given the passed parameters, go calculate the model data.

ats_mcmc::stepAll

This method will step through the rest of the MCMC, from wherever you are now.

ats_mcmc::step

The method makes one "step" (or link) in the MCMC.

result = ats_mcmc::GetProposal()

Get a set of proposal candidates for (ideally) the next step in the chain.

result = ats_mcmc::GetPrior(params [, /GoodStep])

This returns the probability of parameters using the prior distribution.

ats_mcmc::cleanup

Cleanup the object.

result = ats_mcmc::init( [indptVar] [, data] [, modelFct=string] [, Priors=p element array of structure] [, Proposals=p element array of structure] [, nSteps=integer] [, Guess=array])

Initializes the ATS_MCMC object.

ats_mcmc__define

Define the ATS_MCMC object.

Routine details

result = ats_mcmc_model_quadratic(x, param)

This provides an(other) (internal) example of how the model function should be defined. For this example, we're using a polynomial model of order two: $y=A + Bx + Cx^2$.

Return value

A n element array, where n is the number of data points. On error, NaN is returned.

Parameters

x in required type=array

An array of n elements for the independent data.

param in required type=type

An array of p elements for the parameters to be used. For this example, this a particular $A$, $B$, and $C$.

top ats_mcmc_model_linear

result = ats_mcmc_model_linear(x, param)

This provides an (internal) example of how the model function should be defined. For this example, we're using a linear model $y=A + Bx$.

Return value

A n element array, where n is the number of data points. On error, NaN is returned.

Parameters

x in required type=array

An array of n elements for the independent data.

param in required type=type

An array of p elements for the parameters to be used. For this example, this a particular $A$ and $B$.

top ats_mcmc::PlotChains

ats_mcmc::PlotChains [, PNG= boolean] [, XRANGE=array] [, NOERASE=NOERASE]

This method will plot a "time" series of the chain for each parameter and the likelihood (well, $\chi^2$). For visual comparison, the thinned parts of the chain are marked in green.

Keywords

PNG in optional type= boolean default= 0B

If set, a PNG will be created named 'ats_mcmc_chains.png'

XRANGE in optional type=array

Set this to an array of [xmin, xmax] to control the X axis range. If not provided, IDL will handle it.

NOERASE in optional default=0B

If set, the display will not be erased before plotting. No affect if png is set.

Uses:

top ats_mcmc::GetStats

ats_mcmc::GetStats [, MEANS=array] [, STDEVS=array] [, /Print]

This method will get some rundimentary statistics about the chain. You can "get" the means and standard deviations, of both the thinned and unthinned chains.

You can also print some more stats, including the acceptance rate and the average chi squared value.

Keywords

MEANS out optional type=array

If a variable is provided, this output keyword will return a px2 array, where p is the number of parameters, of the means of each chain. The thinned means are in the 0th index (of the second dimension) and the unthinned means are in the 1st index.

STDEVS out optional type=array

If a variable is provided, this output keyword will return a px2 array, where p is the number of parameters, of the standard deviations of each chain. The thinned standard deviations are in the 0th index (of the second dimension) and the unthinned standard deviations are in the 1st index.

Print in optional type=boolean default=0B

If set, statistics are printed to the console.

top ats_mcmc::PlotHistos

ats_mcmc::PlotHistos [, PNG= boolean]

This method will produce a histogram of the PDF for each parameter, using the thinned chain.

Keywords

PNG in optional type= boolean default= 0B

If set, a PNG will be created named 'ats_mcmc_histo.png'

Uses:

top ats_mcmc::PlotFit

ats_mcmc::PlotFit [, PNG= boolean]

This method will plot the data provided and the "best fit" parameters to the model. Currently, it uses the mean of the chain for each parameter.

Keywords

PNG in optional type= boolean default= 0B

If set, a PNG will be created named 'ats_mcmc_best_fit.png'

Uses:

top ats_mcmc::thinChainInd

result = ats_mcmc::thinChainInd( [/NoBurn])

This method will go get the indices of a chain that should effectively "thin" the chain. It does so by selecting every 10th link, excluding (nominally) the burn-in.

Return value

An array of indices. This array can be used to index the original array to effectively thin the chain.

Keywords

NoBurn in optional type=boolean default=0B

Does nothing for now. Might be enabled to allow thining including the burn in, although I can't see why you'd want that.

top ats_mcmc::GetVars

ats_mcmc::GetVars [, LL=array] [, Params=array]

This method "gets" properties of the object.

Keywords

LL out optional type=array

If set, the likelihood at each step is returned as an n element array, where n is the number of steps in the chain.

Params in optional type=array

If set, the parameters at each step in the chain are returned as a nxp element array, where n is the number of steps in the chain and p is the number of parameters.

top ats_mcmc::AcceptStep

result = ats_mcmc::AcceptStep(llhood)

This method decides whether or not to accept the cnadidate link in the chain, based on the passed likelihood and the likelihood at the previous step.

Return value

Returns 1 if we should accept the link, 0 if not.

Parameters

llhood in required type=float

The likelihood you wish to evaluate. This will be compared to the likelihood in the previous step. And although this is called the likelihood, it really is the $\chi^2$ value.

top ats_mcmc::GetLikelihood

result = ats_mcmc::GetLikelihood(modelData)

Given some model data, calcluate the likelihood. For this implementation of a MCMC, "likelihood" really means "$\chi^2$".

Return value

The likelihood, er the $\chi^2$

Parameters

modelData in required type=array

An array of "modelled data."

top ats_mcmc::CalcModelData

result = ats_mcmc::CalcModelData(param)

Given the passed parameters, go calculate the model data. This is done by using CALL_FUNCTION to call the function.

The model data.

Parameters

param in required type=array

The parameters of the model to be used to calculate the model data.

top ats_mcmc::stepAll

ats_mcmc::stepAll

This method will step through the rest of the MCMC, from wherever you are now. Basically, this calls the method Step until we reach the end of the chain.

top ats_mcmc::step

ats_mcmc::step

The method makes one "step" (or link) in the MCMC. This is the heart of the MCMC. It calls several of the other methods in turn:

1) Generate proposal (GetProposal) 2) Asses prior probabilites (GetPrior) 3) Calculate the model data for the candidate parameters (CalcModelData) 4) Get the likelihood for the candidate (GetLikelihood) 5) Decide to accept or reject the candidate (AcceptStep)

If accepted, the candidate is entered into the chain. If not, the last link is repeated.

Finally, the property iStep is incremented.

top ats_mcmc::GetProposal

result = ats_mcmc::GetProposal()

Get a set of proposal candidates for (ideally) the next step in the chain.

Return value

A p element array, where p is the number of parameters. On error, NaN is returned.

top ats_mcmc::GetPrior

result = ats_mcmc::GetPrior(params [, /GoodStep])

This returns the probability of parameters using the prior distribution. Right now, uniform priors are assumed. On error, NaN is returned.

Return value

Since we're assuming uniform priors, 1 is returned if all the parameters are within bounds, 0 if not. Also, goodStep is set to 0 if the parameters are not within bounds.

Parameters

params in required type=array

The array of parameters.

Keywords

GoodStep out optional type=boolean

Whether or not the parameters are within bounds.

top ats_mcmc::cleanup

ats_mcmc::cleanup

Cleanup the object. Really, free all the pointers associated with the object.

top ats_mcmc::init

result = ats_mcmc::init( [indptVar] [, data] [, modelFct=string] [, Priors=p element array of structure] [, Proposals=p element array of structure] [, nSteps=integer] [, Guess=array])

Initializes the ATS_MCMC object. See the file documentation for more.

Return value

1 on success. 0 if something went wrong.

Parameters

indptVar in optional type=n element array

The independent data. If both this and data are not given, I'll initialize using some test data.

data in optional type=n element array

The dependent data. If both this and indptVar are not given, I'll initialize using some test data.

Keywords

modelFct in optional type=string default=ats_mcmc_model_linear

The model you wish to fit the model to. See the full documentation for the description. By default, a linear model is used.

Priors in optional type=p element array of structure

The array of structures that describes the prior distribution for each parameter, assuming p parameters.

Technically, this is optional, but it only gets defined if you don't provide two parameters. So, it's pretty much required.

Proposals in optional type=p element array of structure

The array of structures that describes the proposal distribution for each parameter, assuming p parameters.

Each structure should have one tag, step, which describes the "half-width" of the proposal. distribution. For eaxample:

proposal1 = {step:step1} proposal1 = {step:step1} proposals = [proposal1, proposal2] 
Technically, this is optional, but it only gets defined if you don't provide two parameters. So, it's pretty much required.

nSteps in optional type=integer default=10000

The total number of steps in the MCMC chain.

Guess in optional type=array default=varies

The initial starting point for the chain. It should be a p element array, where p is the number of parameters. If it's not, then the starting point is in the middle of the prior distribution.

top ats_mcmc__define

ats_mcmc__define

Define the ATS_MCMC object.

File attributes

 Modification date: Wed May 14 18:57:19 2014 Lines: 365 Docformat: rst rst