Skip to content

noise_proposal

noise_variance_probability(kernel_gram, current_noise_variance, proposed_noise_variance, data)

Calculate the posterior probability for a proposed noise variance, given the current kernel gram matrix and current noise variance.

See equation 21 in Saad2023.

Parameters:

Name Type Description Default
kernel_gram Float[ndarray, 'D D']

The gram matrix of the kernel.

 required
current_noise_variance ScalarFloat

The current noise variance.

 required
proposed_noise_variance ScalarFloat

The proposed noise variance.

 required
data Dataset

The data containing the observations, as a Dataset object.

 required

Returns:

Type Description
ScalarFloat

The posterior probability for the proposed noise variance.
Source code in gallifrey/moves/noise_proposal.py 
 64
 65
 66
 67
 68
 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
104
105
106
107
108
109
110
111
112
113
114
def noise_variance_probability(
    kernel_gram: Float[jnp.ndarray, "D D"],
    current_noise_variance: ScalarFloat,
    proposed_noise_variance: ScalarFloat,
    data: Dataset,
) -> ScalarFloat:
    """
    Calculate the posterior probability for a
    proposed noise variance, given the current
    kernel gram matrix and current noise variance.

    See equation 21 in Saad2023.

    Parameters
    ----------
    kernel_gram : Float[jnp.ndarray, "D D"]
        The gram matrix of the kernel.
    current_noise_variance : ScalarFloat
        The current noise variance.
    proposed_noise_variance : ScalarFloat
        The proposed noise variance.
    data : Dataset
        The data containing the observations,
        as a Dataset object.

    Returns
    -------
    ScalarFloat
        The posterior probability for the proposed noise variance.
    """

    # calculate Σ⁻¹y | [n,1]
    inv_sigma_r = calculate_inv_sigma_r(
        data.y,
        kernel_gram,
        current_noise_variance,
    )

    # calculate the means of the Gaussian process, µ = Kxx*Σ⁻¹y | [n,1]
    means = kernel_gram @ inv_sigma_r

    # create inverse gamma distribution
    diff = data.y - means
    inv_gamma = InverseGamma(
        concentration=1 + data.y.size / 2,
        scale=1 + 0.5 * diff @ diff,
    )
    # calculate the log probability of the proposal
    log_prob_proposal = inv_gamma.log_prob(proposed_noise_variance)

    return log_prob_proposal

noise_variance_proposal(key, kernel_gram, noise_variance, data)

Sample a new proposal for the noise variance, and the log probability of the proposal.

See equation 21 in Saad et al. 2023.

Parameters:

Name Type Description Default
key PRNGKeyArray

The random key.

 required
kernel_gram Float[ndarray, 'D D']

The gram matrix of the kernel.

 required
noise_variance ScalarFloat

The current noise variance.

 required
data Dataset

The data containing the observations, as a Dataset object.

 required

Returns:

Type Description
ScalarFloat

The new noise variance.
ScalarFloat

The log probability of the proposal.
Source code in gallifrey/moves/noise_proposal.py 
10
11
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
def noise_variance_proposal(
    key: PRNGKeyArray,
    kernel_gram: Float[jnp.ndarray, "D D"],
    noise_variance: ScalarFloat,
    data: Dataset,
) -> tuple[ScalarFloat, ScalarFloat]:
    """
    Sample a new proposal for the noise variance,
    and the log probability of the proposal.

    See equation 21 in Saad et al. 2023.

    Parameters
    ----------
    key : PRNGKeyArray
        The random key.
    kernel_gram : Float[jnp.ndarray, "D D"]
        The gram matrix of the kernel.
    noise_variance : ScalarFloat
        The current noise variance.
    data : Dataset
        The data containing the observations,
        as a Dataset object.

    Returns
    -------
    ScalarFloat
        The new noise variance.
    ScalarFloat
        The log probability of the proposal.
    """

    # calculate Σ⁻¹y | [n,1]
    inv_sigma_r = calculate_inv_sigma_r(
        data.y,
        kernel_gram,
        noise_variance,
    )

    # calculate the means of the Gaussian process, µ = Kxx*Σ⁻¹y | [n,1]
    means = kernel_gram @ inv_sigma_r

    # sample new noise variance from inverse gamma distribution
    diff = data.y - means
    inv_gamma = InverseGamma(
        concentration=1 + data.y.size / 2,
        scale=1 + 0.5 * diff @ diff,
    )
    new_noise_variance = inv_gamma.sample(seed=key, sample_shape=())
    log_prob = inv_gamma.log_prob(new_noise_variance)

    return new_noise_variance, log_prob  # type: ignore