First attempt at a hierarchical model with custom likelihood

Questions
1

Hello, I am trying to implement a simple hierarchical model in pymc3. My data is a population of n=1\dots N
different time series. Each time series \mathbf F^{(n)} is modeled with a Gaussian likelihood with a diagonal covariance matrix (which is known), and the mean function is a non-linear function of 5 parameters:

\boldsymbol\theta^{(n)}=\{\Delta F^{(n)},F_b^{(n)},t_0^{(n)} t_\textrm{eff}^{(n)},t_E^{(n)}\}

The joint posterior for all time series parameters and the hyperparameters of the model is given by:

p\left(\{\boldsymbol\theta_n\},\mu,\sigma\lvert\{\mathbf F_n\}, \{\boldsymbol\sigma_n\}\right) \ \propto \prod_{n=1}^N\mathcal{L}\left(\mathbf F^{(n)},\boldsymbol\sigma^{(n)}\lvert \boldsymbol\theta^{(n)}\right)\,p\left(\boldsymbol\theta^{(n)}\lvert \mu, \sigma\right)\,p(\mu,\sigma)

where \{\sigma_n\} are the diagonal elements of the covariance matrix for the n-th time series, and

p\left(\boldsymbol\theta^{(n)}\lvert \mu, \sigma\right)= p\left(\Delta F^{(n)}\right)\,p\left(F_b^{(n)}\right)\, p\left(t_0^{(n)}\right) \,p\left(t_\textrm{eff}^{(n)}\lvert t_E^{(n)}\right)p\left(t_E^{(n)}\lvert \mu, \sigma\right)

I’ve attempted to implement this model in pymc3 as

import numpy as np
import pymc3 as pm
import theano
import theano.tensor as T

model = pm.Model()

with model:    
    n_population = len(lightcurves)
    
    BoundedNormal = pm.Bound(pm.Normal, lower=0.0) # Half-Normal dist

    # Hyperpriors
    mu_tE = BoundedNormal('mu_tE', mu=100, sd=200.)
    sig_tE = BoundedNormal('sig_tE', mu=20, sd=100.)
        
    # Priors 
    tE = pm.Lognormal('tE', mu=mu_tE, sd=sig_tE, shape=n_population)
    DeltaF = BoundedNormal('DeltaF', mu=[np.max(lc[:, 1]) for lc in lightcurves], 
                           sd=1., shape=n_population)
    Fb = pm.Normal('Fb', mu=0., sd=0.1, shape=n_population)
    t0 = pm.Uniform('t0', [np.min(lc[:, 0]) for lc in lightcurves], 
                    [np.max(lc[:, 0]) for lc in lightcurves], shape=n_population)
    # p(teff|tE)
    def conditional_prior(value):
        teff = T.cast(value, 'float64')
        return -T.log(tE) - (teff/tE)**2/1**2 
    
    teff = pm.DensityDist('teff', conditional_prior, shape=n_population)
    
    # Likelihood 
    def ln_likelihood(lightcurves):
        lnL = T.cast(0., 'float64')

        for n in range(n_population):
            lightcurve = T.cast(lightcurves[n], 'float64')
            t = lightcurve[:, 0]
            F = lightcurve[:, 1]
            sig_F = lightcurve[:, 2]
        
            # Mean model
            u0 = teff/tE
            u = T.sqrt(u0[n]**2 + ((t - t0[n])/tE[n])**2)
            A = lambda u: (u**2 + 2)/(u*T.sqrt(u**2 + 4))

            mu_F = DeltaF[n]*(A(u) - 1)/(A(u0[n]) - 1) + Fb[n]
            
            # Gaussian likelihood
            lnL += T.sum((-(mu_F - F)**2/sig_F**2))
        
        return lnL
    
    likelihood = pm.DensityDist('likelihood', ln_likelihood, observed={'lightcurves':lightcurves})

model.check_test_point()
model.vars

where lightcurves is a list containing N time series, each with a different length.

I have no idea if this is implemented correctly. I am getting an error: ValueError: setting an array element with a sequence..

2

You need to cast n within the forloop of ln_likelihood into theano using theano.shared (or there is a theano range as well you can try). But more generally I think you need to rewrite the ln_likelihood using theano.scan to make this efficient.