

In [1]:

    
%load_ext autoreload
%autoreload 2
pylab.rcParams['figure.figsize'] = (10.0, 8.0)

Pseudo Marginal MCMC Example

This notebook gives an example of performing Poisson regression with a Gaussian Process Prior. The model for this example is \begin{align} \mathbf{f} = [f(t_1), \dots, f(t_T)] &\sim \mathcal{GP}(\mu_0, K_\theta) \\ \lambda &\equiv \exp( \mathbf{f} ) \\ y(t_i) &\sim \text{Pois}(\lambda(t_i)) && i \in \{1, \dots, T\} \end{align}

where $\mathbf{f}$ is some latent function evolving through time, $\theta$ are GP kernel hyper parameters that define it's 'smoothness', and $y_i \equiv y(t_i)$ are count observations at certain points in time (or corresponding to certain tiles in time, depending on the interpretation).

The rest of the notebook shows

generating a periodic 1-d 'ground truth' GP
Fitting a laplace approximation to the GP posterior CONDITIONED on knowing the hypers
Using the laplace approx as an internal mechanism to compute an unbiased approximation to the marginal likelihood
Performing pseudo-marginal MCMC with that marginal likelihood approximation
Comparing it to 'whitened' mh (ACF, ESS)



In [2]:

    
# Generate random data
np.random.seed(45)
from flecks.kernel import PerKernelUnscaled, SQEKernelUnscaled, LinearKernel
from numpy.linalg import inv, cholesky, solve
from scipy.stats import multivariate_normal

# generate poisson counts with a periodic kernel
sig_noise = 1e-6
pkern = PerKernelUnscaled()
skern = SQEKernelUnscaled()
lkern = LinearKernel()
N     = 200                      # sample size
T     = 50                       # number of time points
ell   = 1                        # 'true' length scale
per   = 20                       # 'true' period
tgrid = np.linspace(0, T, N)     # time grid
dt    = tgrid[1]-tgrid[0]        # tile length

# generate latent periodic GP
K     = pkern.K(tgrid, tgrid, hypers=[per, ell]) + sig_noise*np.eye(N) #+ .01*skern.K(tgrid,tgrid)
f_gt  = cholesky(K).dot(np.random.randn(N))
f0_gt = 3   # bias - total offset

# compute latent intensity function
lam_gt = np.exp(f_gt + f0_gt)
y      = np.random.poisson(lam=lam_gt) 

#visualize data
plt.plot(tgrid, y, 'ro')
plt.plot(tgrid, lam_gt)
plt.title("Synthetic Data and True Latent Function")









    Out[2]:





<matplotlib.text.Text at 0x10719f910>

Laplace Approximation for $p(\mathbf{f} | y, \theta)$



In [7]:

    
# Demonstrate laplace approximation to f_mu, f0, and f_Sig
from flecks.inference.laplace_approx import laplace_approx
from flecks.likelihoods.poisson import poisson_log_posterior

# Kernel integrates out f0 by adding 25*np.ones(N,N) - equivalent to f0 ~ N(0,5^2)
def kern_func(th): 
    return pkern.K(tgrid, tgrid, hypers=th) + 25*np.ones((N,N)) + 1e-6*np.eye(N) 
K_theta = kern_func([per, ell])

# define log like (implicitly carries around observations y)
#   args: 
#     f      = latent function (N length vector) (thing we want the laplace posterior for)
#     cK_inv = inverse of the cholesky decomposition of the gram matrix that governs f
#     grad   = bool to return gradient vector instead of likelihood
#     hessp  = bool to return the hessian product instead of the likelihood (needs p)
#     p      = length N vector to be multiplied by Hessian 
def log_like(f, cK_inv, grad=False, hessp=False, p=None): 
    return poisson_log_posterior(y=y, f=f, cK_inv=cK_inv, grad=grad, hessp=hessp, p=p)

# compute laplace approx - this optimizes the likelihood using Scipy, then computes the inverse hessian
f_mu, f_Sig, _, _ = laplace_approx(K_theta, log_like=log_like)

# plot samples, laplace posterior mean and samples
cf_Sig = cholesky(f_Sig)
#for n in range(10): 
#    fi = cf_Sig.dot(np.random.randn(N)) + f_mu
#    plt.plot(tgrid, np.exp(fi), c='grey')
pmu,   = plt.plot(tgrid, np.exp(f_mu), c='blue', linewidth=3)        # WHY DO I NEED A COMMA HERE?
ptrue, = plt.plot(tgrid, np.exp(f_gt+f0_gt), c='green', linewidth=3)
plt.plot(tgrid, np.exp(f_mu + 2*sqrt(np.diag(f_Sig))), c='grey', linewidth=2)
plt.plot(tgrid, np.exp(f_mu - 2*sqrt(np.diag(f_Sig))), c='grey', linewidth=2)

