4

I'm still a noob in PyMC3, so the question might me naive, but I don't know how to translate this pymc2 code in pymc3. In particular it's not clear to me how to translate the R function.

beta = pymc.Normal('beta', mu=0, tau=1.0e-4)
s = pymc.Uniform('s', lower=0, upper=1.0e+4)
tau = pymc.Lambda('tau', lambda s=s: s**(-2))

### Intrinsic CAR
@pymc.stochastic
def R(tau=tau, value=np.zeros(N)):
    # Calculate mu based on average of neighbors
    mu = np.array([sum(W[i]*value[A[i]])/Wplus[i] for i in xrange(N)])

    # Scale precision to the number of neighbors
    taux = tau*Wplus
    return pymc.normal_like(value, mu, taux)

@pymc.deterministic
def M(beta=beta, R=R):
    return [np.exp(beta + R[i]) for i in xrange(N)]

obsvd = pymc.Poisson("obsvd", mu=M, value=Y, observed=True)
model = pymc.Model([s, beta, obsvd])

Code from https://github.com/Youki/statistical-modeling-for-data-analysis-with-python/blob/945c13549a872d869e33bc48082c42efc022a07b/Chapter11/Chapter11.rst, and http://glau.ca/?p=340

Can you help me? Thanks

merv
  • 67,214
  • 13
  • 180
  • 245
marcodena
  • 550
  • 2
  • 5
  • 15
  • The pymc3 doc is indeed a pain... I spent and lost days to adapt some codes from pymc2 to pymc3... are u still interested ? – Stéphane Aug 04 '17 at 08:48
  • 1
    Yes, I have the same concern. I don't want to start another question, since this *exactly* captures my concern. So I am instead offering up a bounty. My specific issue: I don't understand how to convert `@deterministic` functions which can have returned results auto-correlated (like the `[np.exp(beta + R[i]) for i in xrange(N)]`) into the PyMC3 equivalent. How can `def M` be written in PyMC3? – Mike Williamson Nov 20 '19 at 17:13
  • 1
    @MikeWilliamson There are implementations of CAR in PyMC documentations. I added one here. – Nathan Nov 26 '19 at 20:12
  • 1
    @Nathan Thank you!! Unfortunately, my bounty expired... I had to work on another task at work and could not revisit this. I wish I could award you the bounty. :( – Mike Williamson Dec 11 '19 at 10:39
  • @MikeWilliamson that’s fine! I’m glad your problem is solved. – Nathan Dec 11 '19 at 13:22

1 Answers1

2

In PyMC3, you can implement the CAR model using the scan function of Theano. There is a sample code in their documentation. There are two implementations for CAR in the linked document. Here is the first one [Source]:


from theano import scan
floatX = "float32"

from pymc3.distributions import continuous
from pymc3.distributions import distribution

class CAR(distribution.Continuous):
    """
    Conditional Autoregressive (CAR) distribution

    Parameters
    ----------
    a : list of adjacency information
    w : list of weight information
    tau : precision at each location
    """
    def __init__(self, w, a, tau, *args, **kwargs):
        super(CAR, self).__init__(*args, **kwargs)
        self.a = a = tt.as_tensor_variable(a)
        self.w = w = tt.as_tensor_variable(w)
        self.tau = tau*tt.sum(w, axis=1)
        self.mode = 0.

    def get_mu(self, x):

        def weigth_mu(w, a):
            a1 = tt.cast(a, 'int32')
            return tt.sum(w*x[a1])/tt.sum(w)

        mu_w, _ = scan(fn=weigth_mu,
                       sequences=[self.w, self.a])

        return mu_w

    def logp(self, x):
        mu_w = self.get_mu(x)
        tau = self.tau
        return tt.sum(continuous.Normal.dist(mu=mu_w, tau=tau).logp(x))

with pm.Model() as model1:
    # Vague prior on intercept
    beta0 = pm.Normal('beta0', mu=0.0, tau=1.0e-5)
    # Vague prior on covariate effect
    beta1 = pm.Normal('beta1', mu=0.0, tau=1.0e-5)

    # Random effects (hierarchial) prior
    tau_h = pm.Gamma('tau_h', alpha=3.2761, beta=1.81)
    # Spatial clustering prior
    tau_c = pm.Gamma('tau_c', alpha=1.0, beta=1.0)

    # Regional random effects
    theta = pm.Normal('theta', mu=0.0, tau=tau_h, shape=N)
    mu_phi = CAR('mu_phi', w=wmat, a=amat, tau=tau_c, shape=N)

    # Zero-centre phi
    phi = pm.Deterministic('phi', mu_phi-tt.mean(mu_phi))

    # Mean model
    mu = pm.Deterministic('mu', tt.exp(logE + beta0 + beta1*aff + theta + phi))

    # Likelihood
    Yi = pm.Poisson('Yi', mu=mu, observed=O)

    # Marginal SD of heterogeniety effects
    sd_h = pm.Deterministic('sd_h', tt.std(theta))
    # Marginal SD of clustering (spatial) effects
    sd_c = pm.Deterministic('sd_c', tt.std(phi))
    # Proportion sptial variance
    alpha = pm.Deterministic('alpha', sd_c/(sd_h+sd_c))

    trace1 = pm.sample(1000, tune=500, cores=4,
                       init='advi',
                       nuts_kwargs={"target_accept":0.9,
                                    "max_treedepth": 15})

The M function is written here as:

mu = pm.Deterministic('mu', tt.exp(logE + beta0 + beta1*aff + theta + phi))
Nathan
  • 471
  • 2
  • 12