Getting wildly inaccurate MAP estimates in a Poisson Model

Questions version agnostic
1

I am trying to fit the following model:

\begin{align*} \alpha_1&\sim \mathcal N(0, 4^2)\\ \sigma &\sim \Gamma(0.5, 1)\\ \iota_0&\sim \mathcal N(1, 0.5^2)\\ \iota_n &\sim \mathcal N(\iota_{n-1}, \sigma^2)\\ \lambda_n &= \iota_n\exp(\alpha_1E_n)\\ g_n &\sim \mathrm{Poisson}(λ_n)\\ \end{align*}

where 50\ge n\ge 1. Here is my dataset:

import numpy as np
import matplotlib.pyplot as plt
plt.style.use("ggplot")

import theano
import pymc3 as pm
import arviz as az
az.style.use("arviz-darkgrid")

import scipy as sp
from scipy import stats

class FakeDataset:
    def __init__(self, dataset_size, alpha1, sigma):
        self.alpha1 = alpha1 
        self.sigma = sigma
        self.elo_diffs, self.iotas, self.lams, self.goals = self.generate_fake_dataset(dataset_size)

    def generate_fake_dataset(self, dataset_size):
        alpha1=self.alpha1
        sigma=self.sigma
        
        elo_diffs = np.random.choice(range(-100, 100, 20), dataset_size).astype(float)
        iotas = [1.]
        lams = []
        for i in range(dataset_size):
            iota_next = np.random.normal(iotas[-1], sigma)
            lam_next = iota_next * np.exp(elo_diffs[i] * alpha1)

            iotas.append(iota_next)
            lams.append(lam_next)
        
        iotas = np.array(iotas[1:])
        lams = np.array(lams)
        goals = np.random.poisson(lams).astype(float)
        return elo_diffs, iotas, lams, goals

dataset = FakeDataset(100, -1/200, 0.05)

print("ELO diff\tIota\tLambda\tGoals")
for i in range(15):
    print(f"{dataset.elo_diffs[i]}\t\t{dataset.iotas[i]:.3f}\t{dataset.lams[i]:.3f}\t{dataset.goals[i]}")

Here is my model:

np.random.seed(12)
model = pm.Model()
run_size = 50
with model:
    alpha1 = pm.Normal('alpha1', mu=0, sd=4)
    sigma = pm.Gamma('sigma', alpha=0.5, beta=1)
    iota0 = pm.Normal('iota0', mu=1, sd=0.5)
    iota = iota0
    for i in range(run_size):
        iota_i = pm.Normal(f'iota_{i}', mu=iota, sd=sigma)
        f = np.exp(alpha1 * dataset.elo_diffs[i])
        lam_i = pm.Deterministic(f"lam_{i}", iota_i * f)
        g_i = pm.Poisson(f"g_{i}", mu=lam_i, observed=dataset.goals[i])
        iota = iota_i
    
    map_est = pm.find_MAP()
    for i in range(run_size):
        print(f"MAP est: {map_est[f'lam_{i}']:.3f}",  f"true: {dataset.lams[i]:.3f}")

The results I am getting are completely off. I suspect this is probably because of this line

f = np.exp(alpha1 * dataset.elo_diffs[i])

Is there a way to reparametrize this model to fit better with PyMC3?

2

The parameterization np.exp(elo_diffs[i] * alpha1) is a little nonstandard because a product inside the exponentiation can lead to some very abrupt behavior. What if you parameterized as elo_diffs[i] + alpha1 instead? Is there a process-based reason that you want to use the multiplication?