plt.legend([pmu, ptrue], ["$\hat \mathbf{f}$", "$\mathbf{f}_{gt}$"])
plt.title("Laplace Posterior and +/- 2 posterior stdevs")









    Out[7]:





<matplotlib.text.Text at 0x108bfb550>

Using the Laplace Approximation for an Unbiased Marginal Likelihood Estimate



In [8]:

    
from flecks.inference.pseudo_marginal import approx_log_marg_like

# test on above dataset
pers = np.sort(np.concatenate( [np.linspace(2*dt, 2*per, 40), np.arange(2,3*per+2)] ) )
log_pers = []
for i,per_i in enumerate(pers): 
    #approximate marginal likelihood, p(y | theta)
    if i%10==0: 
      print ".. computing %d of %d"%(i,len(pers))
    theta   = [per_i, ell]
    log_py  = approx_log_marg_like(kern_func(theta), log_like)
    log_pers.append(log_py)
    
# visualize likelihood function
plt.plot(pers, log_pers)
plt.vlines(per, ymin=np.min(log_pers), ymax=np.max(log_pers), colors='green', linewidth=5)
plt.title("Approximate log marginal likelihood, $\\log  p(Y | \\theta)$")









    



.. computing 0 of 100
.. computing 10 of 100
.. computing 20 of 100
.. computing 30 of 100
.. computing 40 of 100
.. computing 50 of 100
.. computing 60 of 100
.. computing 70 of 100
.. computing 80 of 100
.. computing 90 of 100






    Out[8]:





<matplotlib.text.Text at 0x108e497d0>

Pseudo Marginal MCMC over Periodic Kernel params with Laplace Approximation

Now actually sample periodic covariance hyperparameters given this approximation - simple pseudomarginal MH!



In [21]:

    
# Run pseudomarginal MH
from flecks.inference.pseudo_marginal import pseudo_gp_mh, jump_proposal_scalar

# starting state
th = np.array([3, 2.228*dt])

# track all samples
Nsamps   = 1000
th_samps = np.zeros((Nsamps, len(th)))
log_py_samps = np.zeros(Nsamps)

# define two simple proposals (will be randomly swapped)
def proposal_jump(th): 
    return np.array([jump_proposal_scalar(th[0]), th[1]])    

def proposal_noise(th, prop_scale = np.array([.1, .01])):
    return th + prop_scale*np.random.randn(2)

log_marg_like = None
num_jump_accepts = 0
num_noise_accepts = 0
for ii in range(Nsamps): 
    
    # report some stats for sanity
    if ii%10==0: 
        accept_ratio = float(num_jump_accepts + num_noise_accepts) / (ii+1)
        info = (ii, Nsamps, accept_ratio, num_jump_accepts, num_noise_accepts)
        print "sample %d of %d (accept ratio: %2.2f, %d jump accepts, %d noise accepts)"%info
        print th

    # randomly choose a proposal function
    if np.random.rand() < .5: 
        jumping = True
        gen_proposal = proposal_jump
    else:
        jumping = False
        gen_proposal = proposal_noise    
        
    # MH accept/reject
    th, accept, log_marg_like = pseudo_gp_mh(th, 
                                             kern_func     = kern_func, 
                                             like_func     = log_like, 
                                             ln_prior      = lambda(th): pkern.hyper_prior_lnpdf(th), 
                                             curr_marg_ll  = log_marg_like, 
                                             proposal_func = gen_proposal ) 
    th_samps[ii, ] = th
    log_py_samps[ii] = log_marg_like
    
    # track literally the minimal chain stats
    num_jump_accepts += jumping and accept
    num_noise_accepts += (not jumping) and accept









    



sample 0 of 1000 (accept ratio: 0.00, 0 jump accepts, 0 noise accepts)
[ 3.          0.55979899]
 th_per prop:  5.96621400103
 th_per prop:  11.9380603778
 th_per prop:  23.8645276856
 th_per prop:  47.6840046432
 th_per prop:  11.8554122487
sample 10 of 1000 (accept ratio: 0.36, 4 jump accepts, 0 noise accepts)
[ 47.68400464   0.55979899]
 th_per prop:  143.232745614
 th_per prop:  95.5739394135
 th_per prop:  191.202887145
 th_per prop:  475.301586915
