Statistics

Figure 1 is an example of a notation for graphical models.

Consider a posterior distribution where the mode is at \(\theta^{*}\) and the Hessian matrix at this point is \(h''(\theta^{*})\). If we consider a multivariate Gaussian density, we can see that it will have its mode at the mean, \(\mu\), and, assuming a covariance matrix, \(\Sigma\), the density will have a Hessian of \(-\Sigma^{-1}\) at the mode. This motivates the Gaussian approximation at the mode of the posterior. Also see Fisher information.

See Chib and Greenberg (1995) for a timeless introduction to the Metropolis-Hasting algorithm.

For a random variable, \(X\) with distrtibution \(f(X;\theta)\) the score function is

\[ s(X; \theta) = \frac{\partial}{\partial\theta} \log f(X:\theta). \]

Under some mild assumptions, the mean value of \(s(X; \theta)\) (when viewed as a function of \(X\)) is zero.

The Fisher information is defined as the variance of the score and denoted \(\mathcal{I}(\theta)\):

\[ \mathcal{I}(\theta) = \mathbb{E}\left[ \left(\frac{\partial}{\partial\theta} \log f(X:\theta) \right)^{2} \right]. \]

It can be shown that the Fisher information is equal to

\[ -\mathbb{E}\left[ \frac{\partial^{2}}{\partial\theta^{2}} \log f(X:\theta) \right] \]

which can be more convenient to work with.

Let \(X_1,\dots,X_n\sim F\) IID on \(\mathbb{R}\), then the empirical distribution function (EDF) is

\[ \hat{F}_n(x) = \frac{1}{n}\sum I(X_i \leq x) \]

where \(I\) is the indicator function. Linearity shows

\[ \mathbb{E}\hat{F}_n(x) = F(x) \]

and basic properties of variance and the Bernoulli distribution shows

\[ \mathbb{V}\hat{F}_n(x) = \frac{1}{n} F(x) (1 - F(x)). \]

