Vanilla implementation of Gaussian Process in PyMC3?

Questions
13

Hey Jesse, I reran the code code by manually defining the mean function and the covariance function, but for some reason, I kept running into this error:

ValueError: The input matrix must be symmetric positive semidefinite.
Apply node that caused the error: multivariate_normal_rv{"(n),(n,n)->(n)"}(RNG(<Generator(PCG64) at 0x31AB40740>), [], [0. 0. 0. ... 0. 0. 0.], Composite{((i2 * exp((-1.0 * i0 * i1))) + i3)}.0)
Toposort index: 5
Inputs types: [RandomGeneratorType, TensorType(int64, shape=(0,)), TensorType(float64, shape=(10,)), TensorType(float64, shape=(10, 10))]
Inputs shapes: ['No shapes', (0,), (10,), (10, 10)]
Inputs strides: ['No strides', (0,), (8,), (80, 8)]
Inputs values: [Generator(PCG64) at 0x31AB40740, array([], dtype=int64), 'not shown', 'not shown']
Outputs clients: [[output[7](multivariate_normal_rv{"(n),(n,n)->(n)"}.0)], [output[1](alpha), Composite{exp((i0 + i1))}(ExpandDims{axis=0}.0, alpha)]]

The only way to sample without running into errors was to increase the amount of jitter in the covariance matrix:

cov = (eta_squared * pm.math.exp( -(rho_squared * ISLAND_DISTANCES ** 2) ) + (pt.eye(ISLAND_COUNT) * 1e-2))

Here’s my final model definition (for the section that incorporates population size):

with pm.Model(coords=coords) as distance_population_model:

    population = pm.MutableData("log_population", KLINE.logpop.values.astype(float))
    CULTURE_ID_ = pm.MutableData("CULTURE_ID", CULTURE_ID)

    # Priors
    alpha_bar = pm.Exponential("alpha_bar", 1)
    gamma = pm.Exponential("gamma", 1)
    beta = pm.Exponential("beta", 1)
    eta_squared = pm.Exponential("eta_squared", 2)
    rho_squared = pm.Exponential("rho_squared", 0.5)

    # Gaussian Process
    mu = pm.math.zeros(ISLAND_COUNT)
    cov = (eta_squared * pm.math.exp( -(rho_squared * ISLAND_DISTANCES ** 2) )
           + (pt.eye(ISLAND_COUNT) * 1e-2))
    alpha = pm.MvNormal('alpha', mu = mu, cov = cov, dims = 'culture')

    # Likelihood
    lambda_T = (alpha_bar / gamma * population[CULTURE_ID_] ** beta) * pm.math.exp(
        alpha[CULTURE_ID_]
    )
    pm.Poisson("T", lambda_T, observed=TOOLS, dims="culture")
  1. Do you know why I’m getting the error?
  2. Is it okay for me to fix the error by adding any constant c to each element in the diagonal of the covariance matrix?