sample 20 of 1000 (accept ratio: 0.48, 5 jump accepts, 5 noise accepts)
[ 95.05384842   0.37674559]
 th_per prop:  285.59447301
 th_per prop:  47.6517278542
 th_per prop:  475.897009767
 th_per prop:  285.508261322
 th_per prop:  47.5637991769
 th_per prop:  190.621323717
 th_per prop:  381.232440945
 th_per prop:  23.8479894984
sample 30 of 1000 (accept ratio: 0.39, 5 jump accepts, 7 noise accepts)
[ 95.30730809   0.35430797]
 th_per prop:  31.8081948277
 th_per prop:  287.503421765
 th_per prop:  383.269322729
 th_per prop:  191.561488499
 th_per prop:  47.9031082888
 th_per prop:  47.8861944997
sample 40 of 1000 (accept ratio: 0.37, 5 jump accepts, 10 noise accepts)
[ 95.81977128   0.27024632]
 th_per prop:  191.578640377
 th_per prop:  48.0126672211
 th_per prop:  31.9609256755
 th_per prop:  191.610231814
 th_per prop:  47.8422553316
sample 50 of 1000 (accept ratio: 0.29, 5 jump accepts, 10 noise accepts)
[ 95.81977128   0.27024632]
 th_per prop:  191.658090035
 th_per prop:  287.514489996
 th_per prop:  287.55518026
 th_per prop:  47.9136749926
 th_per prop:  48.0890125506
 th_per prop:  577.212038143
 th_per prop:  192.412121718
sample 60 of 1000 (accept ratio: 0.26, 5 jump accepts, 11 noise accepts)
[ 96.21255642   0.2140617 ]
 th_per prop:  19.1636764474
 th_per prop:  57.5593713233
 th_per prop:  9.81004607633
sample 70 of 1000 (accept ratio: 0.30, 6 jump accepts, 15 noise accepts)
[ 19.60828112   0.22283661]
 th_per prop:  39.4527592184
 th_per prop:  9.90991377202
 th_per prop:  4.93390275802
 th_per prop:  39.4337901046
 th_per prop:  138.104049907
 th_per prop:  39.4641561975
sample 80 of 1000 (accept ratio: 0.31, 6 jump accepts, 19 noise accepts)
[ 19.73417051   0.32037592]
 th_per prop:  39.5562495526
 th_per prop:  9.85247336432
 th_per prop:  6.58811218559
 th_per prop:  4.90692467434
 th_per prop:  39.3017851641
 th_per prop:  39.3079510971