The EDF converges to the CDF in probability (which can be observed from Markov's inequality.) The DKW inequality bounds this convergence which allows for the construction of confidence bounds on the EDF. A functional \(T : F \mapsto \theta\) is called a statistical functional. The plug-in estimator of \(\theta\) is simply \(T(\hat{F}_n)\).

The bootstrap is a simulation based method to estimate standard errors and confidence intervals (there is a standard R package called boot.) Unlike the jackknife the bootstrap has access to the distribution of a statistic rather.

Consider a statistic \(T\) which is a function of a sample of \(X_i\sim F\); we want to compute \(\mathbb{V}_F(T)\). The subscript \(F\) is to indicate that this is with respect to the distribution \(F\). Let \(\hat{F}\) be the empirical distribution function (EDF) of the \(X_i\). The bootstrap uses \(\mathbb{V}_F(T)\) to approximate \(\mathbb{V}_{\hat{F}}(T)\) which is then itself estimated via simulation.

The simulation typically is just sampling from the EDF with replacement, and we can always run more simulation to get an arbitrarily good approximation of \(\mathbb{V}_{\hat{F}}(T)\) (there are \(\binom{n + n - 1}{n}\) potential datasets to sample here, for a dataset of size \(n=140\) that is already the number of atoms in the known universe.) The real source of error is how well \(\hat{F}\) approximate \(F\).

In the parametric bootstrap, rather than sampling from the EDF of the data, a sample is generated from the parametric distribution parameterised by the estimated parameter value.

The profile likelihood is a lower dimensional version of the likelihood function. Consider a likelihood function \(\mathcal{L}(\theta)\) where \(\theta = (\psi,\lambda)\) where \(\lambda\) are nuisance parameters. The profile likelihood is

\[ \mathcal{L}_{\text{profile}}(\psi) := \max_{\lambda} \mathcal{L}(\psi,\lambda). \]

Roughly, it plays the role of a nonparametric paired t-test.

Roughly, it plays the role of a nonparametric t-test.

Roughly, it extends the Mann–Whitney U test to more than two groups.

The Laird-Ware form of a linear mixed effect model (LMM) for the \(j\)th observation in the \(i\)th group of measurements is as follows:

\[ Y_{ij} = \beta_1 + \sum_k \beta_k X_{kij} + \sum_{k} \delta_{ki} Z_{kij} + \epsilon_{ij}. \]

• the \(\beta_k\) are the fixed effect coefficients and the \(X_{kij}\) the fixed effect covariates,
• the \(\delta_k\) are the random effect coefficients and the \(Z_{kij}\) the random effect covariates, it is important to note that while the \(beta_k\) are treated as parameters to be estimated, the \(\delta_k\) are treated as random variables and it is their distribution that is estimated.
• the \(\epsilon_{ij}\) is a random variable.
• the distribution of the random effect coefficients is a (zero-mean) multivariate normal distribution parameterised by \(\psi\) and the random noise \(\epsilon\) comes from a zero-mean normal distribution with variance parameterised by \(\sigma^2\).

Consider a random process, \(X_t\) which is a CTMC on \(\mathbb{N}_0\) where the only possible transitions from \(n\) are to \(n\pm 1\), potentially with an absorbing state at zero. Supposing from state \(n\) the process moves to state \(n+1\) at rate \(a_n\) and to state \(n-1\) at rate \(b_n\) the forward equations for the distribution are

\[ \frac{d}{dt} p_n(t) = p_{n-1}(t)a_{n-1} - p_{n}(t)a_n + p_{n+1}(t)b_{n+1} - p_{n}(t)b_{n}. \]

We multiple both sides of the equation through by \(n\) and sum in \(n\) over \(\mathbb{N}_0\) to get everything in terms of expected values of functions of \(X_t\).

\[ \frac{d}{dt} \mathbb{E}[X_t] = \mathbb{E}[((X_t + 1) - X_t) a_{X_t}] + \mathbb{E}[((X_t - 1) - X_t) b_{X_t}] \]

The "trick" is to note that you can write the first and third sums as \(\sum (n+1) p_{n}(t) a_{n}\) and \(\sum (n-1) p_n b_n\). Similarly, for higher order moments this generalises as

\[ \frac{d}{dt} \mathbb{E}[X_t^k] = \mathbb{E}[((X_t + 1)^k - X_t^k) a_{X_t}] + \mathbb{E}[((X_t - 1)^k - X_t^k) b_{X_t}]. \]

Recall the binomial formula,

\[ (x+y)^k = \sum_{i = 0}^{k} \binom{k}{i} x^i y^{k-i}. \]

This can then be used to derive the differential equations for the moments.

\[ \frac{d}{dt} \mathbb{E}[X_t^k] = \sum_{i = 0}^{k-1} \binom{k}{i} \left( \mathbb{E}[ X_t^i a_{X_t} ] + (-1)^{k-i} \mathbb{E}[ X_t^i b_{X_t} ] \right) \]

Note that the sums are from \(0\) to \(k-1\) but the binomial coefficients use the full \(k\). Given this expression, any polynomial \(a_n\) and \(b_n\) are candidates for the application of a moment closure, which essentially just entails either truncating the cumulants, or assuming a distributional form which expresses higher moments in terms of lower ones.

This refers to the equation

\[ \mathbb{E}g(X) = \int g(x) f_{X}(x) dx. \]

The name derives from the fact that this is so often treated as self-evident rather than being viewed as a theorem.

This result is particularly useful when carrying out MCMC in a transformed parameter space as described in the notes on changing parameterisation above.

Given to ways to measure something, a Bland-Altman plot is a way to compare them.

• Probability and Statistics Cookbook from Matthias Vallentin.

Lagaris et al (1998) describe a method to solve differential equations numerically with neural networks. The simplest example they use is the following differential equation:

\[ \frac{d\Psi}{dx} + \left(x + \frac{1 + 3x^{2}}{1 + x + x^{3}}\right) \Psi = x^{3} + 2x+x^{2}\frac{1 +3x^{2}}{1 + x + x^{3}} \]

with the initial condition that \(\Psi(0) = A\) for \(0\leq x\leq 1\). The solution to this equation is available in closed form:

\[ \Psi = \frac{e^{-x^{2} / 2}}{1 + x + x^{3}} + x^{2}. \]

The trial solution proposed is

\[ \Psi_{t} = A + x N(x,p) \]

where \(N\) is a neural network parameterised by \(p\). This trial solution satisfied the initial condition for all values of \(p\) which leads to an unconstrained optimisation problem. The loss function is

\[ \text{loss}(p) = \sum_{i} \left\{ \frac{d\Psi_{t}}{dx_{i}} + \left(x_{i} + \frac{1 + 3x_{i}^{2}}{1 + x_{i} + x_{i}^{3}}\right) \Psi_{t} - \left[ x_{i}^{3} + 2x_{i}+x_{i}^{2}\frac{1 +3x_{i}^{2}}{1 + x_{i} + x_{i}^{3}} \right] \right\}^{2} \]

for some set of trial points \(x_{i}\). There are closed forms available for both the loss function and its gradient which opens the way for multiple optimisation routines to be applied.

Figure 4 shows the result of less than a minute of training with gradient descent with only ten test points.

The script that generated the image in Figure 4 is shown below.

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1)

