## ats606/

# ats_mcmc__define.pro

ATS606, Statistics, Monte Carlo, Markov Chain, MCMC

includes main-level programThis 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
```

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
```

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.`

## Other file information

- Uses:

## Class description for ats_mcmc

### Properties

#### Properties in ats_mcmc

- 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

## top ats_mcmc_model_quadratic

`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.

### Other attributes

- 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'

### Other attributes

- 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'

### Other attributes

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

### Return value

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:

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

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

## File attributes

Modification date: | Wed May 14 18:57:19 2014 |

Lines: | 365 |

Docformat: | rst rst |