Hi there,

I am trying to fit a mixture of Gaussians to my data. I got two questions:

- Is the approach I am using the right one?
- How can assign membership probabilities to each data point? I mean, I would like to know the probability of each data point belonging to one of the Gaussians.

Here is my code.

```
def model(data):
w = jnp.array([0.7, 0.3])
weights = numpyro.sample("weights", dist.Dirichlet(w))
mu_main = numpyro.sample("mu_main", dist.Normal(0, 5))
std_main = numpyro.sample("std_main", dist.Uniform(0, 20)) # consider using dist.LogNormal
mu_side = numpyro.sample("mu_side", dist.Normal(-10, 30))
std_side = numpyro.sample("std_side", dist.Uniform(5, 30))
numpyro.sample("obs", dist.MixtureSameFamily(dist.Categorical(probs=jnp.array(weights)), dist.Normal(loc=jnp.array([mu_main, mu_side]), scale=jnp.array([std_main, std_side]))), obs=data)
```

Thanks in advance.

Medrano.