Hello, I want to migrate a GLM model to Rolling Regression Model, but I’m encountering some issues.
I don’t know how to keep all factors non-negative in MvNormal. My current approach is to use pm.math.exp . This actually changes the distribution of the sample, and I’m not sure if it’s a good idea. At the same time, pm.Truncated is not implemented for multivariate distributions, and MvHalfNormal seems to be unavailable.
Here is my current code, and I’m hoping someone could help me.
signal_num = data.shape[1]-1
COORDS = {'coeffs':data.iloc[:,:-1].columns.tolist(),'time':data.index.values}
model = pm.Model(coords = COORDS)
with model:
sigma = pm.HalfNormal("sigma", sigma=0.1)
X = pm.MutableData("X",data.iloc[:,:-1].values,dims=('time','coeffs'))
beta_sigma = pm.Exponential("beta_sigma", 10.0,dims='coeffs')
beta = pm.MvGaussianRandomWalk("beta",
mu=np.zeros(signal_num),
dims=("time","coeffs"),
cov=pytensor.tensor.diag(beta_sigma),
init_dist=pm.MvNormal.dist(0.1, 0.01*np.eye(signal_num)))
Y = pm.MutableData("Y",data['Y'].values)
mu = pm.math.sum(X*pm.math.exp(beta),axis=1)
Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
idata = jax.sample_numpyro_nuts(chains=1,draws=1000,tune=1000,postprocessing_backend="cpu",)
My original GLM model is
COORDS = {'coeffs':data.iloc[:,:-1].columns.tolist()}
model = pm.Model(coords = COORDS)
with model:
sigma = pm.HalfNormal("sigma", sigma=0.1)
X = pm.MutableData("X",data.iloc[:,:-1].values)
Y = pm.MutableData("Y",data['Y'].values)
beta = pm.TruncatedNormal(name = 'beta',mu=0.1,sigma=0.1,lower=0,dims="coeffs")
mu = pm.math.dot(X,beta)
Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
idata = jax.sample_numpyro_nuts(chains=4,draws=10000,tune=10000)