MATLAB - massima verosimiglianza per la funzione periodica (angoli)

Question

Nel MATLAB Statistiche tutorial v'è una sezione chiamata "Montaggio di una distribuzione più complicato: una miscela di due normali" 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) 

vorrei applicare la stessa metodologia per il montaggio di due o più normali in un insieme univariato di valori che ho, ma all'interno di un dominio periodico. Vale a dire, angoli con valori compresi tra 0 ° e 360 ​​° collegati insieme come un intervallo circolare. Non sono sicuro di come dichiararlo per far capire a MATLAB questo tipo di terminologia.

È possibile modificare questa implementazione per aggiungere il caso di intervallo circolare?

saluti, Ignacio

Answer 1

Sto usando la distribuzione di von Mises http://en.wikipedia.org/wiki/Von_Mises_distribution che è considerato una buona approssimazione alla normale avvolto (vedere le condizioni in Topics in '' statistiche circolari S. Rao Jammalamadaka, A. Sengupta' '). Sfortunatamente, non ho il MATLAB per testarlo, ma penso che il codice sia in esecuzione. Quindi, qualcosa di simile può essere fatto se l'approssimazione è valida:

Funzione principale:

% è necessario fornire il vettore colonna 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); 

funzione ML

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; 

Spero sia d'aiuto! C'è anche una cassetta degli attrezzi in R che si occupa di questo tipo di stima http://cran.r-project.org/web/packages/circular/circular.pdf.

In caso di problemi con la positività di k, per evitare di occuparsi dell'ottimizzazione vincolata, eseguire k = exp (kk) e stima kk.