aecker/moksm — pull request #3 by kqshan (2014-02-26 22:34:15+00:00)
Title: Bugfix/tdist

Unlike the Gaussian case, the MLE for a t-distribution's location and scale parameters are not the sample mean and covariance. This pull request implements a more correct algorithm for fitting the t-distribution parameters. It also contains a test harness to demonstrate the effect of the changes.

In this single-component, univariate, stationary case, we can see that the current MoKsm code overestimates the scale parameter Σ:
![image](https://f.cloud.github.com/assets/6790411/2276219/ea2316ae-9f30-11e3-88f6-fc2beaa5e609.png)
The ground truth (black line) and the parameters given by the stats toolbox's fitdist() function (green line) are given for comparison.

The code in this pull request correctly estimates the parameters:
![image](https://f.cloud.github.com/assets/6790411/2276272/68a29d56-9f31-11e3-8957-04d8815172a6.png)

Overestimating the scale parameter can have unexpected effects in a clustering application. With the current MoKsm code:
![image](https://f.cloud.github.com/assets/6790411/2276353/597ba5f6-9f32-11e3-8320-3d0014e6df25.png)
In this case, about 90% of the spikes have been assigned to cluster 3. The right subplot uses the MoKsm-fitted model to generate the same number of data points; it bears little resemblance to the original training data.

The updated code is more like what we expect:
![image](https://f.cloud.github.com/assets/6790411/2276446/a1d45a68-9f33-11e3-9e3f-23066373a5e5.png)

These examples used very non-Gaussian data (DoF=3) to emphasize the effect the Gaussian assumption in the M-step. In our experience, real data looks more Gaussian than this (DoF ≈ 9), so the effect is not really noticeable.

-----

Description of the algorithm:

One way of attaining the MLE is via EM with the introduction of a new latent variable. This works well for us because we are doing EM anyways for the mixture model. Furthermore, if the degrees-of-freedom parameter is known, then the M-step has closed-form solutions. The idea is that you introduce a gamma-distributed random variable U:

U ~ Gamma(df/2, df/2)

Now if you generate a random variable Y with a normal distribution conditional on U=u:

Y | U=u ~ Normal( mu, C / u )

Then the marginal distribution of Y will be t-distributed with location parameter mu, scale parameter C, and degrees-of-freedom df. There is one u for each observation y. In the EM framework, these u are our latent variables (along with the mixture component indicator variables z). Once the u are "known", then the M-step becomes an MLE over a normal distribution.

McLachlan and Peel (2000), Finite Mixture Models, chapter 7 "Multivariate t Mixtures" gives the following E-step expression for u with respect to component i (eq 7.22, 7.6):

u_i = (df + d) ./ (df + (Ytrain-mu_i)' * inv(C_i) * (Ytrain-mu_i) )

where d is the number of dimensions and Ytrain are the [d x #samples] observations (bsxfun omitted for clarity). Note that u_i → ones(1, #samples) as df → Inf (the Gaussian case). For the M-step, the book gives (in the stationary case) the following closed-form solutions for mu and C for component i (eq 7.32, 7.33):

mu_i = sum(post_i .* u_i .* Ytrain), 2) / sum(post_i .* u_i)
C_i = (post_i .* u_i .* (Ytrain-mu_i) * (Ytrain-mu_i)') / sum(post_i)

where post_i is the posterior expectation that z == i. In our case (where the component means are allowed to drift over time), we first optimize mu_i with C_i fixed (using the RTS smoother), then update C_i with mu_i fixed. This pull request applies these modifications to the MoKsm code.