# drawSamples

**Class: **HamiltonianSampler

Generate Markov chain using Hamiltonian Monte Carlo (HMC)

## Syntax

`chain = drawSamples(smp)`

[chain,endpoint,accratio]
= drawSamples(smp)

[chain,endpoint,accratio]
= drawSamples(___,Name,Value)

## Description

generates
a Markov chain by drawing samples using the Hamiltonian Monte Carlo
sampler `chain`

= drawSamples(`smp`

)`smp`

.

`[`

also returns the
final state of the Markov chain in `chain`

,`endpoint`

,`accratio`

]
= drawSamples(`smp`

)`endpoint`

and
the fraction of accepted proposals in `accratio`

.

`[`

specifies
additional options using one or more name-value pair arguments. Specify
name-value pair arguments after all other input arguments.`chain`

,`endpoint`

,`accratio`

]
= drawSamples(___,`Name,Value`

)

## Input Arguments

`smp`

— Hamiltonian Monte Carlo sampler

`HamiltonianSampler`

object

Hamiltonian Monte Carlo sampler, specified as a `HamiltonianSampler`

object.

`drawSamples`

draws samples from the target
log probability density in `smp.LogPDF`

. Use the `hmcSampler`

function to create a sampler.

### Name-Value Arguments

Specify optional pairs of arguments as
`Name1=Value1,...,NameN=ValueN`

, where `Name`

is
the argument name and `Value`

is the corresponding value.
Name-value arguments must appear after other arguments, but the order of the
pairs does not matter.

*
Before R2021a, use commas to separate each name and value, and enclose*
`Name`

*in quotes.*

**Example: **`'Burnin',500,'NumSamples',2000`

generates
a Markov chain by discarding 500 burn-in samples and then drawing
2000 samples.

`Burnin`

— Number of burn-in samples to discard

`1000`

(default) | positive integer

Number of burn-in samples to discard from the beginning of the Markov chain, specified as a positive integer.

**Example: **`'Burnin',500`

`NumSamples`

— Number of samples to draw

`1000`

(default) | positive integer

Number of samples to draw from the Markov chain using the HMC sampler, specified as a positive integer.

The `drawSamples`

method generates this number
of samples after the burn-in period.

**Example: **`'NumSamples',2000`

`ThinSize`

— Markov chain thinning size

`1`

(default) | positive integer

Markov chain thinning size, specified as a positive integer.

Only one out of the `'ThinSize'`

number of
samples are kept. The rest of the samples are discarded.

**Example: **`'ThinSize',5`

`StartPoint`

— Initial point to start sampling from

`smp`

`.StartPoint`

(default) | numeric column vector

Initial point to start sampling from, specified as a numeric
column vector with the same number of elements as the `StartPoint`

property
of the sampler `smp`

.

**Example: **`'StartPoint',randn(5,1)`

`VerbosityLevel`

— Verbosity level of Command Window output

`0`

(default) | positive integer

Verbosity level of Command Window output during sampling, specified
as `0`

or a positive integer.

With the value set to `0`

, `drawSamples`

displays
no details during sampling.

With the value set to a positive integer, `drawSamples`

displays
details of the sampling. To set the output frequency, use the `'NumPrint'`

name-value
pair argument.

`drawSamples`

displays the output as a table
with these columns.

Heading | Description |
---|---|

`ITER` | Iteration number |

`LOG PDF` | Log probability density at the current iteration |

`STEP SIZE` | Leapfrog integration step size at the current iteration. If the step size is jittered, it can vary between iterations. |

`NUM STEPS` | Number of leapfrog integration steps at the current iteration. If the number of steps is jittered, it can vary between iterations |

`ACC RATIO` | Acceptance ratio, that is, the fraction of proposals that are accepted. The acceptance ratio is calculated from the beginning of sampling, including the burn-in period. |

`DIVERGENT` | Number of times the sampler failed to generate a valid
proposal due to the leapfrog iterations generating |

**Example: **`'VerbosityLevel',1`

`NumPrint`

— Verbose output frequency

`100`

(default) | positive integer

Verbose output frequency, specified as a positive integer.

If the `'VerbosityLevel'`

value is a positive
integer, then `drawSamples`

outputs sampling details
every `'NumPrint'`

iterations.

**Example: **`'NumPrint',200`

## Output Arguments

`chain`

— Markov chain generated using Hamiltonian Monte Carlo

numeric matrix