sample 90 of 1000 (accept ratio: 0.29, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  58.8939400855
 th_per prop:  6.51730515314
 th_per prop:  6.55819177492
sample 100 of 1000 (accept ratio: 0.26, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  2.10532874813
 th_per prop:  9.83993571323
 th_per prop:  39.2530883997
 th_per prop:  78.5639935007
 th_per prop:  6.59306586035
sample 110 of 1000 (accept ratio: 0.23, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  6.48336477725
 th_per prop:  9.90081879187
 th_per prop:  39.3028551561
 th_per prop:  39.3025629718
 th_per prop:  39.2762832439
 th_per prop:  78.6451495983
 th_per prop:  4.89302823007
sample 120 of 1000 (accept ratio: 0.21, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  78.5296682023
 th_per prop:  4.96245650624
 th_per prop:  9.71347230042
 th_per prop:  9.83172863456
 th_per prop:  9.70735664884
 th_per prop:  9.82930200304
 th_per prop:  39.2634471325
 th_per prop:  58.9720055247
sample 130 of 1000 (accept ratio: 0.20, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  39.2539048616
 th_per prop:  6.481059909
 th_per prop:  3.9938010669
 th_per prop:  9.78117720192
 th_per prop:  6.56621954322
sample 140 of 1000 (accept ratio: 0.18, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  9.8352603657
 th_per prop:  4.90623146989
 th_per prop:  9.90161502435
 th_per prop:  9.81075128744
 th_per prop:  78.5915803151
sample 150 of 1000 (accept ratio: 0.17, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  9.93650619194
 th_per prop:  39.3313775557
 th_per prop:  39.2877404159
 th_per prop:  6.52793140063
 th_per prop:  78.5338178512
sample 160 of 1000 (accept ratio: 0.16, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  39.2933861164
 th_per prop:  6.52420150855
 th_per prop:  39.2204303981
 th_per prop:  9.90749214176
 th_per prop:  3.35473946449
sample 170 of 1000 (accept ratio: 0.15, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  58.9874408822
 th_per prop:  39.2918615466
 th_per prop:  98.1713583986
 th_per prop:  3.89931806117
 th_per prop:  58.886588502
 th_per prop:  58.9030240676
sample 180 of 1000 (accept ratio: 0.14, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  137.507742927
 th_per prop:  39.2702419566
 th_per prop:  39.3190168147
 th_per prop:  39.2661950671
 th_per prop:  9.89625756505
sample 190 of 1000 (accept ratio: 0.14, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  58.9818456367
 th_per prop:  39.3719765063
 th_per prop:  39.2660701886
 th_per prop:  6.59483269622
 th_per prop:  9.92663132775
 th_per prop:  9.84806042366
 th_per prop:  4.88722732005
sample 200 of 1000 (accept ratio: 0.13, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  3.27073769719
 th_per prop:  6.50113802526
 th_per prop:  58.9820188399
 th_per prop:  9.8034955963
 th_per prop:  9.83529345539
 th_per prop:  9.78716682808
 th_per prop:  9.8304006522
sample 210 of 1000 (accept ratio: 0.12, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  117.84863655
 th_per prop:  39.24311238
 th_per prop:  6.52164650746
 th_per prop:  39.2280248081
 th_per prop:  39.2957199682
 th_per prop:  9.84144169228
 th_per prop:  59.0163794981
sample 220 of 1000 (accept ratio: 0.12, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  9.82749483555
 th_per prop:  9.77150249651
 th_per prop:  39.2463589806
 th_per prop:  6.49601607538
 th_per prop:  39.3345850906
 th_per prop:  39.2991812992
 th_per prop:  9.80858109157
 th_per prop:  9.84580498043
sample 230 of 1000 (accept ratio: 0.11, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  9.83324129576
 th_per prop:  9.78445631812
 th_per prop:  9.7468025318
 th_per prop:  39.3337981793
 th_per prop:  4.9264709741
 th_per prop:  6.59396836931
sample 240 of 1000 (accept ratio: 0.11, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  117.785580194
 th_per prop:  9.81802577346
 th_per prop:  58.8711025865
 th_per prop:  39.3419687239
 th_per prop:  39.2376672818
 th_per prop:  39.2488950413
sample 250 of 1000 (accept ratio: 0.10, 6 jump accepts, 20 noise accepts)
[ 19.64190183   0.33722515]
 th_per prop:  3.29369362638
 th_per prop:  78.402793205
 th_per prop:  39.1901627284
 th_per prop:  59.684303233
 th_per prop:  39.7978369104
 th_per prop:  9.98378325267
sample 260 of 1000 (accept ratio: 0.10, 6 jump accepts, 21 noise accepts)
[ 19.89822017   0.44302682]
 th_per prop:  9.9296061018
 th_per prop:  39.8691894583
 th_per prop:  39.9902401951
 th_per prop:  6.63912329504
sample 270 of 1000 (accept ratio: 0.11, 6 jump accepts, 23 noise accepts)
[ 20.03008864   0.65144684]
 th_per prop:  6.6868471671
 th_per prop:  6.68370579393
sample 280 of 1000 (accept ratio: 0.11, 6 jump accepts, 26 noise accepts)
[ 20.05182478   0.77403235]
 th_per prop:  40.1079879252
 th_per prop:  39.6406038332
 th_per prop:  3.91293577035
 th_per prop:  5.04474915893
 th_per prop:  40.0198274313
sample 290 of 1000 (accept ratio: 0.12, 6 jump accepts, 28 noise accepts)
[ 20.04188099   0.76942636]
 th_per prop:  40.1352980664
 th_per prop:  80.2188897708
 th_per prop:  4.02084421149
 th_per prop:  9.98129030392
 th_per prop:  40.1278228657
 th_per prop:  100.167788488
 th_per prop:  80.2396042208
sample 300 of 1000 (accept ratio: 0.11, 6 jump accepts, 28 noise accepts)
[ 20.04188099   0.76942636]
 th_per prop:  5.00495567767
 th_per prop:  4.85632013459
 th_per prop:  4.98700554181
 th_per prop:  6.6849928991
 th_per prop:  40.1111950774
 th_per prop:  2.85195905482
 th_per prop:  40.0327480608
sample 310 of 1000 (accept ratio: 0.11, 6 jump accepts, 29 noise accepts)
[ 19.88560151   0.78333889]
 th_per prop:  9.98215186955
 th_per prop:  59.6871697121
 th_per prop:  9.9234471477
 th_per prop:  40.2824759633
sample 320 of 1000 (accept ratio: 0.12, 6 jump accepts, 32 noise accepts)
[ 20.18444776   0.74775885]
 th_per prop:  100.953623083
 th_per prop:  40.3664725041
 th_per prop:  10.0169701083
 th_per prop:  10.0402183095
 th_per prop:  10.0960769004
 th_per prop:  10.0419413965
sample 330 of 1000 (accept ratio: 0.11, 6 jump accepts, 32 noise accepts)
[ 20.18444776   0.74775885]
 th_per prop:  40.4623478585
 th_per prop:  100.294531465
 th_per prop:  3.31146427928
 th_per prop:  60.1929967433
 th_per prop:  40.2050576054
sample 340 of 1000 (accept ratio: 0.11, 6 jump accepts, 33 noise accepts)
[ 20.05894064   0.7460451 ]
 th_per prop:  80.3654722288
 th_per prop:  1.83074392686
 th_per prop:  6.70160299288
 th_per prop:  40.1067960333
 th_per prop:  60.2371055228
 th_per prop:  10.0771399377
sample 350 of 1000 (accept ratio: 0.11, 6 jump accepts, 33 noise accepts)
[ 20.05894064   0.7460451 ]
 th_per prop:  40.1335128226
 th_per prop:  40.1205643412
 th_per prop:  80.2121086402
 th_per prop:  60.1859785437
sample 360 of 1000 (accept ratio: 0.11, 6 jump accepts, 34 noise accepts)
[ 20.14860906   0.72178281]
 th_per prop:  5.03715983369
 th_per prop:  40.2676711808
 th_per prop:  10.0763162858
 th_per prop:  60.511233563
 th_per prop:  60.465210798
 th_per prop:  4.99860861732
sample 370 of 1000 (accept ratio: 0.11, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  4.98426770521
 th_per prop:  9.82525156366
 th_per prop:  39.6891452239
 th_per prop:  39.680828868
 th_per prop:  9.91056690177
 th_per prop:  39.7456853204
sample 380 of 1000 (accept ratio: 0.11, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  59.5167630926
 th_per prop:  9.98537245251
 th_per prop:  59.5024565271
 th_per prop:  5.00224703619
 th_per prop:  79.316457449
 th_per prop:  39.5997747566
 th_per prop:  9.80767641355
 th_per prop:  39.5906716037
sample 390 of 1000 (accept ratio: 0.10, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  39.6963745277
 th_per prop:  9.90411986606
 th_per prop:  39.6793536598
 th_per prop:  9.87071183682
sample 400 of 1000 (accept ratio: 0.10, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  39.6438352137
 th_per prop:  39.7418736189
 th_per prop:  59.4756935658
 th_per prop:  99.2377212142
 th_per prop:  9.90644081163
 th_per prop:  6.6540918483
 th_per prop:  9.92152959911
sample 410 of 1000 (accept ratio: 0.10, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  9.97686073393
 th_per prop:  79.4516621623
 th_per prop:  4.06694315944
 th_per prop:  9.97862552857
 th_per prop:  39.6522033125
sample 420 of 1000 (accept ratio: 0.10, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  39.6635243258
 th_per prop:  6.63942699434
 th_per prop:  6.61821587647
 th_per prop:  59.4541554999
 th_per prop:  39.6903172249
sample 430 of 1000 (accept ratio: 0.10, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  99.187536269
 th_per prop:  5.00633276717
 th_per prop:  6.59378220159
 th_per prop:  6.54371095014
sample 440 of 1000 (accept ratio: 0.09, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  99.1852460399
 th_per prop:  59.5756172253
 th_per prop:  6.58467054671
 th_per prop:  9.8865809907
sample 450 of 1000 (accept ratio: 0.09, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  6.63482388374
sample 460 of 1000 (accept ratio: 0.09, 6 jump accepts, 35 noise accepts)
[ 19.83780759   0.78564106]
 th_per prop:  9.88914026187
 th_per prop:  3.23057384504
 th_per prop:  79.6171670832
sample 470 of 1000 (accept ratio: 0.09, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.7043345747
 th_per prop:  4.96860726747
 th_per prop:  39.7674415777
 th_per prop:  6.5485363054
sample 480 of 1000 (accept ratio: 0.09, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  159.197160734
 th_per prop:  9.82973593637
 th_per prop:  79.6002916343
 th_per prop:  10.0129163445
 th_per prop:  9.96659547015
 th_per prop:  59.6782034054
 th_per prop:  9.97278308618
sample 490 of 1000 (accept ratio: 0.09, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.7932058985
 th_per prop:  59.6779706045
 th_per prop:  59.788009807
 th_per prop:  119.382762276
 th_per prop:  99.4198404438
 th_per prop:  4.964463096
 th_per prop:  39.7327797111
sample 500 of 1000 (accept ratio: 0.08, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  79.5831655795
 th_per prop:  39.8251460457
 th_per prop:  79.6686181463
 th_per prop:  6.71585882424
 th_per prop:  6.57558084481
 th_per prop:  39.8151962617
 th_per prop:  3.99688646389
 th_per prop:  9.95566620804
 th_per prop:  59.7115472696
sample 510 of 1000 (accept ratio: 0.08, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.85380212168
 th_per prop:  10.0479368557
 th_per prop:  119.412643799
 th_per prop:  9.88148741514
 th_per prop:  6.60343546617
 th_per prop:  9.90962561975
 th_per prop:  59.640977235
sample 520 of 1000 (accept ratio: 0.08, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  4.03057929815
 th_per prop:  4.9724341154
 th_per prop:  39.8076864931
 th_per prop:  5.05372999904
 th_per prop:  6.62668879785
 th_per prop:  39.8233644758
sample 530 of 1000 (accept ratio: 0.08, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.91215829625
 th_per prop:  59.7382747838
 th_per prop:  9.97616822787
sample 540 of 1000 (accept ratio: 0.08, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  10.027219109
 th_per prop:  9.89044569861
sample 550 of 1000 (accept ratio: 0.08, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  6.62682940619
 th_per prop:  59.7151692931
 th_per prop:  39.811877078
 th_per prop:  9.95082556082
 th_per prop:  39.8387136332
sample 560 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.6441479082
 th_per prop:  99.411590009
 th_per prop:  79.6089300464
 th_per prop:  6.72938341176
 th_per prop:  39.7890108578
sample 570 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.8166830842
 th_per prop:  39.8461814666
 th_per prop:  39.8745696315
 th_per prop:  4.96581575131
 th_per prop:  99.4832085598
 th_per prop:  39.7456009629
sample 580 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.94282162199
 th_per prop:  10.006429145
 th_per prop:  59.7255698151
sample 590 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  79.5973915312
 th_per prop:  6.66822644647
 th_per prop:  99.4596604275
 th_per prop:  79.5425320489
 th_per prop:  59.701866789
sample 600 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  79.5791121798
 th_per prop:  10.0085488143
 th_per prop:  6.63628867966
 th_per prop:  39.8483992205
 th_per prop:  39.818731553
 th_per prop:  59.6873411607
 th_per prop:  3.39825467451
sample 610 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  6.61240110037
 th_per prop:  2.44037527595
 th_per prop:  2.49592537412
 th_per prop:  99.506893798
 th_per prop:  39.824866753
sample 620 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.92920171665
 th_per prop:  9.86503785571
 th_per prop:  9.97761642151
 th_per prop:  59.7936297811
sample 630 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  6.58017945106
 th_per prop:  59.6654889879
 th_per prop:  59.7117434583
 th_per prop:  39.7332139101
sample 640 of 1000 (accept ratio: 0.07, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  4.97696306853
 th_per prop:  39.8734780005
 th_per prop:  9.94402403713
 th_per prop:  39.8439066195
sample 650 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.6404057135
 th_per prop:  9.90314477487
sample 660 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  159.197106171
 th_per prop:  9.91898434249
 th_per prop:  9.97439809541
 th_per prop:  6.59196810451
 th_per prop:  9.94934916232
sample 670 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.7332938717
 th_per prop:  79.6562756146
 th_per prop:  9.86904701572
sample 680 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.96987064374
 th_per prop:  9.975493712
 th_per prop:  39.8426515209
 th_per prop:  39.8174458554
 th_per prop:  39.793866044
 th_per prop:  79.6250891219
 th_per prop:  6.6851187626
sample 690 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.8068296375
 th_per prop:  59.6603519451
 th_per prop:  9.92433505949
 th_per prop:  5.05252328845
 th_per prop:  59.7003623434
 th_per prop:  39.8109004054
sample 700 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  79.5983882845
 th_per prop:  6.59774345963
 th_per prop:  39.7670819193
 th_per prop:  9.9313095807
 th_per prop:  59.7151376765
 th_per prop:  6.65895046661
 th_per prop:  6.68561428192
sample 710 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  99.5260673087
 th_per prop:  9.93398026961
 th_per prop:  59.6696838861
sample 720 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.91745450294
 th_per prop:  39.7659442264
 th_per prop:  6.70073714911
 th_per prop:  79.5996440317
 th_per prop:  59.6876932765
 th_per prop:  9.8996793081
 th_per prop:  3.36025521253
 th_per prop:  59.6665660033
sample 730 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.9326307838
 th_per prop:  39.8689018572
 th_per prop:  3.39390845943
sample 740 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.7461878083
 th_per prop:  10.0039099758
 th_per prop:  59.7132508469
 th_per prop:  59.7069086725
 th_per prop:  6.53317004181
 th_per prop:  9.97251181053
sample 750 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.8160187997
 th_per prop:  79.6366209467
 th_per prop:  59.6749697091
sample 760 of 1000 (accept ratio: 0.06, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.8277467482
 th_per prop:  6.66882111239
 th_per prop:  39.8468667107
 th_per prop:  79.6172915411
 th_per prop:  39.8025852242
 th_per prop:  39.8200859637
sample 770 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.6463454519
 th_per prop:  6.63284640191
 th_per prop:  39.7498824927
sample 780 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.6621985512
 th_per prop:  10.0620799217
 th_per prop:  39.7742182139
 th_per prop:  3.32955729661
 th_per prop:  39.8088688885
sample 790 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  10.0239562312
 th_per prop:  4.86184859213
 th_per prop:  10.0006894111
 th_per prop:  39.7229094003
sample 800 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  3.30388592819
 th_per prop:  39.8261490115
 th_per prop:  39.7977709793
 th_per prop:  39.7784426051
 th_per prop:  59.6513946237
sample 810 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  79.5785088148
 th_per prop:  4.98609385908
 th_per prop:  4.99333043826
sample 820 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  4.89872254558
 th_per prop:  9.94497214574
 th_per prop:  4.97100927263
 th_per prop:  5.02220647542
 th_per prop:  9.91900655979
 th_per prop:  39.7977710757
sample 830 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.6437307927
 th_per prop:  4.99461642394
 th_per prop:  6.67115380137
 th_per prop:  6.58707008184
sample 840 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  4.99340168793
 th_per prop:  4.96293012431
 th_per prop:  9.98810439643
 th_per prop:  39.7896351972
sample 850 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  79.6117327171
 th_per prop:  9.93247152527
 th_per prop:  9.91045858813
 th_per prop:  9.86623353514
 th_per prop:  79.6047206272
 th_per prop:  59.6729110183
sample 860 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.8299654528
 th_per prop:  5.03551196081
 th_per prop:  39.8186616123
 th_per prop:  4.85593416694
 th_per prop:  59.6609976855
 th_per prop:  39.7406075824
 th_per prop:  6.67000606793
 th_per prop:  4.0343606207
sample 870 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.8859222895
 th_per prop:  39.9196132555
 th_per prop:  39.7140826863
 th_per prop:  9.96910584581
 th_per prop:  9.81616538162
 th_per prop:  6.61288516664
sample 880 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.7702617832
 th_per prop:  4.02874953696
 th_per prop:  39.7554548836
 th_per prop:  9.93497883949
sample 890 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.95342848493
 th_per prop:  59.6859876742
 th_per prop:  39.7649754025
 th_per prop:  6.60975035352
 th_per prop:  39.8185881819
sample 900 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  6.66757610969
 th_per prop:  59.7250496075
 th_per prop:  39.8604202525
sample 910 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.6949757388
 th_per prop:  6.64067197507
 th_per prop:  79.535450403
 th_per prop:  59.7163972356
 th_per prop:  59.8411021351
 th_per prop:  9.98240781034
 th_per prop:  4.93262617363
 th_per prop:  3.90095345416
sample 920 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  4.99417061156
 th_per prop:  3.99002634438
 th_per prop:  10.0290275449
 th_per prop:  2.78501204118
 th_per prop:  9.9226700239
 th_per prop:  59.644930566
sample 930 of 1000 (accept ratio: 0.05, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.6902243186
 th_per prop:  6.57258546376
 th_per prop:  39.7542249912
 th_per prop:  39.6582718787
 th_per prop:  6.53983834404
 th_per prop:  59.6643830329
sample 940 of 1000 (accept ratio: 0.04, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  59.6855450127
 th_per prop:  39.8612305386
 th_per prop:  9.92802331291
 th_per prop:  59.7258609429
 th_per prop:  6.63272987343
 th_per prop:  9.90236619041
 th_per prop:  39.7994281012
 th_per prop:  59.6459073672
sample 950 of 1000 (accept ratio: 0.04, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  39.8124381006
 th_per prop:  4.97990012818
 th_per prop:  39.7864876844
 th_per prop:  39.892836454
 th_per prop:  39.7570487657
sample 960 of 1000 (accept ratio: 0.04, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.91801773205
 th_per prop:  39.7358597421
 th_per prop:  6.64807195808
 th_per prop:  39.6802022754
sample 970 of 1000 (accept ratio: 0.04, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  9.99057316633
 th_per prop:  6.63231717301
 th_per prop:  9.9927636588
 th_per prop:  9.97359278791
 th_per prop:  3.38319275743
 th_per prop:  9.97192095577
 th_per prop:  59.7327925451
sample 980 of 1000 (accept ratio: 0.04, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  10.0359000192
 th_per prop:  39.6999038826
 th_per prop:  10.0247400312
 th_per prop:  39.7804175541
 th_per prop:  5.0422912872
 th_per prop:  10.0166361293
sample 990 of 1000 (accept ratio: 0.04, 6 jump accepts, 36 noise accepts)
[ 19.90081217   0.88789212]
 th_per prop:  79.6297484859
 th_per prop:  6.54753016552
 th_per prop:  39.8013322617
 th_per prop:  59.7951351343
 th_per prop:  79.5399869817
 th_per prop:  39.8016097806
 th_per prop:  119.424353221



In [22]:

    
#plot traces and hists
plt.plot(th_samps[:,0])
plt.plot(th_samps[:,1])
plt.legend(["per trace", "length_scale trace"])
plt.title("MCMC traces")
plt.show()
plt.plot(log_py_samps)
plt.show()
plt.hist(th_samps[Nsamps/2:, 0], normed=True)
plt.hist(th_samps[Nsamps/2:, 1], normed=True)
plt.legend(["per", "ell"])
plt.title("MCMC Hists")
plt.show()



In [37]:

    
#produce posterior samples for each
from lgcp.mcmc.ess import elliptical_slice
plt.plot(tgrid, np.exp(f_gt+f0_gt))
plt.legend(["$\mathbf{f}_{gt}$"])
curr_f = f_mu

for ii in range(Nsamps/2, Nsamps, 10):
    K_theta  = pkern.K(tgrid, tgrid, [th_samps[ii], ell_samps[ii]]) + sig_noise*np.eye(N)
    cK_theta = cholesky(K_theta)
    cK_inv   = inv(cK_theta)
    for ni in range(30): 
        fi = cK_theta.dot(np.random.randn(N))
        curr_f, ll = elliptical_slice(initial_theta = curr_f, 
                                      prior         = fi, 
                                      lnpdf         = lambda(f): poisson_log_like(y,f,cK_inv) )   
    plt.plot(tgrid, np.exp(curr_f), c='grey')
plt.show()



In [627]:



In [23]:

    
from scipy.stats import gamma
xgrid = np.linspace(.01, 10, 100)
plt.plot(xgrid, gamma(2,scale=2).pdf(xgrid))









    Out[23]:





[<matplotlib.lines.Line2D at 0x111025cd0>]



In [27]:

    
plt.hist(np.random.randn(1000)*.1)









    Out[27]:





(array([   7.,   13.,   66.,  148.,  224.,  247.,  155.,   88.,   41.,   11.]),
 array([-0.30824535, -0.25016605, -0.19208674, -0.13400744, -0.07592814,
        -0.01784883,  0.04023047,  0.09830977,  0.15638908,  0.21446838,
         0.27254768]),
 <a list of 10 Patch objects>)



In [ ]:

    
#gen proposal, stop early if it's out of bounds
    jump_prop=False
    if np.random.rand() < .5: 
        jump_prop=True
        factor = np.random.geometric(p=.5) + 1  
        if np.random.rand() < .5: 
            factor = 1./factor
        #print factor
        th_prop = factor * th_curr + .1*np.random.randn()
    else: 
        th_prop = th_curr + .05*np.random.randn()
    ell_prop = ell_curr + .05*np.random.randn()
    log_prior = pkern.hyper_prior_lnpdf([th_prop, ell_prop])
    if log_prior < -1e20: 
        th_samps[ii] = th_curr
        ell_samps[ii] = ell_curr
        log_py_samps[ii] = log_py_curr
        continue 
    
    #compute log marginal like
    K_prop = pkern.K(tgrid, tgrid, [th_prop, ell_prop]) + sig_noise*np.eye(N) + 25*np.ones((N,N))
    log_py_prop = approx_log_marg_like(y, K_prop)
    print "period: ", th_prop
    print "length: ", ell_prop
    print log_py_prop - log_py_curr
    rat = log_py_prop - log_py_curr + \
          pkern.hyper_prior_lnpdf([th_prop, ell]) - pkern.hyper_prior_lnpdf([th_curr, ell])
    if np.log(np.random.rand()) < rat: 
        if jump_prop: 
            Njump_accepts += 1
        Naccept += 1 
        th_curr = th_prop
        ell_curr = ell_prop
        log_py_curr = log_py_prop
    
    #update state
    th_samps[ii] = th_curr
    ell_samps[ii] = ell_curr
    log_py_samps[ii] = log_py_curr