Hi all,
Following the remark of @jessegrabowski, I checked if .dist nodes are updated. But it looks like they are not, and I am stuck on how to proceed.
Specifically, I am motivated by the following problem. Suppose that I want distribute N values according to some distribution with the mean \mu and standard deviation \sigma. Then I count how many of them fall in given bins from 0 to D and repeat such experiment T times. I arrive at multinomial likelihood, from which I want to recover the parameters (\mu, \sigma).
I wrote the following code, but it looks like it does not any inference:
import pymc as pm
import aesara as ae
import aesara.tensor as at
import aesara.tensor.math as tt
with pm.Model() as model:
# Data
T = tt.as_tensor_variable(counts.shape[0], dtype="int64")
D = tt.as_tensor_variable(counts.shape[1] - 1, dtype="int64")
counts_matrix = tt.as_tensor_variable(counts, dtype="int64")
# Priors
μ = pm.HalfNormal('μ', 5.0, initval=2.0)
σ = pm.HalfNormal('σ', 5.0, initval=2.0)
# Likelihood
### Parameters of the lognormal distribution
σlog = tt.sqrt(tt.log((σ / μ)**2 + 1))
μlog = tt.log(μ) - σlog**2 / 2
π, _ = ae.scan(lambda t: tt.exp(pm.logcdf(pm.Lognormal.dist(mu = μlog, sigma = σlog), t + 0.5)) - tt.exp(pm.logcdf(pm.Lognormal.dist(mu = μlog, sigma = σlog), t - 0.5)),
sequences = [tt.arange(0, D + 1)])
llks, _ = ae.scan(lambda t, X: pm.logp(pm.Multinomial.dist(n = tt.sum(X[t, :]), p = π / tt.sum(π)), X[t, :]),
sequences = [tt.arange(0, T)], non_sequences=[counts_matrix])
pm.Potential('llk', llks)
trace = pm.sample(tune=1000, chains=4, init='advi')
the gist can be found here 20221015-pymc-question.ipynb · GitHub. You can see that there is no result at the end.
I have done something quite similar using Stan and it worked there, but I have difficulties implementing this in pymc.
functions {
/* discretized version of gamma, Weibull, and lognormal distributions */
vector dPMF(real param1, real param2, int K) {
vector[K] res;
for (k in 1:K)
res[k] = lognormal_cdf(k - 0.5 | param1, param2);
return append_row(res[1] / res[K], (tail(res, K-1) - head(res, K-1)) / res[K]);
}
}
data {
int<lower = 1> D;
int<lower = 1> T;
array[D + 1, T] int<lower=0> n; // counts
}
transformed data {
array[T] int<lower = 0> n_sum;
for (t in 1:T)
n_sum[t] = sum(n[, t]);
real jitter = 1e-9;
}
parameters {
real<lower=0> mu, sigma;
}
transformed parameters {
real param1, param2,
lps = 0.0;
{
// Lognormal distribution
param2 = sqrt(log(square(sigma / mu) + 1.0));
param1 = log(mu) - square(param2) / 2.0;
vector[T+1] beta = dPMF(param1, param2, T+1);
for (t in 1:T-1) if (n_sum[t] > 0) {
vector[D+1] beta_current;
beta_current[1:D] = head(beta, D);
beta_current[D+1] = max({1.0 - sum(head(beta_current, D)), jitter});
beta_current /= sum(beta_current);
lps += multinomial_lpmf(n[, t] | beta_current);
}
}
}
model {
/* priors */
sigma ~ cauchy(0.0, 5.0);
mu ~ std_normal();
/* likelihood */
target += lps;
}
Many thanks in advance!
Upd1: When I do the same simulation in pymc3, the results are correct. Please see 20221015-pymc3-version.ipynb · GitHub
