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:
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_, 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_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!