Markov chain generated using Hamiltonian Monte Carlo, returned as a numeric matrix.

Each row of `chain`

is a sample, and each
column represents one sampling variable.

`endpoint`

— Final state of Markov chain

numeric column vector

Final state of the Markov chain, returned as a numeric column
vector of the same length as `smp.StartPoint`

.

`accratio`

— Acceptance ratio

numeric scalar

Acceptance ratio of the Markov chain proposals, returned as a numeric scalar. The acceptance ratio is calculated from the beginning of sampling, including the burn-in period.

## Examples

### Draw Samples Using HMC Sampler

Create MCMC chains for a multivariate normal distribution using a Hamiltonian Monte Carlo (HMC) sampler.

Define the number of parameters to sample and their means.

NumParams = 100; means = randn(NumParams,1); standevs = 0.1;

First, save a function `normalDistGrad`

on the MATLAB® path that returns the multivariate normal log probability density and its gradient (`normalDistGrad`

is defined at the end of this example). Then, call the function with arguments to define the `logpdf`

input argument to the `hmcSampler`

function.

logpdf = @(theta)normalDistGrad(theta,means,standevs);

Choose a starting point of the sampler. Create the HMC sampler and tune its parameters.

startpoint = randn(NumParams,1); smp = hmcSampler(logpdf,startpoint); smp = tuneSampler(smp);

Draw samples from the posterior density, using a few independent chains. Choose different, randomly distributed starting points for each chain. Specify the number of burn-in samples to discard from the beginning of the Markov chain and the number of samples to generate after the burn-in. Set the `'VerbosityLevel'`

to print details during sampling for the first chain.

NumChains = 4; chains = cell(NumChains,1); Burnin = 500; NumSamples = 2000; for c = 1:NumChains if c == 1 showOutput = 1; else showOutput = 0; end chains{c} = drawSamples(smp,'Burnin',Burnin,'NumSamples',NumSamples,... 'Start',randn(size(startpoint)),'VerbosityLevel',showOutput,'NumPrint',500); end

|==================================================================================| | ITER | LOG PDF | STEP SIZE | NUM STEPS | ACC RATIO | DIVERGENT | |==================================================================================| | 500 | 8.450463e+01 | 4.776e-01 | 5 | 9.060e-01 | 0 | | 1000 | 8.034444e+01 | 4.776e-01 | 9 | 8.810e-01 | 0 | | 1500 | 9.156276e+01 | 4.776e-01 | 2 | 8.867e-01 | 0 | | 2000 | 8.027782e+01 | 2.817e-02 | 6 | 8.890e-01 | 0 | | 2500 | 9.892440e+01 | 4.648e-01 | 2 | 8.904e-01 | 0 |

After obtaining a random sample, investigate issues such as convergence and mixing to determine whether the samples represent a reasonable set of random realizations from the target distribution. To examine the output, plot the trace plots of the samples for the first few variables using the first chain.

A number of burn-in samples have been removed to reduce the effect of the sampling starting point. Furthermore, the trace plots look like high-frequency noise, without any visible long-range correlation between the samples. This indicates that the chain is well mixed.

for p = 1:3 subplot(3,1,p); plot(chains{1}(:,p)); ylabel(smp.VariableNames(p)) axis tight end xlabel('Iteration')

The `normalDistGrad`

function returns the logarithm of the multivariate normal probability density with means in `Mu`

and standard deviations in `Sigma`

, specified as scalars or columns vectors the same length as the `startpoint`

. The second output argument is the corresponding gradient.

function [lpdf,glpdf] = normalDistGrad(X,Mu,Sigma) Z = (X - Mu)./Sigma; lpdf = sum(-log(Sigma) - .5*log(2*pi) - .5*(Z.^2)); glpdf = -Z./Sigma; end

## Tips

After creating an HMC sampler using the

`hmcSampler`

function, you can compute MAP (maximum-a-posteriori) point estimates, tune the sampler, draw samples, and check convergence diagnostics using the methods of the`HamiltonianSampler`

class. For an example of this workflow, see Bayesian Linear Regression Using Hamiltonian Monte Carlo.

## Version History

**Introduced in R2017a**

## See Also

### Functions

### Classes

## Comando MATLAB

Hai fatto clic su un collegamento che corrisponde a questo comando MATLAB:

Esegui il comando inserendolo nella finestra di comando MATLAB. I browser web non supportano i comandi MATLAB.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)