Skip to content

Instantly share code, notes, and snippets.

@fonnesbeck

fonnesbeck/Simple GP.ipynb

Created February 18, 2016 02:15
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save fonnesbeck/f2da55c0e2f6d8f35a25 to your computer and use it in GitHub Desktop.
Save fonnesbeck/f2da55c0e2f6d8f35a25 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
Simple GP.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Example of simple GP fit, adapted from Stan's [example-models repository](https://github.com/stan-dev/example-models/blob/master/misc/gaussian-process/gp-fit.stan)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"sns.set()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"x = np.array([-5, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4, \n",
"-3.9, -3.8, -3.7, -3.6, -3.5, -3.4, -3.3, -3.2, -3.1, -3, -2.9, \n",
"-2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1, -2, -1.9, -1.8, \n",
"-1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1, -0.9, -0.8, -0.7, \n",
"-0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, \n",
"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, \n",
"1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, \n",
"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, \n",
"4.5, 4.6, 4.7, 4.8, 4.9, 5])\n",
"\n",
"y = np.array([1.04442478194401, 0.948306088493654, 0.357037759697332, 0.492336514646604, \n",
"0.520651364364746, 0.112629866592809, 0.470995468454158, -0.168442254267804, \n",
"0.0720344402575861, -0.188108980535916, -0.0160163306512027, \n",
"-0.0388792158617705, -0.0600673630622568, 0.113568725264636, \n",
"0.447160403837629, 0.664421188556779, -0.139510743820276, 0.458823971660986, \n",
"0.141214654640904, -0.286957663528091, -0.466537724021695, -0.308185884317105, \n",
"-1.57664872694079, -1.44463024170082, -1.51206214603847, -1.49393593601901, \n",
"-2.02292464164487, -1.57047488853653, -1.22973445533419, -1.51502367058357, \n",
"-1.41493587255224, -1.10140254663611, -0.591866485375275, -1.08781838696462, \n",
"-0.800375653733931, -1.00764767602679, -0.0471028950122742, -0.536820626879737, \n",
"-0.151688056391446, -0.176771681318393, -0.240094952335518, -1.16827876746502, \n",
"-0.493597351974992, -0.831683011472805, -0.152347043914137, 0.0190364158178343, \n",
"-1.09355955218051, -0.328157917911376, -0.585575679802941, -0.472837120425201, \n",
"-0.503633622750049, -0.0124446353828312, -0.465529814250314, \n",
"-0.101621725887347, -0.26988462590405, 0.398726664193302, 0.113805181040188, \n",
"0.331353802465398, 0.383592361618461, 0.431647298655434, 0.580036473774238, \n",
"0.830404669466897, 1.17919105883462, 0.871037583886711, 1.12290553424174, \n",
"0.752564860804382, 0.76897960270623, 1.14738839410786, 0.773151715269892, \n",
"0.700611498974798, 0.0412951045437818, 0.303526087747629, -0.139399513324585, \n",
"-0.862987735433697, -1.23399179134008, -1.58924289116396, -1.35105117911049, \n",
"-0.990144529089174, -1.91175364127672, -1.31836236129543, -1.65955735224704, \n",
"-1.83516148300526, -2.03817062501248, -1.66764011409214, -0.552154350554687, \n",
"-0.547807883952654, -0.905389222477036, -0.737156477425302, -0.40211249920415, \n",
"0.129669958952991, 0.271142753510592, 0.176311762529962, 0.283580281859344, \n",
"0.635808289696458, 1.69976647982837, 1.10748978734239, 0.365412229181044, \n",
"0.788821368082444, 0.879731888124867, 1.02180766619069, 0.551526067300283])\n",
"\n",
"N = len(y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Original Stan model:\n",
"\n",
"```\n",
"// Fit a Gaussian process's hyperparameters\n",
"// for squared exponential prior\n",
"\n",
"data {\n",
" int<lower=1> N;\n",
" vector[N] x;\n",
" vector[N] y;\n",
"}\n",
"transformed data {\n",
" vector[N] mu;\n",
" for (i in 1:N) \n",
" mu[i] <- 0;\n",
"}\n",
"parameters {\n",
" real<lower=0> eta_sq;\n",
" real<lower=0> rho_sq;\n",
" real<lower=0> sigma_sq;\n",
"}\n",
"model {\n",
" matrix[N,N] Sigma;\n",
"\n",
" // off-diagonal elements\n",
" for (i in 1:(N-1)) {\n",
" for (j in (i+1):N) {\n",
" Sigma[i,j] <- eta_sq * exp(-rho_sq * pow(x[i] - x[j],2));\n",
" Sigma[j,i] <- Sigma[i,j];\n",
" }\n",
" }\n",
"\n",
" // diagonal elements\n",
" for (k in 1:N)\n",
" Sigma[k,k] <- eta_sq + sigma_sq; // + jitter\n",
"\n",
" eta_sq ~ cauchy(0,5);\n",
" rho_sq ~ cauchy(0,5);\n",
" sigma_sq ~ cauchy(0,5);\n",
"\n",
" y ~ multi_normal(mu,Sigma);\n",
"}\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from pymc3 import Model, MvNormal, HalfCauchy, sample, traceplot\n",
"import theano.tensor as T"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Applied log-transform to η_sq and added transformed η_sq_log to model.\n",
"Applied log-transform to ρ_sq and added transformed ρ_sq_log to model.\n",
"Applied log-transform to σ_sq and added transformed σ_sq_log to model.\n"
]
}
],
"source": [
"with Model() as gp_fit:\n",
" \n",
" mu = np.zeros(N)\n",
" \n",
" η_sq = HalfCauchy('η_sq', 5)\n",
" ρ_sq = HalfCauchy('ρ_sq', 5)\n",
" σ_sq = HalfCauchy('σ_sq', 5)\n",
" \n",
" Sigma = T.stack([[η_sq * np.exp(-ρ_sq * (x[i] - x[j])**2) * (σ_sq*int(i==j) or 1) \n",
" for i in range(N)] \n",
" for j in range(N)])\n",
" \n",
" obs = MvNormal('obs', mu, Sigma, observed=y)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Assigned NUTS to η_sq_log\n",
"Assigned NUTS to ρ_sq_log\n",
"Assigned NUTS to σ_sq_log\n"
]
},
{
"ename": "RecursionError",
"evalue": "maximum recursion depth exceeded",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mRecursionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-16-87dffea1edf3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mgp_fit\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mgp_trace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/sampling.py\u001b[0m in \u001b[0;36msample\u001b[0;34m(draws, step, start, trace, chain, njobs, tune, progressbar, model, random_seed)\u001b[0m\n\u001b[1;32m 122\u001b[0m \u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodelcontext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 123\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 124\u001b[0;31m \u001b[0mstep\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0massign_step_methods\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstep\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 125\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 126\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mnjobs\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/sampling.py\u001b[0m in \u001b[0;36massign_step_methods\u001b[0;34m(model, step, methods)\u001b[0m\n\u001b[1;32m 67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0;31m# Instantiate all selected step methods\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 69\u001b[0;31m \u001b[0msteps\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mselected_steps\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ms\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mselected_steps\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mselected_steps\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 70\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 71\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msteps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/sampling.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0;31m# Instantiate all selected step methods\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 69\u001b[0;31m \u001b[0msteps\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvars\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mselected_steps\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ms\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mselected_steps\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mselected_steps\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 70\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 71\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msteps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/step_methods/nuts.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, vars, scaling, step_scale, is_cov, state, Emax, target_accept, gamma, k, t0, model, profile, **kwargs)\u001b[0m\n\u001b[1;32m 67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscaling\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdict\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 69\u001b[0;31m \u001b[0mscaling\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mguess_scaling\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mPoint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscaling\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvars\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvars\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 70\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 71\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/tuning/scaling.py\u001b[0m in \u001b[0;36mguess_scaling\u001b[0;34m(point, vars, model)\u001b[0m\n\u001b[1;32m 107\u001b[0m \u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodelcontext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 108\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 109\u001b[0;31m \u001b[0mh\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfind_hessian_diag\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpoint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvars\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 110\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 111\u001b[0m \u001b[0mh\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfixed_hessian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpoint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvars\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/tuning/scaling.py\u001b[0m in \u001b[0;36mfind_hessian_diag\u001b[0;34m(point, vars, model)\u001b[0m\n\u001b[1;32m 101\u001b[0m \"\"\"\n\u001b[1;32m 102\u001b[0m \u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodelcontext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 103\u001b[0;31m \u001b[0mH\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfastfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhessian_diag\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlogpt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvars\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 104\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mH\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mPoint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpoint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 105\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/memoize.py\u001b[0m in \u001b[0;36mmemoizer\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 13\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mkey\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 14\u001b[0;31m \u001b[0mcache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 16\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py\u001b[0m in \u001b[0;36mhessian_diag\u001b[0;34m(f, vars)\u001b[0m\n\u001b[1;32m 101\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 102\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mvars\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 103\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mhessian_diag1\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mvars\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 104\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 105\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mempty_gradient\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 101\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 102\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mvars\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 103\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mhessian_diag1\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mvars\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 104\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 105\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mempty_gradient\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py\u001b[0m in \u001b[0;36mhessian_diag1\u001b[0;34m(f, v)\u001b[0m\n\u001b[1;32m 92\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mgradient1\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 93\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 94\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mtheano\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhess_ii\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0midx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 95\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 96\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/scan_module/scan_views.py\u001b[0m in \u001b[0;36mmap\u001b[0;34m(fn, sequences, non_sequences, truncate_gradient, go_backwards, mode, name)\u001b[0m\n\u001b[1;32m 67\u001b[0m \u001b[0mgo_backwards\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgo_backwards\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 69\u001b[0;31m name=name)\n\u001b[0m\u001b[1;32m 70\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 71\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/scan_module/scan.py\u001b[0m in \u001b[0;36mscan\u001b[0;34m(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict)\u001b[0m\n\u001b[1;32m 743\u001b[0m \u001b[0;31m# and outputs that needs to be separated\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 744\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 745\u001b[0;31m \u001b[0mcondition\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mupdates\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscan_utils\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_updates_and_outputs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 746\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mcondition\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 747\u001b[0m \u001b[0mas_while\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py\u001b[0m in \u001b[0;36mhess_ii\u001b[0;34m(i)\u001b[0m\n\u001b[1;32m 90\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 91\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mhess_ii\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 92\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mgradient1\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 93\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 94\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mtheano\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhess_ii\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0midx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/Repositories/pymc3/pymc3/theanof.py\u001b[0m in \u001b[0;36mgradient1\u001b[0;34m(f, v)\u001b[0m\n\u001b[1;32m 42\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mgradient1\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 43\u001b[0m \u001b[0;34m\"\"\"flat gradient of f wrt v\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 44\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgrad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdisconnected_inputs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'warn'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 45\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 46\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gradient.py\u001b[0m in \u001b[0;36mgrad\u001b[0;34m(cost, wrt, consider_constant, disconnected_inputs, add_names, known_grads, return_disconnected, null_gradients)\u001b[0m\n\u001b[1;32m 458\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 459\u001b[0m var_to_app_to_idx = _populate_var_to_app_to_idx(\n\u001b[0;32m--> 460\u001b[0;31m outputs, wrt, consider_constant)\n\u001b[0m\u001b[1;32m 461\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 462\u001b[0m \u001b[0;31m# build a dict mapping var to the gradient of cost with respect to var\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gradient.py\u001b[0m in \u001b[0;36m_populate_var_to_app_to_idx\u001b[0;34m(outputs, wrt, consider_constant)\u001b[0m\n\u001b[1;32m 885\u001b[0m \u001b[0;31m# add all variables that are true ancestors of the cost\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 886\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0moutput\u001b[0m \u001b[0;32min\u001b[0m \u001b[0moutputs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 887\u001b[0;31m \u001b[0maccount_for\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 888\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 889\u001b[0m \u001b[0;31m# determine which variables have elements of wrt as a true\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gradient.py\u001b[0m in \u001b[0;36maccount_for\u001b[0;34m(var)\u001b[0m\n\u001b[1;32m 881\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0midx\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 882\u001b[0m \u001b[0midx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 883\u001b[0;31m \u001b[0maccount_for\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mipt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 884\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 885\u001b[0m \u001b[0;31m# add all variables that are true ancestors of the cost\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"... last 1 frames repeated, from the frame below ...\n",
"\u001b[0;32m/Users/fonnescj/anaconda3/lib/python3.5/site-packages/theano/gradient.py\u001b[0m in \u001b[0;36maccount_for\u001b[0;34m(var)\u001b[0m\n\u001b[1;32m 881\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0midx\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 882\u001b[0m \u001b[0midx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 883\u001b[0;31m \u001b[0maccount_for\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mipt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 884\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 885\u001b[0m \u001b[0;31m# add all variables that are true ancestors of the cost\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mRecursionError\u001b[0m: maximum recursion depth exceeded"
]
}
],
"source": [
"with gp_fit:\n",
" \n",
" gp_trace = sample(1000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"traceplot(gp_trace)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.1"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment