Computational Statistics 计算统计

0 What is Statistics?

There exists two main philosophic ideas for Statistics: Frequentist, and Bayesian. Frequentists assume there is a fixed, exact value for a parameter (same as classical mechanics); Bayesian assume the "value" is jumping, having a probability everywhere, and doesn't have a fixed value (same as quantum mechanics). The real world tells us, that generally, Bayesian is better than Frequentist.

More interestingly, Frequentists believe in "god", but Bayesian believes there's no "god".

What is Statistical Inference? samples -> parameters + variance. From a higher view, variance is a more important idea than mean (any statistics can have a variance).

Three types of statistical inference classification:

• Frequentists & Bayesian

• Parametric & Non-parametric

• Parameter Estimation & Hypothesis Testing

1 Formula

1.1 Expectation

An expectation is defined as the weighted sum of the whole variable inside the operator \(\mathbb E[\cdot]\).

The first law is

\[\mathbb E_{f(X)\sim \phi(X)}[g(f(X))]\equiv \mathbb E_{X\sim\psi(X)} [g(f(X))], \]

or

\[\mathbb E[g(X)]=\mathbb E_X[g(X)], \]

which is called Law of the unconscious statistician (LOTUS). Oftenly use this to eliminate the function in the subscript.

My favorite example is

\[\mathbb E_{X^2}[X^2]=\mathbb E_{X}[X^2]. \]

Another example of \(\chi^2_k\) generating function \(M_Y(t)\) is

\[\mathbb E_Y[\text e^{tY}]=\left(\mathbb E_{X_i^2}[\text e^{tX_i^2}]\right)^k=\left(\mathbb E_{X_i}[\text e^{tX_i^2}]\right)^k. \]

1.2 Moment Generating Function (MGF)

A generating function is a function, whose coefficients of series/intergral expansion is the sequence, by introducing a new variable. It defines a probability.

For example, a moment generating function is in the form of

\[M_X(t)=\mathbb E_X[\text e^{tX}]=\text{oftenly a closed form function} \]

