Example of simple GP fit, adapted from Stan's example-models repository.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
In [2]:
x = np.array([-5, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4, 
-3.9, -3.8, -3.7, -3.6, -3.5, -3.4, -3.3, -3.2, -3.1, -3, -2.9, 
-2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2, -1.9, -1.8, 
-1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1, -0.9, -0.8, -0.7, 
-0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 
0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 
1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 
3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 
4.5, 4.6, 4.7, 4.8, 4.9, 5])

y = np.array([1.04442478194401, 0.948306088493654, 0.357037759697332, 0.492336514646604, 
0.520651364364746, 0.112629866592809, 0.470995468454158, -0.168442254267804, 
0.0720344402575861, -0.188108980535916, -0.0160163306512027, 
-0.0388792158617705, -0.0600673630622568, 0.113568725264636, 
0.447160403837629, 0.664421188556779, -0.139510743820276, 0.458823971660986, 
0.141214654640904, -0.286957663528091, -0.466537724021695, -0.308185884317105, 
-1.57664872694079, -1.44463024170082, -1.51206214603847, -1.49393593601901, 
-2.02292464164487, -1.57047488853653, -1.22973445533419, -1.51502367058357, 
-1.41493587255224, -1.10140254663611, -0.591866485375275, -1.08781838696462, 
-0.800375653733931, -1.00764767602679, -0.0471028950122742, -0.536820626879737, 
-0.151688056391446, -0.176771681318393, -0.240094952335518, -1.16827876746502, 
-0.493597351974992, -0.831683011472805, -0.152347043914137, 0.0190364158178343, 
-1.09355955218051, -0.328157917911376, -0.585575679802941, -0.472837120425201, 
-0.503633622750049, -0.0124446353828312, -0.465529814250314, 
-0.101621725887347, -0.26988462590405, 0.398726664193302, 0.113805181040188, 
0.331353802465398, 0.383592361618461, 0.431647298655434, 0.580036473774238, 
0.830404669466897, 1.17919105883462, 0.871037583886711, 1.12290553424174, 
0.752564860804382, 0.76897960270623, 1.14738839410786, 0.773151715269892, 
0.700611498974798, 0.0412951045437818, 0.303526087747629, -0.139399513324585, 
-0.862987735433697, -1.23399179134008, -1.58924289116396, -1.35105117911049, 
-0.990144529089174, -1.91175364127672, -1.31836236129543, -1.65955735224704, 
-1.83516148300526, -2.03817062501248, -1.66764011409214, -0.552154350554687, 
-0.547807883952654, -0.905389222477036, -0.737156477425302, -0.40211249920415, 
0.129669958952991, 0.271142753510592, 0.176311762529962, 0.283580281859344, 
0.635808289696458, 1.69976647982837, 1.10748978734239, 0.365412229181044, 
0.788821368082444, 0.879731888124867, 1.02180766619069, 0.551526067300283])

N = len(y)

Original Stan model:

// Fit a Gaussian process's hyperparameters
// for squared exponential prior

data {
  int<lower=1> N;
  vector[N] x;
  vector[N] y;
}
transformed data {
  vector[N] mu;
  for (i in 1:N) 
    mu[i] <- 0;
}
parameters {
  real<lower=0> eta_sq;
  real<lower=0> rho_sq;
  real<lower=0> sigma_sq;
}
model {
  matrix[N,N] Sigma;

  // off-diagonal elements
  for (i in 1:(N-1)) {
    for (j in (i+1):N) {
      Sigma[i,j] <- eta_sq * exp(-rho_sq * pow(x[i] - x[j],2));
      Sigma[j,i] <- Sigma[i,j];
    }
  }

  // diagonal elements
  for (k in 1:N)
    Sigma[k,k] <- eta_sq + sigma_sq; // + jitter

  eta_sq ~ cauchy(0,5);
  rho_sq ~ cauchy(0,5);
  sigma_sq ~ cauchy(0,5);

  y ~ multi_normal(mu,Sigma);
}
In [15]:
from pymc3 import Model, MvNormal, HalfCauchy, sample, traceplot
import theano.tensor as T
In [14]:
with Model() as gp_fit:
    
    mu = np.zeros(N)
    
    η_sq = HalfCauchy('η_sq', 5)
    ρ_sq = HalfCauchy('ρ_sq', 5)
    σ_sq = HalfCauchy('σ_sq', 5)
    
    Sigma = T.stack([[η_sq * np.exp(-ρ_sq * (x[i] - x[j])**2) * (σ_sq*int(i==j) or 1) 
                              for i in range(N)] 
                                 for j in range(N)])
    
    obs = MvNormal('obs', mu, Sigma, observed=y)
Applied log-transform to η_sq and added transformed η_sq_log to model.
Applied log-transform to ρ_sq and added transformed ρ_sq_log to model.
Applied log-transform to σ_sq and added transformed σ_sq_log to model.
In [16]:
with gp_fit:
    
    gp_trace = sample(1000)
