Skip to content

Full-run helpers

Thin loops over the step functions, for when you do not need control of the recursion. Both interfaces run identical code.

from mcmckit import ram

result = ram(log_post, x0=[0.0, 0.0], n_samples=10_000)

print(result.mean(), result.std())
result.discard(1000).plot_corner()

Each helper accepts either a bare log_post callable or a Problem, and returns a Result.

metropolis 

metropolis(log_post, x0, n_samples, proposal_cov, param_names=None, rng=None)

Run random-walk Metropolis-Hastings to completion.

Parameters:

Name Type Description Default
log_post callable or Problem

theta -> log posterior, or a :class:~mcmckit.Problem.

 required
x0 (array - like, shape(d))

Starting point.

 required
n_samples int

Number of steps to take.

 required
proposal_cov float or array - like

Fixed proposal covariance.

 required
param_names sequence of str
None
rng Generator
None

Returns:

Type Description
Result
Source code in mcmckit/runners.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def metropolis(log_post, x0, n_samples, proposal_cov, param_names=None, rng=None):
    """Run random-walk Metropolis-Hastings to completion.

    Parameters
    ----------
    log_post : callable or Problem
        ``theta -> log posterior``, or a :class:`~mcmckit.Problem`.
    x0 : array-like, shape (d,)
        Starting point.
    n_samples : int
        Number of steps to take.
    proposal_cov : float or array-like
        Fixed proposal covariance.
    param_names : sequence of str, optional
    rng : numpy.random.Generator, optional

    Returns
    -------
    Result
    """
    log_post, param_names = _resolve(log_post, param_names)
    x = np.asarray(x0, dtype=float)
    cov = as_cov(proposal_cov, x.shape[0])
    logp = log_post(x)

    chain = np.empty((n_samples, x.shape[0]))
    logps = np.empty(n_samples)
    n_acc = 0
    for k in range(n_samples):
        x, logp, accepted = mh_step(log_post, x, logp, cov, rng=rng)
        chain[k] = x
        logps[k] = logp
        n_acc += bool(accepted)

    return _finish(chain, logps, param_names, n_acc / n_samples)

ram 

ram(log_post, x0, n_samples, initial_cov=None, gamma=0.51, target_rate=0.234, param_names=None, rng=None)

Run Robust Adaptive Metropolis to completion.

Self-tunes the proposal covariance, so a rough initial_cov is fine.

Returns:

Type Description
Result
Source code in mcmckit/runners.py 
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def ram(log_post, x0, n_samples, initial_cov=None, gamma=0.51,
        target_rate=0.234, param_names=None, rng=None):
    """Run Robust Adaptive Metropolis to completion.

    Self-tunes the proposal covariance, so a rough ``initial_cov`` is fine.

    Returns
    -------
    Result
    """
    log_post, param_names = _resolve(log_post, param_names)
    x = np.asarray(x0, dtype=float)
    d = x.shape[0]
    S = np.linalg.cholesky(as_cov(0.1**2 if initial_cov is None else initial_cov, d))
    logp = log_post(x)

    chain = np.empty((n_samples, d))
    logps = np.empty(n_samples)
    n_acc = 0
    for k in range(n_samples):
        x, logp, S, accepted = ram_step(
            log_post, x, logp, S, k + 1, gamma=gamma,
            target_rate=target_rate, rng=rng,
        )
        chain[k] = x
        logps[k] = logp
        n_acc += bool(accepted)

    return _finish(chain, logps, param_names, n_acc / n_samples)

dram 

dram(log_post, x0, n_samples, initial_cov=None, dr_scale=0.1, adapt_start=100, adapt_interval=10, regularization=1e-06, param_names=None, rng=None)

Run Delayed Rejection Adaptive Metropolis to completion.

Returns:

Type Description
Result
Source code in mcmckit/runners.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
def dram(log_post, x0, n_samples, initial_cov=None, dr_scale=0.1,
         adapt_start=100, adapt_interval=10, regularization=1e-6,
         param_names=None, rng=None):
    """Run Delayed Rejection Adaptive Metropolis to completion.

    Returns
    -------
    Result
    """
    log_post, param_names = _resolve(log_post, param_names)
    x = np.asarray(x0, dtype=float)
    state = init_dram_state(x, initial_cov)
    logp = log_post(x)

    chain = np.empty((n_samples, x.shape[0]))
    logps = np.empty(n_samples)
    n_acc = 0
    for k in range(n_samples):
        x, logp, state, accepted = dram_step(
            log_post, x, logp, state, dr_scale=dr_scale,
            adapt_start=adapt_start, adapt_interval=adapt_interval,
            regularization=regularization, rng=rng,
        )
        chain[k] = x
        logps[k] = logp
        n_acc += bool(accepted)

    return _finish(chain, logps, param_names, n_acc / n_samples)

