Skip to content

resampling

MC_resample(metamodel, nresamps=3000)

Use Monte Carlo resampling to estimate whether ML parameters are accurate or not. See section 5 here for details: https://github.com/sievers/phys512-2025/blob/main/notes/mcmc.pdf

Parameters:

Name Type Description Default
metamodel MetaModel

Metamodel containing our fit model/parameters

 required
nresamps int, default; 3000

Number of points to resample at.

 3000

Returns:

Name Type Description
samps ndarray

Array of resampled parameters
likelihood_chain ndarray

Array of likelihood ratios
Source code in witch/resampling.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def MC_resample(
    metamodel: MetaModel,
    nresamps: int = 3000,
) -> (np.ndarray, np.ndarray):
    """
    Use Monte Carlo resampling to estimate whether ML parameters are accurate or not.
    See section 5 here for details:
    https://github.com/sievers/phys512-2025/blob/main/notes/mcmc.pdf

    Parameters
    ----------
    metamodel : MetaModel
        Metamodel containing our fit model/parameters
    nresamps : int, default; 3000
        Number of points to resample at.

    Returns
    -------
    samps : np.ndarray
        Array of resampled parameters
    likelihood_chain : np.ndarray
        Array of likelihood ratios
    """
    if not np.any(metamodel.cov):
        raise RuntimeError(
            "Error: cov is not computed. You must fit data before doing resampling"
        )
    likelihood_chain = np.zeros(
        nresamps
    )  # Contains the likelihood_ratio, L_exact/L_gaussian
    samps = np.ones((nresamps, len(metamodel.to_fit)))
    samps *= np.array(metamodel.parameters)
    temp_samps = draw_samps(
        par_means=metamodel.parameters[metamodel.to_fit],
        cov=metamodel.cov[metamodel.to_fit].T[metamodel.to_fit],
        nsamps=nresamps,
    ).T
    samps[..., metamodel.to_fit] = temp_samps

    inv_cov = invsafe(metamodel.cov)
    base_loglike, _, _ = joint_objective(
        metamodel=metamodel, do_loglike=True, do_grad=False, do_curve=False
    )
    for i in range(nresamps):
        cur_pars = copy(metamodel.parameters).at[:].set(samps[i])
        cur_meta = copy(metamodel).update(
            pars=cur_pars, errs=metamodel.errs, cov=metamodel.cov, chisq=metamodel.chisq
        )
        loglike, _, _ = joint_objective(
            metamodel=cur_meta, do_loglike=True, do_grad=False, do_curve=False
        )
        l_exact = np.exp(-1 / 2 * (loglike - base_loglike))
        delta_p = samps[i] - metamodel.parameters
        l_gauss = np.exp(-1 / 2 * np.dot(delta_p.T, np.dot(inv_cov, delta_p)))
        likelihood_chain[i] = l_exact / l_gauss

    return samps, likelihood_chain

draw_samps(par_means, cov, nsamps=3000)

Given a maximum likelihood computed parameter means and covariance, draw a sample assuming gaussian distribution under that mean, cov. See this stack exchange on how to generate samples from a given cov, p: https://stats.stackexchange.com/questions/120179/generating-data-with-a-given-sample-covariance-matrix

Parameters:

Name Type Description Default
par_means Array

Array of parameters means

 required
cov Array

Covariance of parameters

 required
nsamps int, default; 3000

Number of samples to draw

 3000

Returns:

Name Type Description
X Array

Sample of p drawn from par_means, cov.
Source code in witch/resampling.py 
 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
104
105
106
107
108
109
110
def draw_samps(par_means: jax.Array, cov: jax.Array, nsamps: int = 3000) -> jax.Array:
    """
    Given a maximum likelihood computed parameter means and covariance,
    draw a sample assuming gaussian distribution under that mean, cov.
    See this stack exchange on how to generate samples from a given
    cov, p: https://stats.stackexchange.com/questions/120179/generating-data-with-a-given-sample-covariance-matrix

    Parameters
    ----------
    par_means : jax.Array
        Array of parameters means
    cov : jax.Array
        Covariance of parameters
    nsamps : int, default; 3000
        Number of samples to draw

    Returns
    -------
    X : jax.Array
        Sample of p drawn from par_means, cov.
    """
    key = jax.random.key(seed=42)
    X = jax.random.multivariate_normal(key=key, mean=par_means, cov=cov, shape=nsamps).T
    for n in range(X.shape[0]):
        X = X.at[n].set(X[n] - X[n].mean())

    # Make each variable in X orthogonal to one another
    L_inv = jnp.linalg.cholesky(jnp.cov(X))
    L_inv = jnp.linalg.inv(L_inv)
    X = jnp.dot(L_inv, X)

    # Rescale X to exactly match Sigma
    L = jnp.linalg.cholesky(cov)
    X = jnp.dot(L, X)

    # Add the mean back into each variable
    for n in range(X.shape[0]):
        X = X.at[n].set(X[n] + par_means[n])

    return X