and by definition, this will yield \(M'_X(0)=\mathbb E[X]\), \(M''_X(0)=\mathbb E[X^2]\), and \(M^{(n)}_X(0)=\mathbb E[X^n]\) so on.

1.3 Bootstrap

Bootstrap is a method to estimate statistics, when only provided a fixed set. Note the real world distribution is not a fixed set.

The core is believe the bootstrap sample represents real data.

Take CLT for example, the \(\bar{X}_n\) means sampling \(n\) times from real world distribution (the set is changing). This is different from sampling a bootstrap set \(\{x_1,\dots,x_n\}\) (the set is fixed).

The first \(\hat \theta=\bar X\), the second \(\bar\theta^*=\dfrac1B\sum \bar X^{*(b)}\). The bootstrap bias is defined as \(\hat{\text{bias}}= \bar{\theta}^*-\hat\theta\).

After \(B\) resamplings, choose percentile from these \(\hat\theta^{*(i)}\)'s.

Bootstrap residuals and pairs are all non-parametric (sometimes we call residual half-parametric, or just non-parametric). The criteria is: do we assume the model? And the distribution of residuals? If all yes, then parametric.

It can be shown that \(\text{Var}(\hat{\text{se}}_N)\approx \dfrac{c_1}{n^2}+\dfrac{c_2}{nN}\). \(N\) is boostrap times/iterations, \(n\) is sample size.

1.4 Fisher Information

Fisher Information measures how much information related to \(\theta\) is in the samples. A higher Fisher Information yields a preciser estimate of the parameter \(\theta\). We notate \(\mathcal I(\theta)\) as a single sample.

Calc: Var of Score Function / neg-expectation Hessian.

1.5 Standard Error

Standard errors \(\hat{s.e.}\) are generalised standard deviation \(\sigma\) (let the statistic be the variable itself).

1.6 GLM

Def: Randomness Distribution + Mapping Function. Eg: Logistic Regression / Poisson Regression. GLM is always in expo-family.

1.7 Expo-family

All distributions can be categorized into two catrgories: expo-family, and others.

\[f(x\mid \theta)=h(x)\exp(\eta(\theta)T(x)-A(\theta)). \]

Canonical sufficient estimates \(T(x)\) only come from expo-family.

1.8 Projection Matrix

For OLS, we have \(\hat y=Py\), and residual vector \(e=y-\hat y=(I-P)y\) which is \(n-p\) dimension.

Variance is defined as mean fluctuation of every direction of \(e\). So the estimated denominator should substract \(p\) dimensions.

And note that \(n\) only occurs in estimating.

1.9 Asymptotic Theory

CLT -> Normal(mean)
EVT -> GEV(extreme)

2 Stability

2.1 Median

For discrete median functionals, there exists two problems: non-uniqueness, and "jump".

A median is defined as \(P(X\leq m)\geq0.5, P(X\geq m)\geq 0.5\). Non-uniqueness means there are multiple medians exist; "jump" means there's no exact \(F^{-1}\left(\dfrac12\right)\) point.

2.2 Confidence Interval

A useful CI doesn't only mean an accurate "coverage". The CI should have information in data itself, and short when satisfying coverage, and nominal coverage close to actual coverage, robustness to assumptions.

A typical example is Uniform Distribution: you can't sample the max \(X_i\) so the actual CI is shorter than nomial CI; and there's high prob to bootstrap sample the highest value \(X_i'\) so will yield a fake spike.

General formula of CI's (Wald's style):

\[[\text{estimate}\pm \text{quantile}\times\hat{\text{SE}}] \]

But not all of them are symmetric.

Note: for monotonic \(g(\lambda)\)'s, computing 95%-CI for \(\lambda\in[a,b]\) and \(g(\lambda)\in[g(a),g(b)]\) are the same.

2.3 Independent Sample Generator

For easy CDF's, use Inverse Function method. Otherwise, use Rejection Sampling.

For a real \(f(x)\), we first try a proposal \(g(x)\), if:

\[\frac{f(x)}{Mg(x)}\geq U, \quad U\sim\text{Unif}(0,1)\text{ and }M=\sup\frac{f(x)}{g(x)} \]

accept; otherwise reject. \(\dfrac1M\) is acceptance rate. Note that \(M\geq1\). We also define acceptance rate at each \(Y\) as \(h(Y)=\dfrac{f(Y)}{Mg(Y)}\in(0,1)\). \(Y\)'s come from \(g\).

Furthermore, treat a \(Y\in(y,y+\delta y)\), then \(P(Y\in(y,y+\delta y))=g(y)\delta y+o(\delta y)\), and \(P(A\mid Y\in(y,y+\delta y))=h(y)+o(1)\), multiply together \(P(Y\in(y,y+\delta y),A)=g(y)h(y)\delta y+o(\delta y)\). This is the prob that \(Y\) is approximately \(y\) and accepted. Therefore, the density of those accepted \(Y\)'s satisfy \(\propto g(y)h(y)\). Don't forget to normalize by dividing the whole acceptance rate.

Note: if we only know a multiplicative constant, the rejection method still works.

Good envelopes: small \(M\) and easy \(g\).

Dependent Sample Generator refers to Markov Chain Monte Carlo/MCMC.

2.4 Profile Likelihood

Only for frequentist method (Bayesian use an integral).

For nuisance parameters, treat their values as "already maximized" likelihood. This is approximately a nested likelihood. Its a Upper Envelope.

2.5 Wilks Theorem

Used for LRT. Know how to prove it (Taylor + Cholesky Decom.)

3 Simulation

3.1 MC

The biggest advantage of simulation study is: we know how the data is generated! (which means we know the true value of parameters, so that to estimate errors)

Common aims: efficiency (seen in rejection-sampling's \(M\)), precision, large/small sample bias, Standard Error, CI coverage (also seen in rejection-sampling's \(M\)), robustness to misspecification/wrong assumption.

MCSE formula:

\[\text{MCSE}=\sqrt{\frac{\text{Var}(T)}{N}} \]

or \(\displaystyle{\sqrt{\frac{\hat p(1-\hat p)}{N}}}\).

If we don't know \(\text{Var}(T)\) then substitute with sample variance. We often use \(3\times \text{MCSE error}\).

MC error is regardless of dimension (always \(\mathcal O(n^{-\frac12})\)) (even higher).

3.2 LLN

WLLN (in probability, the prob of sequence converges) and SLLN (almost surely, the sequence itself converges).

3.3 Classical Inference

For normals, use t-test (don't use z-test); otherwise, bootstrap (non-parametric).

If sample numbers too small, bootstrap (very sensitive to lack of outliers) will underestimate the variance (so the nominal coverage is higher than actual); for t-test (very sensitive to existence of outliers) will yield a high actual coverage (overestimate variance by outliers).

Be aware of time-series.

3.4 EM

EM is for those MLE with some latent variables. GMM, HMM, K-means are all EM algorithms.

ZIP is a perfect idea to understand EM. Think of \(\rho\) as a prior, indicator function as likelihood.

3.5 Tests

z-test is also called Wald test. It's not better than LRT (likelihood ratio test, supported by Wilks' Theorem).

t-test is also called Studentized t-test.

3.6 General Multivariate Strategy

First draw \(Y_1\), then do \(Y_k\mid Y_{1:k-1}=y_{1:k-1}\) in order. For example a bivariate normal, first do \(X\sim\mathcal N(0,1)\), then do \(Y\mid X\sim\mathcal N(\rho x,1-\rho^2)\).

3.7 Conjugate Priors

A flat prior is often improper(very wide/doesn't converge).

Binomial <-> Beta
Exponential <-> Gamma
Poisson <-> Gamma
Multinomial <-> Dirichlet
Normal(\(\mu\)) <-> Normal
Normal(\(\sigma^2\)) <-> Inverse Gamma (Or we often use precision \(\tau=\sigma^{-2}\))
Uniform <-> Pareto

3.8 Random Effect Models

\(Y_i=\mu+\alpha_i+\varepsilon_i\), \(\varepsilon\sim\mathcal N(0,\sigma_1^2)\), \(\alpha_i\sim \mathcal N(0,\sigma^2)\),

Here we have 16 observations, so we have 19 parameters together (\(\mu, \{\alpha\}_i, \sigma_1,\sigma^2\)).

More generally,

\[Y_{ij}\sim \mathcal N(\mu+\alpha_i,\sigma_1^2),\quad \alpha_i\sim\mathcal N(0,\sigma_2^2),\quad \mu\sim\mathcal N(\mu_0,\tau^2\sigma_1^2),\quad \dfrac{1}{\sigma_1^2}\sim\text{Gamma}(a_1,b_1),\quad \dfrac{1}{\sigma_2^2}\sim\text{Gamma}(a_2,b_2) \]

Prior Scaling: We define the variance of \(\mu\) scales to \(Y\)'s. (to avoid a \(\sigma_3\) which is hard to compute)

4 R

4.1 nlm and nlminb

They are all for non-linear minimization. nlminb is always better (quasi-Newton). glm's are based on lower architectures like nlm functions.

4.2 boot

Common function is sample(x, length(x), replace = TRUE). For example:

stat_fun <- function(x) {
  mean(x)
}

B <- 1000
theta_star <- numeric(B)
for (b in 1:B) {
  x_star <- sample(x, length(x), replace = TRUE)
  theta_star[b] <- stat_fun(x_star)
}

We can also use boot library:

my_mean_func <- function(data, indices) {
  resampled_data <- data[indices]
  return(mean(resampled_data))
}

my_data <- rnorm(100, mean = 50, sd = 10)
boot_results <- boot(data = my_data, statistic = my_mean_func, R = 1000)

Though in the course we wrote:

bootstrap <-function(x,nboot,theta,...)
{
  data <- matrix(sample(x, size = length(x) * nboot, replace = T ), nrow = nboot)
  answer <- apply(data, 1, theta, ...)
  answer
}

Apply a function to the dimension designated: apply(data, 1, theta, ...). 1 stands for the dimension 1 (row), means the row is the variable here, so the theta is calculated for each row.

Same methodology for cbind, rbind, where cbind adds column numbers, otherwise row numbers. (rows represent samples, columns represent feature dimensions.)

4.3 Monte Carlo Experiments

For coverage:

cover <- numeric(N)
for (r in 1:N) {
  dat <- rnorm(n)
  ci <- t.test(dat)$conf.int
  cover[r] <- as.integer(ci[1] <= 0 && 0 <= ci[2])
}
p_hat <- mean(cover)
mcse <- sqrt(p_hat * (1 - p_hat) / N)

The $ symbol is the extraction operator.

For a Beta(0.5, 2) rejection-sampling example:

sim.fn2 <- function(a = 0.5, b = 2, nvals = 1000) {
  X.accepted <- NULL; i <- 0
  repeat {
    Y <- runif(1)^(1/a)
    X <- runif(1, 0, Y^(a - 1))
    if (X <= Y^(a - 1) * (1 - Y)^(b - 1)) {
      i <- i + 1; X.accepted[i] <- Y
    }
    if (i >= nvals) break
  }
  X.accepted
}

Repeat until nvals sample generated.

For a Bivariate Normal Simulation (MC Variability):

sim.bvn <- function(nvals = 1000, rho = 0.25)
{
  X <- rnorm(nvals)
  Y <- rnorm(nvals, mean = rho * X, sd = sqrt(1 - rho^2))
  cbind(X, Y)
}

monte.carlo.fn1 <- function(nvals = 1000, rho = 0.25)
{
  simvals <- sim.bvn(nvals, rho)
  X <- simvals[,1]
  Y <- simvals[,2]
  mean(X^2 + 2 * Y^3 > 3)
}

monte.carlo.fn2 <- function(nreps = 100, nvals, rho)
{
  est <- numeric(nreps)
  for (i in 1:nreps)
  {
    est[i] <- monte.carlo.fn1(nvals, rho)
  }
  est
}

simvals <- monte.carlo.fn2(100, 25000, 0.25)
hist(simvals)

The sim.bvn here shows how to generate bivariate normals with rho.

Note if we want to evaluate the variability of a "mean", we need to repeat many "means". This is a triple loop.

4.4 Profile Lik

plik.plot <- function(range)
{
  n <- length(range)
  plik <- numeric(n)
  for(i in 1:n){
    psi <<- range[i]
    plik[i] <- -1 * (nlm(norm.loglik.chi, 1)$minimum)
  }
  plot(range, plik, type=’l’, xlab="psi", ylab="log-likelihood")
}

Here the <<- is assigning the variable outside.

Plot: The first two parameters are vectors (x and y), then type, and xlab, ylab.

4.5 JAGS

model {
  for (i in 1:N) {
    X[i] ~ dnorm(theta, tau)
  }
  theta ~ dnorm(1, 0.001)
  sigma2 <- 1.0/sqrt(tau)
  tau ~ dgamma(0.001, 0.001)
}

5 Statistical Functionals

5.1 Robust Statistics

In frequentists, we assume statistic a functional of CDF's. This is not Bayesian because for a given \(F\), the \(p(\theta\mid F)\) isn't fixed (also depends on prior \(p(\theta)\)).

For mean \(t(F) = \mathbb E_F[X]\),

\[\begin{aligned} t(\tilde F)=&\int s\text d \tilde F(s)\\ =&\int s((1-\varepsilon)\text d F+\varepsilon s\delta_x(s)\text ds)\\ =&(1-\varepsilon)\int s\text dF+\varepsilon x\int_x^x \delta_x(s) \text ds\\ =&(1-\varepsilon)\int s\text dF+\varepsilon x \end{aligned} \]

Hence

\[\frac{\delta t}{\delta F}=x-\int s\text dF=x-\mathbb E_F[X]. \]

Mean is Frechet differentiable.

For median \(t(F) = \text{median}(F)\),

Denote

\[t(\tilde F)=\text{median}(\tilde F)=s(\varepsilon)=\tilde F^{-1}(\frac12) \]

\[t(F)=\text{median}(F)=s(0)=F^{-1}(\frac12) \]

Hence

\[\lim_{\varepsilon\rightarrow 0} \frac{t(\tilde F)-t(F)}{\varepsilon}=\lim_{\varepsilon\rightarrow 0} \frac{s(\varepsilon)-s(0)}{\varepsilon}=\frac{\partial s}{\partial \varepsilon}\bigg|_{\varepsilon=0} \]

Expand \(\tilde F^{-1}\) and take derivative of \(\varepsilon\),

\[\begin{aligned} (1-\varepsilon)F(s(\varepsilon))+\varepsilon \Delta_x(s(\varepsilon))=&\frac12\\ \frac{\partial s}{\partial \varepsilon}\bigg|_{\varepsilon=0}=&\frac{\text{sgn}(x-s(0))}{2\frac{\partial F}{\partial s}\big|_{s=s(0)}}\\ \frac{\delta t}{\delta F}=&\frac{\text{sgn}(x-F^{-1}(\frac12))}{2f(F^{-1}(\frac12))}. \end{aligned} \]

Or，we use \(\dfrac{\delta t}{\delta F}\) in statistics, or \(U(x,F)\) /Score.

The variance is functional Fisher Information.

Median is a typical non-Frechet differentiable example (not satisfying linear sum property).

6 Random Process

We learn these as Bayesian Approach. Why do we do sampling? Because we can easily know a value for each point(P), but don't know the closed formula for the entire distribution(NP).

6.1 Markov Chain Monte Carlo

MCMC is a sampling method, for finding prob dist. MCMC is P algorithm.

Markov property is the main idea. (Transient Analysis/C-K eq + Asymptotic Analysis/Stationary eq)

Time Homogeneity: \(\Pr(X_n = j|X_0 = i) = \Pr(X_{n+k} = j|X_k = i)\). Based on this, we have two key defs: Stationary/Equilibrium Distribution and (Detailed) Balance. NOTE: Stationary is NOT Convergence!

\[\pi P=\pi\quad\text{or}\quad \pi=\lim_{n\rightarrow\infty} \bar v P^n \]

\[\pi_i P_{ij}=\pi_j P_{ji} \]

Solving stationary is the same as computing eigenvalues. We denote sub as states, sup as timestep(which isn't displayed).

Solving detailed balance is the same as reversibility.

We often note \(\bar \pi,\pi^{(0)}\) as the initial dist, and \(\pi\) is the stationary dist.

More properties:

• A \(\lambda=1\) ensures stationary.

• A \(|\lambda|<1\) ensures convergence.

• If \(\lambda<0\) then bounce, other wise monotone.

For continuous states:

\[\pi(y)=\int \pi(x)k(y\mid x)\text dx \]

We focus on discrete time Markov chains more (discrete or continuous state). The matrix here is ALWAYS square matrix.

\[P(X_{m+n+1}=j\mid X_m=i)=\sum_k P(X_{m+n+1}=j\mid X_{m+n}=k)=\sum_k P_{kj}(P^n)_{ik}=(P^{n+1})_{ij}. \]

This is called Chapman-Kolmogorov Equation(C-K Equation).

How to think an MCMC process?

timestep -> state -> components of a single state

burning & thinning

6.2 Metropolis-Hastings (JAGS)

JAGS used for discrete.

\[\alpha=\min (1,\dfrac{P(x')q(x\mid x')}{P(x)q(x'\mid x)}) \]

the \(\dfrac{q(x\mid x')}{q(x'\mid x)}\) is called Hastings ratio/correction factor (eliminates bias, for example punishing a high \(x'\) and rewards a low 'x').

then compare runif \(u\) and \(\alpha\).

So \(P(x,y)=\alpha\times q(y\mid x)\). (Metropolis Rule)

• IF set \(q\) is symmetric, this is Metropolis algorithm.

• Gibbs Sampler is the simplest MH (\(q(x'\mid x)=p(x_i\mid x_{-i})\)), only choose one each turn and set others constant. Never reject, but sometimes slow.

• Simulated Annealing is using MCMC. (the craziest part is it doesn't require computing full \(P\), note that Bayesian's main problem is the denominator computing)

Problems with tuning parameters: too high -> hard to accept; too low -> step too small.

6.3 Hamiltonian Mechanics (STAN)

STAN used for continuous.

For huge datasets, posterior acts like probability in frequentists. (Bernstein-von Mises Theorem, Bayesian CLT)

If regular: \(N\rightarrow\infty\), converge to Normal, \(\mu\) defined by MLE(so sometimes can use flat prior), \(\sigma^2\) defined by Fisher Information(divided by datasize).

• What is Regular?

Identifiable, inside parameter space, proper priors etc.

Otherwise, multiple modes etc.

Regularize = Prior = Punishment

6.4 Spatial Point Processes

Poisson Process: Every point is independent from each other. (looks not random)

Note that micro-settings they are not random, but in macro-settings the density is converged.

Strauss Process (a kind of Gibbs Process): Close points will be punished. (looks random)

7 Sampling and Optimizing Summary

7.1 Sampling

Want a probability.

• Classic: Inverse, Rejection-sampling.

• MCMC: Gibbs, Metropolis-Hasting, STAN etc.

7.2 Optimizing

Want a best value.

• 0-order: Without derivative information, eg. Hill Climbing, Simulated Annealing, Genetic Algorithms, Random Search etc.

• 1-order: All kinds or variations of GD (PGD, SGD, CGD, Adam etc).

• 2-order: Newton or Quasi-Newton Methods.

7.3 Variational Inference (Sampling+Optimizing)

Use a computational-friendly dist to approximate original dist. (ELBO etc)