mala 

mala(log_post_and_grad, x0, n_samples, step_size, param_names=None, rng=None)

Run MALA to completion.

Parameters:

Name Type Description Default
log_post_and_grad callable or Problem

theta -> (log posterior, gradient), or a Problem carrying gradients.

 required

Returns:

Type Description
Result
Source code in mcmckit/runners.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def mala(log_post_and_grad, x0, n_samples, step_size, param_names=None, rng=None):
    """Run MALA to completion.

    Parameters
    ----------
    log_post_and_grad : callable or Problem
        ``theta -> (log posterior, gradient)``, or a Problem carrying
        gradients.

    Returns
    -------
    Result
    """
    f, param_names = _resolve_grad(log_post_and_grad, param_names)
    x = np.asarray(x0, dtype=float)
    logp, grad = f(x)

    chain = np.empty((n_samples, x.shape[0]))
    logps = np.empty(n_samples)
    n_acc = 0
    for k in range(n_samples):
        x, logp, grad, accepted = mala_step(f, x, logp, grad, step_size, rng=rng)
        chain[k] = x
        logps[k] = logp
        n_acc += bool(accepted)

    return _finish(chain, logps, param_names, n_acc / n_samples)

adaptive_mala 

adaptive_mala(log_post_and_grad, x0, n_samples, initial_step_size=0.1, gamma=0.51, target_rate=0.574, param_names=None, rng=None)

Run MALA with log-space step-size adaptation to completion.

Returns:

Type Description
Result
Source code in mcmckit/runners.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def adaptive_mala(log_post_and_grad, x0, n_samples, initial_step_size=0.1,
                  gamma=0.51, target_rate=0.574, param_names=None, rng=None):
    """Run MALA with log-space step-size adaptation to completion.

    Returns
    -------
    Result
    """
    f, param_names = _resolve_grad(log_post_and_grad, param_names)
    x = np.asarray(x0, dtype=float)
    logp, grad = f(x)
    log_step = np.log(initial_step_size)

    chain = np.empty((n_samples, x.shape[0]))
    logps = np.empty(n_samples)
    n_acc = 0
    for k in range(n_samples):
        x, logp, grad, log_step, accepted = adaptive_mala_step(
            f, x, logp, grad, log_step, k + 1,
            gamma=gamma, target_rate=target_rate, rng=rng,
        )
        chain[k] = x
        logps[k] = logp
        n_acc += bool(accepted)

    return _finish(chain, logps, param_names, n_acc / n_samples)

gibbs 

gibbs(log_post, x0, n_samples, blocks=None, proposal_std=0.5, param_names=None, rng=None)

Run Metropolis-within-Gibbs to completion.

Parameters:

Name Type Description Default
blocks sequence of sequence of int

Index groups updated together. Defaults to one block per parameter.

 None

Returns:

Type Description
Result

acceptance_rate is the mean over blocks. Per-block rates are on the returned object as block_acceptance_rates.
Source code in mcmckit/runners.py 
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def gibbs(log_post, x0, n_samples, blocks=None, proposal_std=0.5,
          param_names=None, rng=None):
    """Run Metropolis-within-Gibbs to completion.

    Parameters
    ----------
    blocks : sequence of sequence of int, optional
        Index groups updated together. Defaults to one block per parameter.

    Returns
    -------
    Result
        ``acceptance_rate`` is the mean over blocks. Per-block rates are on
        the returned object as ``block_acceptance_rates``.
    """
    log_post, param_names = _resolve(log_post, param_names)
    x = np.asarray(x0, dtype=float)
    d = x.shape[0]
    if blocks is None:
        blocks = [[i] for i in range(d)]
    logp = log_post(x)

    chain = np.empty((n_samples, d))
    logps = np.empty(n_samples)
    n_acc = np.zeros(len(blocks))
    for k in range(n_samples):
        x, logp, accepted = gibbs_step(
            log_post, x, logp, blocks, proposal_std, rng=rng
        )
        chain[k] = x
        logps[k] = logp
        n_acc += np.asarray(accepted)

    rates = n_acc / n_samples
    result = _finish(chain, logps, param_names, float(rates.mean()))
    result.block_acceptance_rates = rates
    return result