Statistical modeling¶
Variation and uncertainty¶
- Mathematical modeling
- Interpretations of probability:Frequentist and Bayesian.
Spatial variation¶
- Regression: input vs feature space
- Covariance functions
- Gaussian process regression
- Parametric and non parametric models
- Geostatistics
Bayesian Computing¶
- Sampling
- Other methods
- Maximum likelihood type II
- Laplace approximation
- INLA
The examples are written in Python 2.7 (I'm sorry R users). These are just some packages that we will be using.
from IPython.display import YouTubeVideo, Image
from sklearn import linear_model
from scipy import stats
import numpy as np
import pylab as plt
import GPy
%matplotlib inline
# Plot parameters
plt.rcParams["figure.figsize"] = [12,4]
... but first youtube videos!¶
YouTubeVideo('adKCHpdY6ZM')
YouTubeVideo('Ivo1uUt9QlE')
Variation and Uncertainty¶
How do we catch a moving object in the air?
We do not usually stop to think about it, we simply do. But behind this task, there is a cognitive process to extract the pattern about the trajectory and speed of the object. In some way, it is like the object is "telling" us where and when it is moving next... but it is not the object that is telling this, it is our $model$ of the object (which by the way is being updated in real time!).
Unlike the cat, the toddler in the video lacks of experience. He understands what he needs to do, but he hasn't managed to understand where is the best place to be standing to wait for the ball or when to start swinging the bat. We can think of these as the $parameters$ of his model. The boy has the right model, he just hasn't calibrated the parameters. The calibration of the model parameters will come with experience. In other words, it will come with more data. The more attempts he carries on, the more repeated examples he will have available and better parameters he will use.
That is pretty much statistical modeling. We formulate an abstract representation of some aspect of our reality. Then we calibrate this representation using the available data and we use it. Afterwards, we might go back an re-adjust our model with what we learnt by using it.
Our interests usually go beyond catching birds of hitting home-runs. We build models for aspects of reality variables more difficult to understand, or even to observe. Think of another example: What kind of clothes should we pack on a travel? Is it going to be cold? is it going to rain? We usually don't want to carry too many things on a travel, but we want to carry enough to manage the different weather conditions we might face. We know weather changes (i.e., it has $variation$)... yes, in some places more than others, and during some times of the year more than others, and that information can be part of our model as well. Yet the process by which all the variables that define the weather conditions for the specific dates of our trip is unintelligible for us. This ignorance of the mechanics of the weather is the source of our $uncertainty$. And when things are uncertain, we need $probability$.
Across the web we find consistent definitions of probability:
- "The chance that a particular event (or set of events) will occur expressed on a linear scale from 0 (impossibility) to 1 (certainty)"
http://mathworld.wolfram.com/Probability.html
- "A probability is just the chance that some event will happen."
https://www.brainpop.com/math/probability/basicprobability/
- "Probability is the measure of the likelihood that an event will occur."
https://en.wikipedia.org/wiki/Probability
- "The extent to which an event is likely to occur, measured by the ratio of the favorable cases to the whole number of cases possible." https://en.oxforddictionaries.com/definition/probability
Yet there are different interpretations about what such definitions imply. Here we will review the main differences between the $Frequentist$ and the $Bayesian$ approaches.
Here are a couple of good blogs about these approaches:
http://jakevdp.github.io/blog/2014/03/11/frequentism-and-bayesianism-a-practical-intro/
http://allendowney.blogspot.com/2016/09/bayess-theorem-is-not-optional.html
Frequentist approach¶
Probabilities are associated to frequencies of events:
$P(A) = \frac{Cases\; involving\; A}{Total\;number\;of\; cases}$¶
This means that we need to define our model in terms of a sample space (whose size define the total number of cases).
Either we know the sample space and the cardinality of the subsets of interest or at least we do many repetitions.
Inference:
- Data $y$ are random
- Model parameters $\theta$ are fixed.
- Inference is performed via maximum likelihood $L(\theta|y) = p(y|\theta)$. This means finding $\theta$ such that $p(y|\theta)$ is maximized.
- We cannot make probability statements about the parameters.
- We estimate $\theta$ and define $confidence\; intervals$ for its value.
Things can get confusing in the last two bullets. A 95% confidence interval doesn't really mean that there is a 95% chance of the value of theta being in the interval, it means that if we could repeat the sampling and estimation many times, 95% of the confidence intervals we define would contain the value of $\theta$.
Image(filename='figs/confidence_intervals.png', width = 500)
Bayesian approach¶
Probabilities represent degrees of belief. No need to think in terms of the sample space.
Inference:
- Data $y$ are fixed
- Model parameters $\theta$ have a probability distribution
Inference is performed via the posterior distribution $p(\theta|y)$: ## $p(\theta|y) = \frac{p(y|\theta) p(\theta)}{p(y)}$
We can incorporate prior knowledge/belief about $\theta$, through $p(\theta)$.
- We do not obtain point estimates of $\theta$, but probability distributions of its value.
- We can make probability statements about $\theta$.
- We define $credible\; intervals$.
A credible interval actually means the chances of $\theta$ being within certain range of values.
# Say the coin is not fair ant 30% of cases it falls in tail
theta_tail = .3
# Let's flip a coint 1000 times
num_trials = 15
trial = np.arange(1,1+num_trials)
bernoulli_output = stats.binom.rvs(1, p = theta_tail, size = num_trials)
# Now compute the moments method estimator at each new trial
ac_tails = np.cumsum(bernoulli_output)
ac_trials = trial.astype(float)
mme = ac_tails/ac_trials
ciuf = mme + 1.96 * np.sqrt(mme*(1-mme)/ac_trials)
cilf = mme - 1.96 * np.sqrt(mme*(1-mme)/ac_trials)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hlines(y = .3, xmin=1, xmax=num_trials, colors = 'lightgreen', lw = 2, label = 'actual value')
ax.plot(trial, mme, label = 'frequentist estimator', lw = 2)
#ax.fill_between(trial, cilf, ciuf, facecolor = 'b', alpha = .2, label = '95% confidence interval (proxy)')
ax.legend()
ax.set_ylim(-.2,1.5)
ax.set_xlim(1, num_trials)
Bayesian approach¶
Say we initially think there is no bias in the coin, so that $P(tail) = .5$
We can still express our confidence about this belief in the model. We can do this with a Beta distribution:
$p(\theta) = \frac{\Gamma(a + b)}{\Gamma(a)\Gamma(b)} \theta^{a-1} (1-\theta)^{b-1}$¶
[You can review the theory here: https://www.cs.cmu.edu/~10701/lecture/technote2_betabinomial.pdf]
_x = np.linspace(1,999)/1000
prior_a = 5
prior_b = 2
plt.plot(_x,stats.beta.pdf(_x, prior_a, prior_b))
With the results above we can compute the prior
# Posterior distribution of theta
posterior_a = prior_a + ac_tails
posterior_b = prior_b + ac_trials - ac_tails
expected_theta = posterior = posterior_a/(posterior_a + posterior_b).astype(float)
ciub = np.array([stats.beta.ppf(.975, a_i, b_i) for a_i,b_i in zip(posterior_a, posterior_b)])
cilb = np.array([stats.beta.ppf(.025, a_i, b_i) for a_i,b_i in zip(posterior_a, posterior_b)])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hlines(y = .3, xmin=0, xmax=num_trials, colors = 'lightgreen', lw = 2, label = 'actual value')
ax.plot(trial, mme, label = 'frequentist estimator', lw = 2)
#ax.fill_between(trial, cilf, ciuf, facecolor = 'b', alpha = .2, label = '95% confidence interval (proxy)')
ax.plot(trial, expected_theta, 'r-', lw = 2, label = 'bayesian estimator')
ax.fill_between(trial, cilb, ciub, facecolor = 'r', alpha = .2, label = '95% credible interval')
ax.legend()
ax.set_ylim(-.2,1.5)
ax.set_xlim(1, num_trials)
There is a balance between prior beliefs and evidence.
new_trials = 5
trial = np.arange(1,1+num_trials)
new_prior_a = 2
new_prior_b = 5
new_tails = stats.binom.rvs(new_trials, p = theta_tail, size=1)
new_posterior_a = new_prior_a + new_tails
new_posterior_b = new_prior_b + new_trials - new_tails
_x = np.linspace(.01,.99,100)
prior_dist = stats.beta.pdf(_x, new_prior_a, new_prior_b)
posterior_dist = stats.beta.pdf(_x, new_posterior_a, new_posterior_b)
likelihood = np.array([stats.binom.pmf(new_tails, p= theta_i,n = new_trials) for theta_i in _x])
scale_factor = np.hstack([prior_dist, posterior_dist]).max()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(_x, prior_dist, 'b', label = 'prior belief')
ax.plot(_x, likelihood * scale_factor/likelihood.max(), 'green', label = 'evidence', lw=2)
ax.plot(_x, posterior_dist, 'r', label = 'posterior')
ax.set_xlim(0,1)
ax.legend()
Example 2: linear regression¶
Frequentist approach¶
Linear regression formulation (in matrix notation)
$\textbf{y} = \textbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}$,
with $\textbf{y} \in \mathbb{R}^N$, $\textbf{X} \in \mathbb{R}^{N \times Q}$, $\boldsymbol{\beta} \in \mathbb{R}^Q$ and $\boldsymbol{\epsilon} \in \mathbb{R}^N$. Only the noise is random: $\boldsymbol{\epsilon} \sim N(\textbf{0}, \textbf{I}\sigma^2)$.
We want to minimize the sum of squares: $\boldsymbol{\epsilon}^\top\boldsymbol{\epsilon} = (\textbf{y} - \textbf{X} \boldsymbol{\beta})^\top (\textbf{y} - \textbf{X} \boldsymbol{\beta})$
so we optimize wrt. $\boldsymbol{\beta}$ and we end up with:
$\hat{\boldsymbol{\beta}} = (\textbf{X}^\top \textbf{X})^{-1} \textbf{X}^\top \textbf{Y}$
$\boldsymbol{\beta}$ is not a random variable, but $\hat{\boldsymbol{\beta}}$ is.
Expected value:¶
$\langle \hat{\boldsymbol{\beta}} \rangle = \langle (\textbf{X}^\top \textbf{X})^{-1} \textbf{X}^\top \textbf{Y} \rangle$
$\langle \hat{\boldsymbol{\beta}} \rangle = \langle (\textbf{X}^\top \textbf{X})^{-1} \textbf{X}^\top (\textbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}) \rangle$
$\langle \hat{\boldsymbol{\beta}} \rangle = \boldsymbol{\beta}$
variance:¶
$var( \hat{\boldsymbol{\beta}} ) = var( (\textbf{X}^\top \textbf{X})^{-1} \textbf{X}^\top \textbf{Y})$
$var( \hat{\boldsymbol{\beta}} ) = var( (\textbf{X}^\top \textbf{X})^{-1} \textbf{X}^\top (\textbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}) )$
$var( \hat{\boldsymbol{\beta}} ) = var( (\textbf{X}^\top \textbf{X})^{-1} \textbf{X}^\top \boldsymbol{\epsilon} )$
$var( \hat{\boldsymbol{\beta}} ) = \boldsymbol{\sigma^2} (\textbf{X}^\top \textbf{X})^{-1}$
Confidence intervals are constructed based on the distribution of $\hat{\boldsymbol{\beta}} $.
Predictions of $\textbf{y}$¶
Predictions of $\textbf{y}$ are made using
$\textbf{y} = \textbf{X} \hat{\boldsymbol{\beta}}$.
Bayesian approach¶
[You can check the maths here, in Chapter 2: http://www.gaussianprocess.org/gpml/ (Appendix A is also handy)]
We start with
$\textbf{y} = \textbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}$,
but now we assume the parameters have a distribution, lets say $\boldsymbol{\beta} \sim N(\textbf{0}, (\boldsymbol{\Sigma}_b)$
$ p(\boldsymbol{\beta}|\textbf{y}) \propto p(\textbf{y}|\boldsymbol{\beta}) p(\boldsymbol{\beta})$
$ p(\boldsymbol{\beta}|\textbf{y}) \propto N(\textbf{X}\boldsymbol{\beta}, \textbf{I}\sigma^2) N(\textbf{0},\boldsymbol{\Sigma}_b)$
$ p(\boldsymbol{\beta}|\textbf{y}) \propto \exp(-\frac{1}{2}(\textbf{y} - \textbf{X}\boldsymbol{\beta})^\top \textbf{I} \sigma^{-2} (\textbf{y} - \textbf{X}\boldsymbol{\beta}) ) \exp(-\frac{1}{2}\boldsymbol{\beta}^\top \boldsymbol{\Sigma}_b^{-1}\boldsymbol{\beta} )$
$ p(\boldsymbol{\beta}|\textbf{y}) \propto \exp(-\frac{1}{2} (\boldsymbol{\beta} - \boldsymbol{\mu}_\beta )^\top (\sigma^{-2} \textbf{X}^\top\textbf{X} + \boldsymbol{\Sigma_b}^{-1}) (\boldsymbol{\beta} - \boldsymbol{\mu}_\beta ))$
for $\boldsymbol{\mu}_\beta = \sigma^{-2} (\sigma^{-2} \textbf{X}^\top\textbf{X} + \boldsymbol{\Sigma_b}^{-1})^{-1} \textbf{X}^\top\textbf{y}$
and so $ p(\boldsymbol{\beta}|\textbf{y}) = N(\sigma^{-2} \textbf{W}^{-1} \textbf{X}^\top\textbf{y}, \textbf{W}^{-1} )$
for $\textbf{W} = \sigma^{-2} \textbf{X}^\top\textbf{X} + \boldsymbol{\Sigma_b}^{-1}$
No $\hat{\boldsymbol{\beta}}$ estimators, just a probability distribution of $\boldsymbol{\beta}$.
Predictions of $\textbf{y}$¶
Predictions of $\textbf{y}$ are made computing a predictive distribution for it as
$p(\textbf{y}_{\ast}|\textbf{y}) = \int p(\textbf{y}|\boldsymbol{\beta}) p(\boldsymbol{\beta}|\textbf{y}) d \boldsymbol{\beta}$¶
... and after some algebra it leads to
$p(\textbf{y}_{\ast}|\textbf{y}) = N(\sigma^{-2} \textbf{X}_\ast\textbf{W}^{-1} \textbf{X}^\top \textbf{y}, \textbf{X}_\ast \textbf{W}^{-1} \textbf{X}_\ast^\top )$¶
Algebra is more complicated here, and we don't always get analytical expressions of the posterior probability. E.g., assume a different prior distribution of the parameters. However there are ways to go around, we will go back to this point later.
More coding¶
Now let's make some plots to get the feeling of what we are saying... because "good teachers tackle examples" (or so I read https://twitter.com/JohnDCook/status/790645796630769664).
import pandas as pd
gbp = pd.read_csv('gbp_price.csv')
gbp.set_index(np.arange(gbp.shape[0])[::-1], inplace = True)
gbp = gbp.sort_values('week')
gbp
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(gbp.week, gbp.price, label = 'GBP price')
ax.set_xlabel('week')
ax.set_ylabel('price in USD')
ax.legend()
num_obsv = gbp.shape[0]
Y = gbp.price[:,None]
X = gbp.week[:,None]#np.random.uniform(0, 20, num_obsv)
# Make a design matrix with intercept
X = np.hstack([np.ones_like(Y), X])
print 'Design matrix:\n', X
Here is a good package in Python for machine learning: http://scikit-learn.org/
# Fit a linear regression using scikit learn (great package for machine learning!)
freq_reg = linear_model.LinearRegression(fit_intercept=False)
freq_reg.fit(X, Y)
print 'Coefficients:\n', freq_reg.coef_
#freq_reg.predict
#freq_reg.predict(X_star)
But we can also make it by ourselves
# Now make a prediction
num_star = 100
X_star = np.column_stack([np.ones((num_star,1)),np.linspace(0,25,num_star)])
# Frequentist
XTX_inv = np.linalg.inv(np.dot(X.T,X))
Y_starf = np.dot(np.dot(X_star,XTX_inv), np.dot(X.T,Y))
# Bayesian
sigma = .01 #Assume sigma is known for the moment, we will discuss later how to fit it
Cov_b = np.eye(2)
W = 1/sigma**2 * np.dot(X.T, X) + Cov_b
W_inv = np.linalg.inv(W)
Y_starb = 1/sigma**2 * np.dot(np.dot(X_star, W_inv), np.dot(X.T, Y))
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(X[:,1],Y)
ax.plot(X_star[:,1],Y_starf, 'b-')
ax.plot(X_star[:,1],Y_starb, 'r--')
ax.set_ylim(1.1,1.6)
So far, so good.
Are we happy with this model: $y = \beta_0 + \beta_1 x + \epsilon $?
What about $y = \beta_0 + \beta_1 x + \beta_2 x^2 + \epsilon $?
or $y = \sum_{i=0}^n \beta_i x^i + \epsilon $ for some larger $n$?
# Some functions that will help the demonstration
def phi_X(X_obsv, order):
# This function returns a design matrix for polynomial interplation
return np.column_stack([X_obsv**i for i in range(order+1)])
def poly_reg(poly_order):
# This function fits a polynomial regression and displays a plot
X_poly = phi_X(X[:,1], order = poly_order)
freq_reg.fit(X_poly, Y)
X_new = phi_X(X_obsv=np.linspace(0, 27, 100)[:,None], order=poly_order)
Y_new = freq_reg.predict(X_new)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(X_new[:,1],Y_new)
ax.scatter(X[:,1],Y)
ax.set_ylim(0,10)
ax.set_ylim(1.1,1.6)
poly_reg(2)
poly_reg(4)
poly_reg(10)
Image(filename='figs/Poincare.jpg', width = 200)
Image(filename='figs/scimethod_text.png', width = 1000)
Our data inputs, called $x$ by convention, usually have some clear meaning to us: e.g., date, spatial location, income, etc.
Here Poincare is talking about working on a space of features with a less immediate meaning, but with high explanatory value. Have you heard of
- Principal component analysis,
- Neural networks,
- Deep learning?
Let's go back to the regression example. We had that
$\textbf{y} = \textbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}$.
$\boldsymbol{\beta} = (\beta_1, \beta_2, ...)^\top$ and $\textbf{X} = ( \textbf{x}^0, \textbf{x}^1, ...) $.
So $\textbf{X}$, no matter how many columns it has, it is a function of a meaningful set of features $\textbf{x}$. This is $\textbf{X} = \phi(\textbf{x})$.
Our previous posterior distribution:
$p(\textbf{y}_{\ast}|\textbf{y}) = N(\sigma^{-2} \textbf{X}_\ast\textbf{W}^{-1} \textbf{X}^\top \textbf{y}, \textbf{X}_\ast \textbf{W}^{-1} \textbf{X}_\ast^\top )$¶
can then be written as:
$p(\textbf{y}_{\ast}|\textbf{y}) = N(\sigma^{-2} \phi(\textbf{x}_\ast)\textbf{W}^{-1} \phi(\textbf{x})^\top \textbf{y}, \phi(\textbf{x}_\ast) \textbf{W}^{-1} \phi(\textbf{x}_\ast)^\top )$¶
with
$\textbf{W} = \sigma^{-2} \phi(\textbf{x})^\top\phi(\textbf{x}) + \boldsymbol{\Sigma_b}^{-1}$¶
... after some more algebra, we can re-express the posterior moments as:
$\langle \textbf{y}_{\ast}|\textbf{y} \rangle = \phi(\textbf{x}_\ast) \boldsymbol{\Sigma_b}\phi(\textbf{x})^\top [\phi(\textbf{x})\boldsymbol{\Sigma_b}\phi(\textbf{x})^\top + \textbf{I}\sigma^2]^{-1} \textbf{y} $¶
$var( \textbf{y}_{\ast}|\textbf{y} ) = \phi(\textbf{x}_\ast) \boldsymbol{\Sigma_b}\phi(\textbf{x}_{\ast})^\top - \phi(\textbf{x}_\ast)\boldsymbol{\Sigma_b}\phi(\textbf{x})^\top [\phi(\textbf{x}) \boldsymbol{\Sigma_b}\phi(\textbf{x})^\top + \textbf{I}\sigma^2 ]^{-1} \phi(\textbf{x}) \boldsymbol{\Sigma_b}\phi(\textbf{x}_\ast)^\top$¶
Simple? Not really if you focus only in the algebra. The important thing to take from this is that the moments of the posterior distribution are defined in terms of inner products of $\phi(\textbf{x})$ or $\phi(\textbf{x}_{\ast})$.
$\phi(\textbf{x})\boldsymbol{\Sigma_b}\phi(\textbf{x})^\top $ is an inner product covariance function, and it is everything we need to define the posterior distribution.
We can define other covariance functions and use them instead. This allows us to use more flexible models.
Gaussian processes¶
First of all, forget about the assumption of i.i.d. We want observations $y_i$ and $y_j$ to be correlated, so that we can learn something from them and have some information to later predict a new $y_{\ast}$.
How do we model this? Well, assume that there is a function $(f)$ that whose realizations follow a multivariate Gaussian distribution with mean $\textbf{M}$ and covariance function $\textbf{K}$. Every finite collection $\{ f_{x_1}, f_{x_2}, \dots \}$ will follow this rule. Then this function is a Gaussian process and we denote it as $f \sim \mathcal{GP}(\textbf{M},\textbf{K})$.
Then we can assume that:
$y = f_x + \epsilon$
Noise is i.i.d., but realizations of $y$ are associated through the unobserved process (aka latent) $f$.
Everything is worked in the same way as before, in terms of our covariance function:
$\langle \textbf{y}_{\ast}|\textbf{y} \rangle = \textbf{K}_{\ast x} [\textbf{K}_{xx} + \textbf{I}\sigma^2]^{-1} \textbf{y} $¶
$var( \textbf{y}_{\ast}|\textbf{y} ) = \textbf{K}_{\ast \ast} - \textbf{K}_{\ast x} [\textbf{K}_{xx} + \textbf{I}\sigma^2]^{-1} \textbf{K}_{x \ast}$¶
Assume a covariance function in which the level of association between observations decreases exponentially:
$K(x_1, x_2) = exp(-.5 (x_1 - x_2)^2)$¶
def K_sqexp(X_1, X_2):
return np.exp(-.5*(X_1 - X_2.T)**2/1.)
x_1 = np.array([5]).reshape(1,1)
x_others = np.linspace(0,10)[:,None]
K_xx = K_sqexp(x_1, x_others)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x_others,K_xx.T, 'b-')
ax.vlines(x=x_1, ymin=0, ymax=1.1, color = 'r')
ax.set_ylim(0,1.1)
Now let's use this covariance in the GBP data.
def predictive_mean_Y(X,Y,K,X_new):
sigma = .001 #We are still assuming this
K_xx = K(X,X)
K_newx = K(X_new,X)
W_inv = np.linalg.inv(K_xx + np.eye(Y.size)*sigma)
return np.dot(K_newx,np.dot(W_inv,Y))
def predictive_var_Y(X,Y,K,X_new):
sigma = .001 #We are still assuming this
K_xx = K(X,X)
K_newx = K(X_new,X)
K_newnew = K(X_new,X_new)
W_inv = np.linalg.inv(K_xx + np.eye(Y.size)*sigma)
return K_newnew - np.dot(K_newx,np.dot(W_inv,K_newx.T))
Y = gbp.price[:,None]
X = gbp.week[:,None]
X_new=np.linspace(1, 25, 100)[:,None]
Y_pm = predictive_mean_Y(X,Y,K_sqexp, X_new)
Y_pv = np.diag(predictive_var_Y(X,Y,K_sqexp, X_new))[:,None]
cilb = Y_pm - 1.96*np.sqrt(Y_pv)
ciub = Y_pm + 1.96*np.sqrt(Y_pv)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(X,Y)
ax.plot(X_new, Y_pm)
ax.fill_between(X_new.flatten(), cilb.flatten(), ciub.flatten(), facecolor = 'r', alpha = .2, label = '95% credible interval')
ax.set_ylim(1.1,1.6)
The fit to the data seems good, but the credible intervals not so much. It still can be improved.
We can use a different covariace functions, of use the same one, but add new parameters to it. For instance, we can start by scaling the input differences:
$K(x_1, x_2) = exp(-\frac{(x_1 - x_2)^2}{2\ell^2})$¶
We will go through other examples, but we will be using a specialized package for it https://github.com/SheffieldML/GPy
# These are different covariance functions
cov_cst = GPy.kern.Bias(1)
cov_lmd = GPy.kern.Linear(1)
cov_rbf = GPy.kern.RBF(1)
cov_per = GPy.kern.PeriodicMatern32(1,period=20)
m = GPy.models.GPRegression(X=X, Y=Y, kernel = cov_per, normalizer=True)
m.optimize()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(X,Y)
ax.set_ylim(1.1,1.6)
m.plot(ax = ax, plot_data=False)
The non-linearity of the predictions is achieved not by defining a set of parameters $\beta_i$ as before, but by fitting a parameter corresponding to each location $x_i$.
Parametric and non-parametric models¶
A parametric model like $y = \beta_0 + \beta_1 x$, is a model that summarizes all the observed data in a finite set of parameters. Making predictions for points outside the dataset used for fitting the model is trivial, and is just based on the parameters.
A non-parametric model doesn't mean that it is a model with no parameters, it's actually the opposite. The number of parameters increase with the number of observations, and prediction of new points require the estimation of new parameters. Therefore, the number of parameters needed to describe certain phenomenon is infinite!
While the computation of so many parameters can be a buden, it is also the large number of parameters that allows the model to represent the characteristics of the data much better and in a less rigid way than parametric models.
Geostatistics¶
[Recommended reading: http://www2.stat.duke.edu/~fei/samsi/Readings/DiggTawnMoye1988.pdf]
Geostatistics is about doing statistics in space (not outer space, more like 2-D space). Traditionally it dealt with the case of predicting a linear model of a Gaussian process in space. Something like the model we described before:
$y = f_x + \epsilon$, where $f \sim \mathcal{GP}(\textbf{M},\textbf{K})$.
The methodology for doing this spatial interpolation is generally regarded as kriging. You have have seen it before (last week). Here we will fit this model with the Bayesian approach we have been discussing.
Now a toy example.
With Gaussian processes we assume there is a latent function $f$ that we do not observe but that represents the actual process we are trying to describe. Suppose we were able to see this function for once.
def latent(X):
# Some crazy function
return 1+.2*np.sqrt(X[:,0])**2 + np.sin(X[:,0]) + np.cos(X[:,1]/3)
space = np.meshgrid(np.linspace(0,10,80), np.linspace(0,10,80))
X_space = np.vstack([space[0].flatten(), space[1].flatten()]).T
f_space = latent(X_space)
H = f_space.reshape(space[0].shape, order = 'C')
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
fig = plt.figure(figsize=(8,8))
ax0 = fig.add_subplot(111,projection='3d')
surf = ax0.plot_surface(space[0], space[1], H, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False)
ax0.set_zlim(-.7, 3.5)
ax0.zaxis.set_major_locator(LinearLocator(10))
ax0.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
ax0.set_zlim(0,5)
ax0.set_title('Ground truth')
fig.colorbar(surf, shrink=0.5, aspect=5)
Now, below is something closer to what we are usually able to observe: just a 'small' set of points distorted by noise.
num_obsv = 100
X_obsv = np.hstack([np.random.uniform(0,10,num_obsv)[:,None],np.random.uniform(0,10,num_obsv)[:,None]])
Y_obsv = latent(X_obsv)[:,None] + np.random.normal(0,.5,num_obsv)[:,None]
fig = plt.figure(figsize=(12,8))
ax0 = fig.add_subplot(121,projection='3d')
surf = ax0.plot_surface(space[0], space[1], H, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False, vmin=0, vmax=5)
ax0.set_zlim(0, 5)
ax0.zaxis.set_major_locator(LinearLocator(10))
ax0.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
ax0.set_title('Ground truth')
ax1 = fig.add_subplot(122, projection='3d')
ax1.scatter(X_obsv[:,0], X_obsv[:,1], zs = Y_obsv, c = Y_obsv, lw=0, s = 40, vmin=-1, vmax=4, cmap=cm.coolwarm)
ax1.scatter(X_obsv[:,0], X_obsv[:,1], zs = 0, c = 'gray', marker='x')
ax1.set_zlim(0, 5)
ax1.set_title('Observed data')
And after doing a spatial regression, or a GP regression in a 2-dimensional input space, we get the following:
m = GPy.models.GPRegression(X_obsv, Y_obsv)#, kernel=GPy.kern.RBF(2) + GPy.kern.Bias(2) )
m.optimize()
Y_pm, Y_pv = m.predict(X_space)
Hp = Y_pm.reshape(space[0].shape, order = 'C')
fig = plt.figure(figsize=(12,8))
ax0 = fig.add_subplot(121,projection='3d')
surf = ax0.plot_surface(space[0], space[1], H, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False, vmin=0, vmax=5)
ax0.set_zlim(0, 5)
ax0.zaxis.set_major_locator(LinearLocator(10))
ax0.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
ax0.set_title('Ground truth')
ax1 = fig.add_subplot(122,projection='3d')
ax1.scatter(X_obsv[:,0], X_obsv[:,1], zs = 0, c = 'gray', marker='x')
surf = ax1.plot_surface(space[0], space[1], Hp, rstride=1, cstride=1, cmap=cm.coolwarm,linewidth=0, antialiased=False, vmin=0, vmax=5)
ax1.set_zlim(0, 5)
ax1.zaxis.set_major_locator(LinearLocator(10))
ax1.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
ax1.set_title('Geostatistic estimate')
The probabilistic estimate draws a smooth surface across the 10 x 10 area based on the relative distance from each location $(x_1,x_2)$ to the set of observations.
Traditionally geostatistics was seen as a tool for modeling continuous data. However, that is no longer the case. It is common to see applications where discrete processes are modeled assuming an underlying continuous process. Similar to Generalized Linear Models (https://en.wikipedia.org/wiki/Generalized_linear_model).
Think of a Poisson process $y \sim Poisson(\lambda)$. We have that the number of cases within an area is defined by an intensity function $\lambda$. While $y$ is discrete, $\lambda$ is continuous and can be modeled using a Gaussian process.
The only requirement is that $\lambda > 0$. This can be satisfied by defining: $\lambda = exp(f)$, where $f \sim \mathcal{GP}$.
In general we can find a monotonic transformation on $f$ that is linked to the process we want to model.
Example: clasification problem¶
def classify(X, noise = True):
f = latent(X)
f -= 2
if noise:
f += np.random.normal(0,1,f.size).reshape(f.shape)
c = np.zeros_like(f)
c[f>0] = 1
return c
space = np.meshgrid(np.linspace(0,10,80), np.linspace(0,10,80))
X_space = np.vstack([space[0].flatten(), space[1].flatten()]).T
c_space = classify(X_space, noise=False)
Hc = c_space.reshape(space[0].shape, order = 'C')
num_obsv = 200
X_obsv = np.hstack([np.random.uniform(0,10,num_obsv)[:,None],np.random.uniform(0,10,num_obsv)[:,None]])
C_obsv = classify(X_obsv)[:,None]
C0 = C_obsv.flatten() == 0
C1 = C_obsv.flatten() == 1
fig = plt.figure(figsize=(12,10))
#ax0 = fig.add_subplot(121)
#ax0.set_title('Ground truth')
#img0 = ax0.imshow(Hc, origin='lower', extent = (0,10,0,10), vmin=-.2, vmax=1.2)
ax1 = fig.add_subplot(111)
ax1.set_title('Observed data')
#img = ax1.scatter(X_obsv[:,0], X_obsv[:,1], c = C_obsv, lw=0, s = 80, vmin=-.2, vmax=1.2)
ax1.scatter(X_obsv[C0,0], X_obsv[C0,1], marker = 'o', color = 'blue', s = 90)
ax1.scatter(X_obsv[C1,0], X_obsv[C1,1], marker = 'v', color = 'red', s = 90)
ax1.set_ylim(0,10)
ax1.set_xlim(0,10)
m = GPy.models.GPClassification(X_obsv, C_obsv, kernel=GPy.kern.Matern32(2))
m.optimize()
Y_pm, Y_pv = m.predict(X_space)
Hp = Y_pm.reshape(space[0].shape, order = 'C')
fig = plt.figure(figsize=(8,8))
#ax0 = fig.add_subplot(121)
#ax0.set_title('Ground truth')
#img0 = ax0.imshow(Hc, origin='lower', extent = (0,10,0,10), vmin=-.2, vmax=1.2)
ax1 = fig.add_subplot(111)
ax1.set_title('Probabilistic estimate')
ax1.scatter(X_obsv[C0,0], X_obsv[C0,1], marker = 'o', color = 'blue', s = 90)
ax1.scatter(X_obsv[C1,0], X_obsv[C1,1], marker = 'v', color = 'darkred', s = 90)
Hp = Y_pm.reshape(space[0].shape, order = 'C')
img = ax1.imshow(Hp, origin='lower', extent = (0,10,0,10), vmin=0, vmax=1)
fig.colorbar(img, fraction = 0.046, pad=0.04)
Bayesian computing¶
Now let's think on the 'how to' of Bayesian computing. We are interested in the posterior distribution $p(\theta|y)$:
$p(\theta|y) = \frac{p(y|\theta) p(\theta)}{p(y)}$¶
We can evaluate the numerator pointwise, but the difficult part here (sometimes) is to find
$ p(y) = \int p(y|\theta) p(\theta) d \theta $¶
This quantity normalizes the posterior distribution, i.e. makes it intergrate to 1.
In many cases there is no analycal solution to it. Which means that we cannot come up with simple algebraic formulation to compute it. In these cases we have to be more creative to solve it. The most common approach is to use Markov Chain Monte Carlo (MCMC) sampling. There is no time to see this topic today, but here you can find a detailed tutorial on it: http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/
Rejection sampling¶
Notice that $ \int p(y|\theta) p(\theta) d \theta $ is the expected value of $p(y|\theta)$ wrt. the distribution $p(\theta)$. An easy way of estimating an expected value is just to draw enough samples of the variable of interest and average them:
$\int p(y|\theta) p(\theta) d \theta \approx \frac{1}{N}\sum_{i=1}^N p(y|\theta_i)$¶
Then we have a clear plan. We will try to generate 'samples' form: $p(y|\theta) \times p(\theta)$. And after getting enough samples we can approxiamte the distribution. Easy, right? The only problem now is to sample in a clever way. The space of $\theta$ that we need to explore can be very large!
Suppose that we have this posterior distribution.
def post_dist(y, mu1=7, mu2=3.8):
return .3*np.exp(-.5 * (y-mu1)**2.)/np.sqrt(2.*np.pi) + .7*np.exp(-.5 * (y-mu2)**2.)/np.sqrt(2.*np.pi)
_y = np.linspace(-5,20,200)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(_y, post_dist(_y))
ax.set_ylim(-.05,.35)
ax.plot(-5*np.random.random(10), np.zeros(10), 'rx', ms=8)
ax.plot(11+5*np.random.random(10), np.zeros(10), 'gx', ms=8)
ax.plot(11*np.random.random(10), np.zeros(10), 'kx', ms=8)
If we just had a proxy distribution $q(x)$ that 'covers' our distribution of interest, we could sample from that one.
With cover I mean: $p(x) \leq k q(x)$, for some $k$.
def proxy_dist(y, mu, sigma2, k):
return k * np.exp(-.5 * (y-mu)**2./sigma2)/np.sqrt(2.*np.pi*sigma2)
_y = np.linspace(-5,20,200)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(_y, post_dist(_y))
ax.plot(_y, proxy_dist(_y, mu=4.6, sigma2=4, k=1.6), 'r--', lw=2)
ax.set_ylim(-.05,.35)
ax.plot(stats.norm.rvs(4.6,2,20), np.zeros(20), 'rx', ms=8)
Because we can compute our target distribution pointwise, and we know $k$ and $q(x)$, we can go further and discard the samples that are not appropriate.
Y_sample = stats.norm.rvs(4.6,2,100)
q_sample = proxy_dist(Y_sample, mu=4.6, sigma2=4, k=1.6)
p_sample = post_dist(Y_sample)
U_sample = np.array([stats.uniform.rvs(0,kq_i,1) for kq_i in q_sample]).flatten()
Y_thinned = Y_sample[U_sample<=p_sample]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(_y, post_dist(_y))
ax.plot(_y, proxy_dist(_y, mu=4.6, sigma2=4, k=1.6), 'r--', lw=2)
ax.plot(Y_sample, np.zeros_like(Y_sample), 'rx', ms=8)
ax.plot(Y_thinned, np.zeros_like(Y_thinned), 'bx', ms=8)
ax.hist(Y_thinned, normed=True, color = 'gray', alpha=.1, bins=20)
ax.set_ylim(-.05,.35)
Now the problem is to find $q(x)$ and $k$. While the tequnique is simple and intuitive, it is not very useful for highdimensional problems.
If you are interested in the topic, you might want to have a look at these methods:
-Importance Sampling
-Metropolis-Hasting Algorithm
-Gibbs sampling
Look here https://www.cs.ubc.ca/~murphyk/MLbook/
Other methods for Bayesian computing¶
Maximum likelihood type II¶
Consists of maximizing the marginal likelihood $p(y)$ wrt. the hyperparameters. In the case of GP, these are the parameters of the Covariance function.
It is faster than MCMC, but there are drawbacks: optimization is not easy, and overfitting my occur.
Laplace approximation¶
Approximates the posterior distribution with a Gaussian approximation.
(R-)INLA¶
Uses a Gaussian Markov Field Approximation to the GP in combination with a modified Laplace approximation.