Non-convergence of Conditional Autoregressive model

Questions
3

I already posted the whole definition (?) here you have the main equations of the CAR effect.

31
Screen Shot 2017-08-07 at 09.38.31.png713×151 9.79 KB

Here you find the new sources:

class CAR3(distribution.Continuous):
    """
    Conditional Autoregressive (CAR) with 2nd order neighbors
    """

    def __init__(self, adjacency, tau, rho, adjacency2, rho2, *args, **kwargs):
        super(CAR3, self).__init__(*args, **kwargs)
        n, m = adjacency.shape
        self.n = n

        adjacency_sparse = scipy.sparse.csr_matrix(adjacency)
        self.adjacency = theano.sparse.as_sparse_variable(adjacency_sparse)
        adjacency_sparse2 = scipy.sparse.csr_matrix(adjacency2)
        self.adjacency2 = theano.sparse.as_sparse_variable(adjacency_sparse2)

        self.n_neighbors = tt.as_tensor_variable(adjacency.sum(1))
        self.n_neighbors2 = tt.as_tensor_variable(adjacency2.sum(1))
        self.mean = tt.zeros(n)
        self.median = self.mean
        self.tau = tt.as_tensor_variable(tau)
        self.rho = tt.as_tensor_variable(rho)
        self.rho2 = tt.as_tensor_variable(rho2)

    def logp(self, x):
        priorvardenom = 1 - (self.rho + self.rho2) + self.rho * self.n_neighbors + self.rho2 * self.n_neighbors2
        priorvar = self.tau * priorvardenom
    
        Wx = theano.sparse.dot(self.adjacency, x.reshape((self.n, 1)))
        Wx2 = theano.sparse.dot(self.adjacency2, x.reshape((self.n, 1)))
        mu_w = (self.rho * tt.sum(Wx, axis=1) + self.rho2*tt.sum(Wx2, axis=1)) / priorvardenom
    
        return tt.sum(continuous.Normal.dist(mu=mu_w, tau=priorvar).logp(x))

data_df = pd.read_csv('dataset.csv')
adjacency = np.genfromtxt("adjacency.csv", delimiter=",")
adjacency2 = np.genfromtxt("adjacency2.csv", delimiter=",")

y = data_df['ncrimes'].values.astype('float32')
offset = data_df['population'].astype('float32')


with pm.Model() as pooled_model:
    b0 = pm.Normal('intercept', mu=0, tau=0.001)

    rho = pm.Uniform('rho', 0, 1)
    rho2 = pm.Uniform('rho2', 0, 1)
    tau = pm.HalfCauchy('tau', 2.5)
    mu_phi = CAR3('mu_phi', adjacency=adjacency, rho=rho, adjacency2=adjacency2, rho2=rho2, tau=tau, shape=len(y))
    phi2 = pm.Deterministic('phi', mu_phi - tt.mean(mu_phi))

    mu = tt.exp(b0 + phi2  + tt.log(offset))

    ncrimes = pm.Poisson('ncrimes', mu=mu, observed=y)

    pooled_trace = pm.sample(3000, njobs=3, tune=5000)  # , , tune=1000

    plot_traces(pooled_trace)

Thanks

adjacency2.csv (14.1 KB)
adjacency.csv (14.1 KB)
dataset.csv (608 Bytes)