# Automatically fitting a mixture distribution with two univariate scaled noncentral student's t components

So far my own trials to fit such a mixture distribution to simulated or real data in R were unsuccessful (even if the data was simulated from a two-component t mixture!!!). I'm about to try the same thing in matlab. But several professionals I've been talking to said that it might be in general hard to fit such a distribution because of the heavy tails in the students t and the resulting flatness of the likelihood function.

So I'd be very interested if ANY of you has EVER been able to fit such a distribution successfully - before I try to do the same thing in matlab, stata, C, ...

ALL parameters should be estimated (mixing law, component means, component standard deviations (or alternatively scaling factors), component degrees of freedom). Any algorithm and programming language goes!

For an example of the component's parametrization see e.g. here: http://en.wikipedia.org/wiki/Student%27s_t-distribution#Non-standardized_Student.27s_t-distribution

I agree this can be hard. It may require lots of data and good starting values. Here's an example of how you could try it in MATLAB.

% Generate some data
rng default
n = randn(1000,1); % same as t with mu=0, sigma=1, df=Inf
t = 10 + 3*trnd(5,2000,1);
x = [n;t];

% Define pdf f. To avoid problems with the optimizer trying values out of
% range, work on the logit scale for p and the log scale for the sigma
% and d.f. parameters.
invlogit = @(x) 1./(1+exp(-x));
f = @(x,logitp,mu1,logs1,logdf1,mu2,logs2,logdf2) ...
invlogit(logitp)    *tpdf((x-mu1)/exp(logs1),exp(logdf1))/exp(logs1) + ...
(1-invlogit(logitp))*tpdf((x-mu2)/exp(logs2),exp(logdf2))/exp(logs2);

% Do the fit, allowing lots of iterations and function evaluations
opt = statset('maxiter',1000,'maxfunevals',1000);
p = mle(x,'pdf',f,'options',opt,'start',[-1   0 2 2   5 2 2])

% Convert back to the natural parameters
prob = invlogit(p(1))
mu1 = p(2)
s1 = exp(p(3))
df1 = exp(p(4))
mu2 = p(5)
s2 = exp(p(6))
df2 = exp(p(7))

• Thanks a lot, Tom! Im pretty busy right now and also have a idea how it could work in R. But I`ll definitely have a look at your proposal asap! – Joz May 12 '14 at 15:17