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 theHamiltonianSampler
class. For an example of this workflow, see Bayesian Linear Regression Using Hamiltonian Monte Carlo.
Version History
Introduced in R2017a
See Also
Functions
Classes
Commande MATLAB
Vous avez cliqué sur un lien qui correspond à cette commande MATLAB :
Pour exécuter la commande, saisissez-la dans la fenêtre de commande de MATLAB. Les navigateurs web ne supportent pas les commandes 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)