# Switching between class distributions

**URL:** <https://discourse.pymc.io/t/switching-between-class-distributions/4340>\
**Category:** Questions\
**Created:** [January 10, 2020, 1:39pm UTC](https://discourse.pymc.io/t/switching-between-class-distributions/4340 "2020-01-10T13:39:09Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jnispen](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jnispen/32/2415_2.png) [@jnispen](https://discourse.pymc.io/u/jnispen)\
**Post date:** [January 10, 2020, 1:39pm UTC](https://discourse.pymc.io/t/switching-between-class-distributions/4340/1 "2020-01-10T13:39:09Z")

</div>

Hello,

I am new to probabilistic modeling/programming so please forgive me if I am using the wrong (or weird) terms.  
For my research I am trying to model IR spectra by considering them as a summation of peaks, a baseline model and residual noise. The peak shapes and baseline model will evolve into a more complex model eventually, but I want to start simple and build up complexity as I go along. For now I have Gaussian peaks, a constant (y)-offset and normally distributed noise. Depending on the class to which the spectrum belongs the class y-offset is either from distribution A or B. Currently I am using generated data which looks like this:

 ![data_plot](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/8/89e0328f7859edc95148d3eb695df122af89f806.png)  
The blue lines are from class A and red lines from class B.

The model looks something like:

```
npeaks = 3
nclasses = 2
with pm.Model() as model:
    # priors
    amp = pm.Uniform('amp', 0, x_s.max(), shape=(nclasses,npeaks))
    mu = pm.Normal('mu', mu=np.linspace(x_val.min(), x_val.max(), npeaks), sd=50, 
               shape=(nclasses,npeaks), transform=pm.distributions.transforms.ordered)
    sigma = pm.HalfNormal('sigma', sd=100, shape=(nclasses,npeaks))
    a_ = [pm.Normal('a_%d' % i, mu=0, sd=10, shape=(len(x_s),1)) for i in range(nclasses)]
    
    # spectral lines (per class)
    y_ = [pm.Deterministic('y_%d' % i, (amp[i] * np.exp(-(x_val - mu[i])**2/(2 * sigma[i]**2))).sum(axis=1) 
                       + a_[i]) for i in range(nclasses)]

    # labels (observed from the data, 0 or 1)
    theta = pm.Beta('theta', 1., 1.)
    label = pm.Bernoulli('label', p=theta, observed=Y)

    # switch to spectrum class A/B depending on label
    y_obs = pm.math.switch(pm.math.eq(label, 1), y_[0], y_[1])
    
    sigma_e = pm.Gamma('sigma_e', alpha=1., beta=1.)
    epsilon = pm.HalfNormal('epsilon', sd=sigma_e)

    y_pred = pm.Normal('y_pred', mu=y_obs, sd=epsilon, observed=x_s)

```

So, depending on the label in the training data (Y), the model switches to variable y\_0 or y\_1.

The problem I have with this basic model is this:

- I am now switching between 2 classes but also would like to switch between 4 classes, or N-classes for that matter. Normally I would use a switch case statement, but that is not available. **How can I achieve this in the current context?**

Any help is greatly appreciated!

---

<div class="post-metadata">

**Author:** ![Gon\_F](https://avatars.discourse-cdn.com/v4/letter/g/f14d63/32.png) [@Gon\_F](https://discourse.pymc.io/u/Gon_F)\
**Post date:** [January 11, 2020, 4:58am UTC](https://discourse.pymc.io/t/switching-between-class-distributions/4340/2 "2020-01-11T04:58:26Z")

</div>

Hello. This is definitely possible, and I think you need to utilize a Dirichlet prior to make this work. Look through these two resources for a primer, [link1](https://discourse.pymc.io/t/multivariate-switch-function/2102) and [link2](https://stats.stackexchange.com/questions/234385/how-to-build-a-pymc-model-to-detect-multiple-switch-points).

---

<div class="post-metadata">

**Author:** ![jnispen](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jnispen/32/2415_2.png) [@jnispen](https://discourse.pymc.io/u/jnispen)\
**Post date:** [January 12, 2020, 10:25am UTC](https://discourse.pymc.io/t/switching-between-class-distributions/4340/3 "2020-01-12T10:25:45Z")

</div>

I have looked at the two links, but I think the problems solved are different from the one I am trying to solve.

1. The number of classes is known and fixed (2 or 4 in my case)
2. Per spectrum I know the label Y (to which class the spectrum belongs)

In the end, the goal is to introduce shared variables for Y and the spectral observations (x\_s) so that (after inference) I can:

1. set the variable Y=z and sample from the distribution z
2. set a (unseen) spectrum x\_s[a] and sample from the posterior distribution label, so I can calculate to which class the spectrum belongs.

An ideal solution would be something like this:

```
# labels (observed from the data, 0 or 1)
theta = pm.Beta('theta', 1., 1.)
label = pm.Bernoulli('label', p=theta, observed=Y)

sigma_e = pm.Gamma('sigma_e', alpha=1., beta=1.)
epsilon = pm.HalfNormal('epsilon', sd=sigma_e)

y_pred = pm.Normal('y_pred', mu=y_[label], sd=epsilon, observed=x_s) 

```

In this way, whatever the label code is ([0|1] in this case, but [0|1|2|3] for 4 classes), but that gives the error:

```
`TypeError: list indices must be integers or slices, not ObservedRV`
```