I am adapting and old two dimensional GMM to the new PyMC5
syntax, and I found an strange sampling error similar to this. Without really understanding what the error is about, I try the solution proposed there (to eliminate the shape argument that was redundant) but the error remains. Here goes the code snippet,
data generation,
import pymc as pm
import numpy as np
import matplotlib.pyplot as plt
import arviz as az
from pytensor.tensor.nlinalg import matrix_inverse
## Generate some data for three groups.
# Means, variances, covariances, and proportions
mus = np.array([[-4,2],[0,1],[6,-2]])
variance1 = [1,.4,1.5]
variance2 = [1,.8,5 ]
covariances = [.5,0, -1]
ps = np.array([0.2, 0.5, 0.3])
D = mus[0].shape[0]
# Total amount of data
N = 1000
# Number of groups
K = 3
# Form covariance matrix for each group
sigmas = [np.array([[variance1[i],covariances[i]],[covariances[i],variance2[i]]]) for i in range(K)]
# Form group assignments
zs = np.array([np.random.multinomial(1, ps) for _ in range(N)]).T
xs = [z[:, np.newaxis] * np.random.multivariate_normal(m, s, size=N) for z, m, s in zip(zs, mus, sigmas)]
# Stack data into single array
data = np.sum(np.dstack(xs), axis=2)
model construction,
testvals = [[-2,-2],[0,0],[2,2]]
# Model structure
with pm.Model() as mvgmm:
# Prior over component weights
p = pm.Dirichlet('p', a=np.array([1.]*K))
# Prior over component means
mus = [pm.MvNormal('mu_%d' % i, mu=pm.floatX(np.zeros(D)), tau=pm.floatX(0.1 * np.eye(D)), shape=(D,), initval=pm.floatX(testvals[i])) for i in range(K)]
# Cholesky decomposed LKJ prior over component covariance matrices
# chol, corr, stds = [pm.LKJCholeskyCov('packed_L_%d' % i, n=D, eta=2., sd_dist=pm.HalfCauchy.dist(1)) for i in range(K)]
L = [pm.LKJCholeskyCov('packed_L_%d' % i, n=D, eta=2., sd_dist=pm.HalfCauchy.dist(1))[0] for i in range(K)]
# Unpack packed_L into full array
# L = [pm.expand_packed_triangular(D, packed_L[i]) for i in range(K)]
# Convert L to sigma and tau for convenience
sigma = [pm.Deterministic('sigma_%d' % i ,L[i].dot(L[i].T)) for i in range(K)]
tau = [pm.Deterministic('tau_%d' % i,matrix_inverse(sigma[i])) for i in range(K)]
# Specify the likelihood
mvnl = [pm.MvNormal.dist(mu=mus[i],chol=L[i]) for i in range(K)]
Y_obs = pm.Mixture('Y_obs',w=p, comp_dists=mvnl, observed=data)
And this is the sampling,
with pm.Model() as mvgmm:
trace = pm.sample(chains=8)
And the error is
SamplingError: Cannot sample from the model, since the model does not contain any free variables.