maximum likelihood estimation of triple integral equation using fminsearch

Hi, I would like help with debugging my code which seems not to be working when i try to get some maximum likelihood estimate of a function like this one:L=-0.5log2π-logσ+logγ-log(L^(-γ)-U^(-γ) )+∫_(-∞)^∞▒log{∫_L^U▒〖x^(-γ-1) exp[-0.5((y_i-x)/σ)^2 ] 〗 dx} {1/√2π γ_*/(L_*^(〖-γ〗_* )-U_*^(γ_* ) ) ∫_(L_*)^(U_*)▒〖x^(〖-γ〗_*-1)/ax exp[(y_i-x)/ax]^2 〗 dx}dy .

In the m-file for the function above, where we need to find the maximum likelihood estimate of the parameters L, U, gamma, sigma, when we assume we know the other four parameters L_star, U_star, a, and gamma_star.

I'am not used to putting functions in this forum i hope the function is clear.I have in clear questions in notepad and pdf format that i can send or post on this forum but i dont how to attach the documents to my question.Any suggestions are welcome:my e-mail is:[e-mail].

thanks, Ik

[Walter Roberson]

There are some characters there that did not come through correctly. Possibly they are sigma characters? The first one is immediately after ∫_(-∞)^∞

[Silibelo Kamwi]

Thanks Walter, those three gaps amazed me as well, but they dont represent anything missing in the likelihood function above.I wrote the function using insert math equation in MS word 2010.

Thanks for having a look.

Silibelo

[Walter Roberson]

What do you have so-far ?

Answer 1

[Silibelo Kamwi]

hi walter , here is my code,thanks for taking your time to read it.

function [paramsEst,Funval,exitflag,output] = HHMaxmispecLike2011_3(params0)
%**************************************************************************
% fminsearch = minimize the scalar function HHloglik, starting at initial (params).
% Funval is the value of the function loglikfnct at the solution paramsEst.
% exitflag = describe the exit condition of fminsearch:
% 1 fminsearch converge to a solution paramsEst
% 0 Maximum number of function evaluations or iterations was reached
% -1 Algorithm was terminated by the output function
% output = returns structure output that contains information about the optimization
%**************************************************************************
%z2=10;
%warning off all
%clear all
%params0=[3; 2; 1; 4];
tic
options=optimset('Display','on','MaxFunEvals',15e3,'MaxIter',15e3);
[paramsEst,Funval,exitflag,output] = fminsearch(@HHloglik,params0,options);%,params1);
toc
return
%**************************************************************************
function [loglikhod ] = HHloglik(params)
%,params1)
% loglikfnct(params,data) function returns the negative log-likelihood value
%**************************************************************************
%**************************************************************************
% initial estimates
L = params(1); % Lower limit
%L1=params1(1);
U = params(2); % Upper limit
%U1=params1(2);
a = params(3); % exponent
%a1=params1(3);
%A1=params1(4);%slope parameter for unspecified variance
s=params(4);% variance for the homoscedastic powerlaw
%**************************************************************************
% Terminate if the conditions are not met
if (a < 0||L < 0||U < L||s < 0)%||a1<0 ||L1<0||U1<L1||s<0)%|| beta1 < 0)
loglikhod = 1.0e+20;
return
end
y11=-100;%-inf;
y22=100;%inf;
%n=length(z2);
tol=1.e-8;
D1 =quadl(@integrand2,y11,y22,tol,[],params);
B = log(a) - log(L^(-a) - U^(-a));
A = -(0.5)*log(2*pi)-log(s);
loglike = (A + B)+D1 ;
loglikhod = -loglike; % negative log-likelihood
return
%**************************************************************************
function z4 = integrand2(x2,params)
% integral part of the log-likelihood function.
L= params(1);
%L1=params1(1);
U= params(2);
% U1=params1(2);
a= params(3);
% a1=params1(3);
s=params(4);
% A1=params1(4);
n=length(x2);
z3=zeros(size(x2));
z4=zeros(size(x2));
for k=1:n
y1=L-0.6.*x2(k);
y2=U+0.6.*x2(k);
z5=@integrand4;
y4 =@(z2)(x2(k).^(-a-1)).*(exp(-0.5.*(((z2 - x2(k))./(s)).^2)));%.*z1(k);
z3(k)=quadgk(y4,y1,y2);
z4(k)=(log(z3(k))).*z5(k);
end
return
%****************************************************

%**************************************************************************
function z1 = integrand4(x2)
% integral part of the log-likelihood function.
params1=[3; 20; 1.5; 0.2];
%L1=params1(1);U1=params1(2);a1=params1(3);A1=params1(4);
L1= params1(1);
U1= params1(2);
a1= params1(3);
A1=params1(4);
n=length(x2);
z1=zeros(size(x2));
for k=1:n
%L11=L1-0.6.*x2(k);U11=U1+0.6.*x2(k);
%y11=-inf;
%y22=inf;
%L=3;U=20;a=1.5;A1=0.2;
%s=(A1.*x2(k));
y3 =@(z2)(1/sqrt(2*pi)).*(a1/(L1^(-a1)-U1^(-a1))).* (x2(k).^(-a1-2)/A1).*exp(-0.5.*(((z2 - x2(k))./(A1.*x2(k))).^2));
z1(k)=quadl(y3,L1,U1);
end
return