My model predicts correctly the p value for a geometric distribution i use for convolution.
BUT, the convergence statistics are bad / wrong and sampling takes forever unless i set tune/draws to 100.
I feel i miss something very basic here, maybe someone can point me into the correct direction:
The model
model = pm.Model()
base_vector = pt.as_tensor_variable(([1.0]+[0.0]*9))
# corresponds to a convolution of base_vector with pm.Geometric.dist(0.25)
geometric_convolution = pt.as_tensor_variable(
[0.26491849, 0.19868887, 0.14901665, 0.11176249, 0.08382187,
0.0628664 , 0.0471498, 0.03536235, 0.02652176, 0.01989132])
with model:
p= pm.Beta('p', alpha=1, beta=1)
distribution = pm.Geometric.dist(p)
convoluted_series = convolve1d(base_vector, distribution, 10)
sigma = pm.HalfNormal('sigma',0.1)
observed_convolution = pm.Normal('observed_convolution', mu=convoluted_series, sigma=sigma, observed=geometric_convolution)
trace = pm.sample(tune=100, draws=100)
The trace statistics:
p: mean=0.25 , sd=0 , r_hat = 1.43
sigma: mean=0, sd=0 , r_hat = 2.36
The convolution function:
import pytensor.tensor as pt
import pytensor.tensor.conv as conv
def convolve1d(x, distribution, length):
dist_pdf = pt.exp(pm.logp(distribution, pt.arange(1, length+1, 1)))
dist_pdf_normalized = dist_pdf / pt.sum(dist_pdf)
convolved = conv.conv2d(input=x[np.newaxis, np.newaxis,: ,np.newaxis],
filters=dist_pdf_normalized[np.newaxis, np.newaxis, :, np.newaxis],
border_mode='full')
return convolved[0, 0, :-length+1 ].T[0]