Including Potentials of Neighbours

Questions
1

I’m trying to model a 2D field of likelihoods of some occurance of X event, which is dependent on a value in the 2D space, and there are also observations of the event. I’d also like the likelihood to be dependent on the likelihoods of the neighbouring distributions. Here is a minimal example:

import numpy as np
import pymc3 as pm
import matplotlib.pyplot as plt

#data
x = np.zeros((10, 10))
#observations
observations= np.zeros((10, 10))

indexgen =pm.DiscreteUniform.dist(lower=0, upper=9)
xlocs  = indexgen.random(size=20)
ylocs= indexgen.random(size=20)
print(xlocs)
sense =pm.Bernoulli.dist(p=0.2)
isObserved = sense.random(size=20)
print(isObserved)
for i, obs  in enumerate(isObserved):
    x[ylocs[i], xlocs[i]]=1
    if obs:
        observations[ylocs[i], xlocs[i]]=1

#add some random observations
xlocs  = indexgen.random(size=4)
ylocs= indexgen.random(size=4)
observations[ylocs, xlocs]=1

fig=plt.figure(figsize=(12, 12))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.set_title('Dataset')
ax1.imshow(x)
ax2.set_title('Observations')
ax2.imshow(observations)
plt.show()
with pm.Model() as model:
    sd = pm.Exponential('sd',lam=3.0)
    mean = pm.Normal('mu' ,mu=0., sd=10.)
    
    theta = pm.Normal('theta', mu=0.01, sd=sd, shape=2)
    
    link = pm.Deterministic('link', 1/(1 - pm.math.exp(-(theta[0] + theta[1]*x))) )
    
    prob = pm.Bernoulli('Prob', p=link, shape=(10, 10), observed=observations)
    pm.Potential('neighbours', prob[i-1]+prob[i+1])
    trace = pm.sample(1000, step=pm.Metropolis(), init='adapt_diag')

I cant get the indexing to work and I get out of bounds errors. How would I add padding to the pm.Potential or limit the application of the pm.Potential to only an area where it is valid.

A second question is how do I index the prob distribution in the y direction.

2

I’d like to learn a bit more about the model you’re working on. It looks like you’re making the following assumptions:

  • You have binary-valued data on a 2D grid which is modeled using a combination of logistic regression and cross-neighbor correlation
  • You want each grid cell’s observation to be correlated with the immediate neighbors.

This sounds quite a bit like the autologistic model which (I think) is synonymous with binary conditional autoregression. Are these terms relevant to your problem?

3

The Autologistic is not a term I’ve heard except in terms of autoregressive neural models. The most appropriate one I think that is relevant is Markov Random Fields. The way I like to think of it is more of a binary segmentation problem where only some examples are observed. I have implemented a version using a binary glm without the condition on neighbours. Which worked fairly well as there was no real spatial information encoded. e.g.

with pm.Model() as model:
    sd = pm.Exponential('sd',lam=3.0)
    mean = pm.Normal('mu' ,mu=0., sd=10.)
    
    data={
    'x' : x.flatten(),
    'y' : observations.flatten()
    }
    prior= {
        'Intercept' : pm.Normal.dist(mu=mean, sd=sd),
        'x' : pm.Normal.dist(mu=mean, sd=sd)
    }
    pm.GLM.from_formula('y ~ x', data, family=pm.glm.families.Binomial(), priors=prior)
    trace = pm.sample(1000,init='adapt_diag')

This results in a very stark classification, spatially.
I was trying to implement something that now included spatial information, such as the neighbouring potentials so that neighbouring pixels would be a bit more likely to be uncertain as to the classification.

4

I see. I can offer a more complicated solution to your problem, but I am definitely interested to hear of a potentially simpler answer if anyone has one.

The MRF you are referring to is the same thing as the terms I mentioned. There’s a very good tutorial on implementing a model with an MRF (aka CAR) component here but it’s not trivial to understand/use if you have no exposure to spatial statistics. The good news is that the MRF model from that link uses some tricks to make sampling much faster than if you tried to construct a Potential term with some ad hoc adjustments.

The challenge is most Bayesian models for gridded data allowing for spatial correlation can’t be specified in terms of sequential correlations (i.e. a directed graphical model) but instead have to be an undirected graphical model requiring definition of the joint likelihood over all sites simultaneously. This is unfortunately more complicated than simply specifying sequential correlations. You can think of the conditional autoregression (CAR) model from the above link as a spatially correlated intercept to be added to the other regression component. The CAR model is also a special case of a Gaussian process (GP) for which PyMC3 has very robust modeling tools.

If you want a quick sketch of what the CAR or GP model might look like then I am happy to show that.

2 Likes
5

Okay, I wasn’t aware of CAR and had no understanding of Gaussian Processes, but CAR and GPs seem to fit exactly the type of models I am interested in. I’ll do some more research and post back a solution for critique and sharing if I get something working. Thanks