Understanding implementation of initial values in the usage with GMMs

Questions v5
1

Hi,

I am currently studying the usage of GMMs in pymc.

I wanted to get the log likelihood of a GMM model and had problems to retrieve it by how I implemented it.

Lets look into this very simple example, and this will work:


### some data:
rng = np.random.default_rng(123)
N = 1000
W = np.array([0.2, 0.5, 0.3])
MU = np.array([-1, 0, 1])
SIGMA = np.array([0.5, 1, 0.5])
component = rng.choice(MU.size, size=N, p=W)
x = rng.normal(MU[component], SIGMA[component], size=N)

### the model that works
K = 3
with pm.Model() as model:
    w = pm.Dirichlet("w", a=np.ones(K))
    mu = pm.Normal("mu", mu=np.linspace(x.min(), x.max(), K), sigma=10, shape=K, transform=pm.distributions.transforms.univariate_ordered)
    sigma = pm.HalfNormal("sigma", sigma=1,)
    y = pm.NormalMixture("y", w=w, mu=mu, sigma=sigma, observed=x)
    trace = pm.sample(idata_kwargs={"log_likelihood": True})
    pm.sample_posterior_predictive(trace, extend_inferencedata=True)

if I change mu to be:

mu = pm.Normal("mu", mu=0, sigma=10, shape=K, transform=pm.distributions.transforms.univariate_ordered, initval=np.linspace(x.min(), x.max(), K))

I will get the

NotImplementedError: Cannot convert models with non-default initial_values

What are the implications of the different use cases. Shouldn’t it be the same?

2

When do you get the error? Is it when computing log_likelihood?

3

Exactly :slight_smile:

4

In the meantime you can get around by not setting initval in the distribution, and instead define in initvals when calling pm.sample

5

I believe univariate_ordered is deprecated as well. You should use ordered (and upgrade to the latest version if you’re not getting a nagging message about it)

6

Thank you for the help.

I checked and adding initvals to sample solves the issue, using ordered is not working and the sampling seems not to converge to the same values for the different chains.

I am using pymc version 5.16.1

Thats the image for ordered

image
image950×526 104 KB

Thats the image for univariate_ordered

image
image950×526 75.3 KB

7

You are experiencing run-to-run variation. Here is the code for univariate_ordered:

    if name in ("univariate_ordered", "multivariate_ordered"):
        warnings.warn(f"{name} has been deprecated, use ordered instead.", FutureWarning)
        return ordered