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?