Assigned NUTS to η_sq_log
Assigned NUTS to ρ_sq_log
Assigned NUTS to σ_sq_log
---------------------------------------------------------------------------
RecursionError                            Traceback (most recent call last)
<ipython-input-16-87dffea1edf3> in <module>()
      1 with gp_fit:
      2 
----> 3     gp_trace = sample(1000)

/Users/fonnescj/Repositories/pymc3/pymc3/sampling.py in sample(draws, step, start, trace, chain, njobs, tune, progressbar, model, random_seed)
    122     model = modelcontext(model)
    123 
--> 124     step = assign_step_methods(model, step)
    125 
    126     if njobs is None:

/Users/fonnescj/Repositories/pymc3/pymc3/sampling.py in assign_step_methods(model, step, methods)
     67 
     68     # Instantiate all selected step methods
---> 69     steps += [s(vars=selected_steps[s]) for s in selected_steps if selected_steps[s]]
     70 
     71     if len(steps)==1:

/Users/fonnescj/Repositories/pymc3/pymc3/sampling.py in <listcomp>(.0)
     67 
     68     # Instantiate all selected step methods
---> 69     steps += [s(vars=selected_steps[s]) for s in selected_steps if selected_steps[s]]
     70 
     71     if len(steps)==1:

/Users/fonnescj/Repositories/pymc3/pymc3/step_methods/nuts.py in __init__(self, vars, scaling, step_scale, is_cov, state, Emax, target_accept, gamma, k, t0, model, profile, **kwargs)
     67 
     68         if isinstance(scaling, dict):
---> 69             scaling = guess_scaling(Point(scaling, model=model), model=model, vars = vars)
     70 
     71 

/Users/fonnescj/Repositories/pymc3/pymc3/tuning/scaling.py in guess_scaling(point, vars, model)
    107     model = modelcontext(model)
    108     try:
--> 109         h = find_hessian_diag(point, vars, model=model)
    110     except NotImplementedError:
    111         h = fixed_hessian(point, vars, model=model)

/Users/fonnescj/Repositories/pymc3/pymc3/tuning/scaling.py in find_hessian_diag(point, vars, model)
    101     """
    102     model = modelcontext(model)
--> 103     H = model.fastfn(hessian_diag(model.logpt, vars))
    104     return H(Point(point, model=model))
    105 

/Users/fonnescj/Repositories/pymc3/pymc3/memoize.py in memoizer(*args, **kwargs)
     12 
     13         if key not in cache:
---> 14             cache[key] = obj(*args, **kwargs)
     15 
     16         return cache[key]

/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py in hessian_diag(f, vars)
    101 
    102     if vars:
--> 103         return -t.concatenate([hessian_diag1(f, v) for v in vars], axis=0)
    104     else:
    105         return empty_gradient

/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py in <listcomp>(.0)
    101 
    102     if vars:
--> 103         return -t.concatenate([hessian_diag1(f, v) for v in vars], axis=0)
    104     else:
    105         return empty_gradient

/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py in hessian_diag1(f, v)
     92         return gradient1(g[i], v)[i]
     93 
---> 94     return theano.map(hess_ii, idx)[0]
     95 
     96 

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/scan_module/scan_views.py in map(fn, sequences, non_sequences, truncate_gradient, go_backwards, mode, name)
     67                      go_backwards=go_backwards,
     68                      mode=mode,
---> 69                      name=name)
     70 
     71 

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/scan_module/scan.py in scan(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict)
    743     # and outputs that needs to be separated
    744 
--> 745     condition, outputs, updates = scan_utils.get_updates_and_outputs(fn(*args))
    746     if condition is not None:
    747         as_while = True

/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py in hess_ii(i)
     90 
     91     def hess_ii(i):
---> 92         return gradient1(g[i], v)[i]
     93 
     94     return theano.map(hess_ii, idx)[0]

/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py in gradient1(f, v)
     42 def gradient1(f, v):
     43     """flat gradient of f wrt v"""
---> 44     return t.flatten(t.grad(f, v, disconnected_inputs='warn'))
     45 
     46 

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gradient.py in grad(cost, wrt, consider_constant, disconnected_inputs, add_names, known_grads, return_disconnected, null_gradients)
    458 
    459     var_to_app_to_idx = _populate_var_to_app_to_idx(
--> 460         outputs, wrt, consider_constant)
    461 
    462     # build a dict mapping var to the gradient of cost with respect to var

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gradient.py in _populate_var_to_app_to_idx(outputs, wrt, consider_constant)
    885     # add all variables that are true ancestors of the cost
    886     for output in outputs:
--> 887         account_for(output)
    888 
    889     # determine which variables have elements of wrt as a true

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gradient.py in account_for(var)
    881                 if i not in idx:
    882                     idx.append(i)
--> 883                 account_for(ipt)
    884 
    885     # add all variables that are true ancestors of the cost

... last 1 frames repeated, from the frame below ...

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gradient.py in account_for(var)
    881                 if i not in idx:
    882                     idx.append(i)
--> 883                 account_for(ipt)
    884 
    885     # add all variables that are true ancestors of the cost

RecursionError: maximum recursion depth exceeded
In [ ]:
traceplot(gp_trace)