Interpretation of posterior predictive checks for a Gaussian Process

Questions v5
, ,
6

Hi Bill and others,

I was on a different project so I didn’t follow this for a while, but I wonder if you could offer any further help to resolve this issue with the following PPC plot (the code for this plot is shown below):

image
image1804×1414 164 KB

My updated code is shown below. I think that this model with MODEL = ‘GP_Latent_TS’ is the correct use of the GP model in PyMC.

Thank you.

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}))