I am doing a generalized linear regression with Negative Binomial likelihood. This problem is different from a typical linear model where the x and y values are separate from each other. I am looking at timeseries data, where for each values in the series, I regress on the previous W values. I would like to define a custom distribution that can take in an array of counts, perform this manipulation, and determine the likelihood. This is my code:
class MyNegBin(pm.Discrete): def __init__(self, W, intercept, beta, beta_recent, alpha, *args, **kwargs): super(MyNegBin, self).__init__(*args, **kwargs) self.W = W self.intercept = intercept self.beta = beta self.beta_recent = beta_recent self.alpha = alpha def logp(self, counts): W = self.W N = len(counts) - W x = np.vstack([counts[i:i+W] for i in range(N)]) y = counts[W:] mu = self.intercept for i in range(W-1): mu += self.beta[i] * x[:,i] mu += self.beta_recent * x[:,W-1] lik = pm.NegativeBinomial.dist(mu=mu, alpha=self.alpha) return lik.logp(y) with pm.Model() as model: intercept = pm.Normal('intercept', mu=0, sd=10) beta = pm.Normal('beta', mu=0, sd=1, shape=W-1) beta_recent = pm.Normal('beta_recent', mu=1, sd=1) alpha = pm.TruncatedNormal('alpha', mu=2000, sigma=500, lower=0) Y_obs = MyNegBin('Y_obs', W, intercept, beta, beta_recent, alpha, observed=counts)`
The problem is that logp() is expecting univariate distribution, so it’s required to take in a constant rather than an array. Is there any way to work around this?