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,
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.shape # 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)
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)) 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.