%
% TruncNormMMLE Matlab file, v.1.0.2011
%
%  TruncNormMMLE  calculates the modified maximum likelihood estimators of
%  a random variable with truncated Gaussian distribution with truncation
%  points A and B. This function was used in Bouzas et al. [1] where N(t),
%  a doubly stochastic Poisson process (DSPP), is analyzed under the
%  assumption that its mean, E[N(t)], has a truncated Gaussian distribution
%  at any instant of time.

% Syntax: function [Mu,sd,i,aux]=MMLE(A,B,n,M,S,tol)
%
%  Inputs:
%      A         - Lower truncation point of the truncated Gaussian distribution.
%      B         - Upper truncation point of the truncated Gaussian distribution.
%      n         - Size of the data set.
%      M         - Mean of the sample data. 
%      S         - Standard deviation of the sample.
%      tol       - Maximum error to assume if approximate solutions are given.
%  Output:
%      The function returns Mu and sd, that is, the modified estimations of
%      the mean and of the standard deviation of the truncated Gaussian
%      distribution. It also returns the number of iterations until the
%      solution is adopted, i. Finally, it returns aux=1 if an exact
%      solution is obtained or aux=2 if an approximated solution is given
%      at iteration i or when a maximum of iterations is reached.
%
%  *****NOTE*****
%      This function needs an external one, G, defined below.   
%           function g=G(k,par)
%                   syms x;
%                   Mu=par(1);
%                   Sigma=par(2);
%                   G_=inline('x.^k.*exp(-(x-Mu).^2/(2.*Sigma^2))','k','Mu','Sigma','x');
%                   g=eval(int(G_(k,Mu,Sigma,x),par(3),par(4)));
%  *****NOTE*****
%
%  Created by (1)P.R. Bouzas (paula@ugr.es) and (2)N. Ruiz-Fuentes (nfuentes@ujaen.es)
%
%  Department of Statistics and Operations Research
%  (1) Universidad de Granada. Campus Universitario de Cartuja. 18071, Granada. Spain
%  (2) Universidad de Jaen. Campus Universitario de Las Lagunillas. 23071, Jaen. Spain  
%
%  Version 1.0. November 2011.
%
%  To cite this file, this would be in an appropriate format:
%  Bouzas, P.R. and Ruiz-Fuentes, N. (2011). TruncNormMMLE Matlab file. Free access DIGIBUG
%
%  References:
%  [1] Bouzas, P.R., A.M. Aguilera and Valderrama, M.J., Forecasting a class of 
%              doubly stochastic Poisson process, Statistical Papers, 43, 
%              507--523, 2002.
function [Mu,sd,i,aux]=TruncNormMMLE(A,B,n,M,S,tol)
clc
disp('                                                                            ')
disp('                                                                            ')
m = 40; aux = 0;
par=[M,S,A,B];
for i=1:m
    Mu=par(1);sd=par(2);
    F1=G(1,par)/G(0,par)-M;
    F2=G(2,par)/G(0,par)-(G(1,par)/G(0,par))^2+2*sd^2/n+S^2-1/n;
    F=[F1 F2];
    F11=G(2,par)/(sd^2*G(0,par))-Mu*G(1,par)/(sd^2*G(0,par))-G(1,par)^2/(sd^2*G(0,par)^2)+Mu*G(1,par)/(sd^2*G(0,par));
    F12=G(2,par)/(sd^3*G(0,par))-G(2,par)/(sd^3*G(0,par))-G(1,par)*G(2,par)/(sd^3*G(0,par)^2)+2*G(1,par)^2/(sd^3*G(0,par))-Mu^2*G(1,par)/(sd^3*G(0,par));
    F21=G(3,par)/(sd^2*G(0,par))-Mu*sd^2/(sd^2*G(0,par))-G(1,par)*G(2,par)/(sd^2*G(0,par))+Mu*G(2,par)/(sd^2*G(0,par))-2*F11;
    F22=G(4,par)/(sd^3*G(0,par))-2*Mu*G(3,par)/(sd^3*G(0,par))+sd*Mu^2*G(2,par)/(sd^3*G(0,par))-G(2,par)^2/(sd^3*G(0,par)^2)+2*Mu*G(1,par)*G(2,par)/(sd^3*G(0,par)^2)-Mu^2*G(2,par)/(sd^3*G(0,par))-2*F12;
    DF=[F11,F12;F21,F22];    
    [mu1,sd1]=solve('DF(1,1)*(mu1-Mu)+DF(1,2)*(sd1-sd)=F(1,1)','DF(2,1)*(mu1-Mu)+DF(2,2)*(sd1-sd)=F(1,2)','mu1,sd1');
    par(1)=eval(mu1);par(2)=eval(sd1);mu1=par(1);sd1=par(2);
    F1=G(1,par)/G(0,par)-M;
    F2=G(2,par)/G(0,par)-(G(1,par)/G(0,par))^2+2*sd1^2/n+S^2-1/n;
    F=[F1 F2];
    if F==[0 0], 
        aux = 1; 
        break;
    elseif sqrt((mu1-Mu)^2)+sqrt((sd1-sd)^2)<=tol, 
        aux = 2;
        break;
    end
    Mu=mu1;sd=sd1;par=[Mu,sd,A,B];
end
disp('------------------------------'),
    if aux == 1, 
        disp('Exact solution:  '),
        disp([Mu sd]),
    elseif  aux == 2, 
        disp('Approximate solution: '),
        disp([Mu, sd]),
    end
disp('Number of iterations:    '),disp(i),