H = 10
mesh_size = 10
# We need to select some points to train the model on and some to test it on.
# It looks like you need to have this as a mutable value if you want to compute
# gradients.
train_mesh = tf.Variable(np.linspace(0, 1, mesh_size), dtype=np.float64)
test_mesh = tf.Variable(np.linspace(0.01, 0.99, mesh_size), dtype=np.float64)
# We need to select some points to plot the resulting approximation on.
plt_mesh = np.linspace(0, 2, 100)
# We need to record the solution at multiple points of the training loop to
# test that it is infact converging so something sensible.
psi_val_stash = []
iter_stash = []
loss_stash = []
vldt_stash = []

@tf.function
def psi_truth(x):
    """Solution of the differential equation."""
    return tf.exp(-0.5 * x**2) / (1 + x + x**3) + x**2


class MyModel(tf.Module):
    """This class respresents the approximation to the solution."""

    def __init__(self, psi_0, num_units, **kwargs):
        """Construct approximate solution with a hard-coded IC."""
        super().__init__(**kwargs)
        self.IC_A = psi_0
        self.w = tf.Variable(np.random.randn(num_units))
        self.u = tf.Variable(np.random.randn(num_units))
        self.v = tf.Variable(np.random.randn(num_units))

    def __call__(self, x):
        """Evaluate the approximation."""
        # matrix with one datum per row
        l1 = tf.math.sigmoid(tf.tensordot(x, self.w, 0) + self.u)
        # We have encoded the IC using the method from Lagaris et al (1998),
        # although there are other formulations that are used in contemporary
        # implementations. Attempting to use a 1-exp(-x) function to get the IC
        # constraint did not yield any improvement.
        l2 = self.IC_A + x * tf.tensordot(l1, self.v, 1)
        return l2


psi_model = MyModel(tf.constant(1.0, dtype=np.float64), H)

num_iters = 2000
learning_rate = 0.050
x = train_mesh


def loss_fn(psi_dash, psi, x):
    return tf.reduce_sum((psi_trial_dash + (x + (1 + 3 * x**2) / (1 + x + x**3)) * psi_trial - (x**3 + 2 * x + x**2 * (1 + 3 * x**2) / (1 + x + x**3)))**2)


for iter_num in range(num_iters):
    # The gradient evaluations need to be within to the loop for TF to
    # understand how they work.
    with tf.GradientTape() as t2:
        with tf.GradientTape() as t1:
            psi_trial = psi_model(x)
        psi_trial_dash = t1.gradient(psi_trial, x)
        # The loss function needs to know the differential equation.
        loss = loss_fn(psi_trial_dash, psi_trial, x)
        vldt_val = loss_fn(psi_trial_dash, psi_trial, test_mesh)
    loss_w_dash, loss_u_dash, loss_v_dash = t2.gradient(loss, [psi_model.w, psi_model.u, psi_model.v])
    psi_val_stash.append(psi_model(plt_mesh).numpy())
    psi_model.w.assign_sub(learning_rate * loss_w_dash)
    psi_model.u.assign_sub(learning_rate * loss_u_dash)
    psi_model.v.assign_sub(learning_rate * loss_v_dash)
    if iter_num % 10 == 0:
        iter_stash.append(iter_num)
        loss_stash.append(loss.numpy())
        vldt_stash.append(vldt_val.numpy())
        print("iteration: {i}\tloss: {l}\tlearning rate: {r}".format(i=iter_num, l=loss.numpy(), r=learning_rate))
    if iter_num % 1000 == 0 and iter_num > 0:
        learning_rate /= 1.5

fig, axs = plt.subplots(1, 2, figsize=(10, 5))
axs[0].plot(plt_mesh, psi_truth(plt_mesh), label="True solution")
for ix in range(1, num_iters, num_iters // 3):
    axs[0].plot(plt_mesh, psi_val_stash[ix], label="Iteration {n}".format(n=ix+1))
axs[0].scatter(train_mesh, psi_truth(train_mesh), label="Training data")
axs[0].legend(loc="upper left")
axs[0].set_title("Approximate solutions of ODE")

axs[1].plot(iter_stash, loss_stash, label="Training")
axs[1].plot(iter_stash, vldt_stash, label="Validation")
axs[1].set_title("Loss functions")
axs[1].set_xlabel("Iteration")
axs[1].legend(loc="upper right")
# fig.show()
fig.savefig("out/lagaris-1998.png", dpi=300)

• When initialising the weights of a neural network, keep in mind that we want the linear map of the input to take a value that has values roughly compatible with a standard Gaussian. As a result, if you have \(n\) (standardised) input variables, you might choose to initialise the weights from \(N(0,1\sqrt{n})\).

https://github.com/ankitrohatgi/WebPlotDigitizer

This tool can be used to extra data from a plot in a figure in paper! If you want to extract a figure from a PDF of a paper, see my Inkscape notes here.