MATLAB - Maximum likelihood for periodic function (angles)

Question

In the MATLAB stats tutorials there is a section called "Fitting a More Complicated Distribution: A Mixture of Two Normals" http://www.mathworks.com/help/stats/examples/fitting-custom-univariate-distributions.html

pdf_normmixture = @(x,p,mu1,mu2,sigma1,sigma2) ...
                     p*normpdf(x,mu1,sigma1) + (1-p)*normpdf(x,mu2,sigma2);
lb = [0 -Inf -Inf 0 0];
ub = [1 Inf Inf Inf Inf];
start = [pStart muStart sigmaStart sigmaStart];
paramEsts = mle(x, 'pdf',pdf_normmixture, 'start',start, 'lower',lb, 'upper',ub)

I would like to apply the same methodology for fitting two or more normals to a univariate set of values that I have, but within a periodic domain. That is, angles that have values of 0° to 360° linked together as a circular range. I am not sure how to declare it in order to make MATLAB understand this kind of terminology.

Is it possible to change this implementation to add the circular range case?

Regards, Ignacio

Answer 1

I am using the von Mises distribution http://en.wikipedia.org/wiki/Von_Mises_distribution which is considered a close approximation to the wrapped normal (see the conditions in Topics in ''circular statistics S. Rao Jammalamadaka, A. SenGupta''). Unfortunately, I do not have the matlab to test it, but I think the code is running. Thus, something like this can be done if the approximation is valid:

Main function:

% you should provide the column vector theta 0-2pi

n=size(theta,1);

mu=0;

k=1;

theParameters=[mu;k];
options = optimset('TolFun',0.01);
outputPar = fminsearch('ml',theParameters,options,n,theta);

ML function

function mLike=ml(theParameters,n,theta)

mu=theParameters(1,1);
k=theParameters(2,1);

theSum=0;
for i=1:n
   theSum=theSum+k*cos(theta(i,1)-mu);
end    
mLike=-n*log(2*pi*besselj(0,k)) + theSum;
mLike=-mLike;

I hope it help! There is also a toolbox in R that deals with this kind of estimation http://cran.r-project.org/web/packages/circular/circular.pdf.

If you have problems with the positivity of k, in order to avoid dealing with constrained optimization, do k=exp(kk) and estimate kk instead.