Interpretation of posterior predictive checks for a Gaussian Process

Questions v5
, ,
10

Dear Jesse,

Please accept my apologies for the late reply!

Due to a combination of responsibilities across different projects, coupled with my vacation, it took me a while to return to pymc.

I must say, your response is one of the best answers I have received from any online community. Your example was particularly helpful in aiding my understanding of the process.

Based on your advice, I intend to use the posterior values for Lambda simply.

Regarding your comment “I’d show an example with your code but it’s not copy-pastable,” just out of curiosity, if possible, I’d be interested in seeing an example from my code. I’m including (copy and paste from my previous post without styling) it below for your reference (if you are referring to a different code of mine, please let me know).

Again, thank you so much for your fantastic help!

– Seongeun

LIKELIHOOD = ‘truncated_normal’
MODEL = ‘GP_Latent_TS’

with pm.Model() as model3:

if LIKELIHOOD == 'lognormal':
    Lambda = pm.LogNormal("Lambda", mu_logx, sigma_logx, shape=m)
    # Truncated not working ...
    #Lambda = pm.Truncated("Lambda", pm.LogNormal.dist(mu = mu_logx, sigma = sigma_logx), shape=m, upper=10)
elif LIKELIHOOD == 'truncated_normal':
    normal_dist = pm.Normal.dist(mu=mu_x, sigma=sigma_x)
    Lambda = pm.Truncated("Lambda", normal_dist, shape=m, lower=0, upper=10)
    #Lambda = pm.Normal("Lambda", 1, 0.5, shape=m)

#===============================================================================
# Mean function
#===============================================================================
if MODEL == 'MvNormal':
    print ('>>>>> Using my own Cov model ...')
    mu = pm.math.dot(np.array(K, dtype=np.float32), Lambda)
else:
    print ('>>>>> Using GP model ...')
    if MODEL in ['GP', 'GP_Latent_TS']:
        mu = LinearMean(K = K, Lambda = Lambda)

# Stat Rethinking, P. 470
eta_sq = pm.Exponential("eta_sq", 2)
#rho_sq = pm.Normal("rho_sq", 3.0, 0.25)
rho_sq = pm.Exponential("rho_sq", 0.5)

if SIGMA_MODEL == 'halfcauchy':
    sigma = pm.HalfCauchy ("sigma", 1)
else:
    sigma = pm.Exponential("sigma", 2.0)

#===============================================================================
# specify the covariance function:
#===============================================================================
if MODEL == 'MvNormal':
    cov = cov_fun_temporal (eta_sq, rho_sq, sigma)
else:
    if MODEL == 'GP':
        cov = eta_sq * pm.gp.cov.Exponential(input_dim=len (Y), ls=rho_sq)
    else:
        cov = eta_sq * pm.gp.cov.Exponential(input_dim=1, ls=rho_sq)
        #cov = eta_sq * pm.gp.cov.ExpQuad(input_dim=1, ls=rho_sq)
        # For stability
        cov += pm.gp.cov.WhiteNoise(1e-6)
        
#===============================================================================
# specify the GP:
#===============================================================================
if MODEL == "GP":
    gp = pm.gp.Marginal(mean_func=mu, cov_func=cov)
elif MODEL == 'GP_Latent':
    gp = pm.gp.Latent(cov_func=cov)
    # @Dmat: temporal distance matrix 
    f = gp.prior("f", X=Dmat)
elif MODEL == "GP_Latent_TS":
    gp = pm.gp.Latent(mean_func=mu, cov_func=cov)
    f = gp.prior("f", X=X)

#===============================================================================
# Likelihood
#===============================================================================
if MODEL == 'MvNormal':
    print ('>>>>> Using MvNormal specifiction ...')
    Y_ = pm.MvNormal ("Y_", mu=mu, cov=cov, observed=Y)
elif MODEL == 'GP':    
    print ('>>>>> Using GP marginal model specifiction ...')
    Y_ = gp.marginal_likelihood("Y_", X=Dmat, y=Y, sigma = sigma)
elif MODEL == 'GP_Latent':
    mu = pm.math.dot(np.array(K, dtype=np.float32), Lambda)
    Y_ = pm.Normal("Y_", mu=mu, sigma = sigma, observed=Y)
elif MODEL == 'GP_Latent_TS':
    if True:
        Y_ = pm.Normal("Y_", mu=f, sigma = sigma, observed=Y)
    else:
        SIGMA = np.diag (np.ones (len(Y))) * sigma**2
        Y_ = pm.MvNormal("Y_", mu=f, cov = SIGMA, observed=Y)
   
if MODEL in ['GP_Latent', 'GP_Latent_TS']:
    trace = pm.sample(NUM_SAMPLE, tune=TUNE, target_accept=0.9, random_seed=RANDOM_SEED,\
                    idata_kwargs = {'log_likelihood': True})
else:
    trace = pm.sample_prior_predictive()

    trace.extend(pm.sample(NUM_SAMPLE, tune=TUNE, target_accept=0.9, random_seed=RANDOM_SEED,\
                    idata_kwargs = {'log_likelihood': True}))