fminsearch fails to return the optimal parameters to approximate a function my MLE
조회 수: 2 (최근 30일)
이전 댓글 표시
Hello,
I'm trying to solve an optimization problem. I'm approximating an unknown function (a CDF) with a hermite series by maximum likelihood but I'm getting a message saying:
Exiting: Maximum number of function evaluations has been exceeded
- increase MaxFunEvals option.
Current function value: -1.106099
How can I improve fit?
The problem is as follows:
- The MLE to be optimized is:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/147643/image.png)
where y and x are two different vectors of same length T, f is a PDF and F is its CDF (c = 0).
- I need to approximate f with the following Hermite series:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/147645/image.png)
where a, mu and sigma are the unknowns, K is the length of the Hermite series, and phi(.) is a standard normal pdf (with mean mu, and standard deviation sigma) evaluated at x, and truncated at 0 (c = 0).
I would like to do evaluate f(.) 7 different series (i.e. k = 1, ... ,7) and then pick the one with the largest likelihood.
So, in my script I have the following code:
% Initialize a matrix of values for theta and define a matrix for the estimated thetas
K = 7;
initial = ones(K+2,1);
thetahat = zeros(K+2,K);
% Estimate f using Hermite series
for k = 1:K
f = @(theta)LhatSemiParametric2(X,Y,theta,k);
options = optimset('MaxFunEvals',1000);
theta = fminsearch(f,initial,options);
thetahat(:,k) = theta;
end
and Lhatsemiparametric is the following function:
function [ loglikelihood ] = LhatSemiParametric( X,Y,theta,k )
% This file finds the parameter that optimize the approximation of the
% unknown pdf
mu = theta(k+1,1);
sig = theta(k+2,1);
M = length(X);
% construct the numerator of fhatx and fhaty
% 1) construct k polynomials
%polx = zeros(M,k);
poly = zeros(M,k);
for i = 1:M
for l = 1:k
poly(i,l) = ((Y(i,1)- mu) / sig)^(l);
end
end
% 2) compute the numerators for fhatx and fhaty
numy = zeros(M,1);
for i = 1:M
for l = 1:k
numy(i,1) = numy(i,1) + poly(i,l) * theta(l,1);
end
end
% 3) add the truncated normal distribution evaluated at Y
numy(:,1) = ((1 + numy(:,1)).^2) .* ((1/sig) .* ...
(normpdf(((Y(:,1) - mu) ./sig),mu,sig) ./ ...
(1 - normcdf((mu / sig),mu,sig))));
% construct the denominator for fhaty
% 1) define an integration grid for the numerical integration
grid = 0:0.01:1;
grid = transpose(grid);
N = length(grid);
% 2) numerical integration
denominatory = ones(M,1);
for i = 1:M
j = 2;
contribution = zeros(N,1);
while grid(j,1) <= Y(i,1) && j < N
sum = 1; % 1 + polynomials
for l = 1:k
sum = sum + ((((0.5 * (grid(j-1,1) + (grid(j,1))))...
- mu) / sig) ^ (l)) * theta(l,1);
end
contribution(j,1) = (sum ^ 2) * ((1/sig) * ...
(normpdf((((0.5 * (grid(j-1,1) + (grid(j,1)))) - mu) ...
/sig),mu,sig) / (1 - normcdf((mu / sig),mu,sig)))) * ...
(grid(j,1) - grid(j-1,1));
j = j + 1;
end
for j = 1:N
denominatory(i,1) = denominatory(i,1) + contribution(j,1);
end
end
% 3) find fhaty
% fhatx(:,1) = numx(:,1) ./ denominatorx(:,1);
% fhaty = zeros(M,1);
fhaty(:,1) = numy(:,1) ./ denominatory(:,1);
% find Fhatx and Fhaty (CDFs) by numerically integrate fhatx and fhaty
% 1) recompute the general pdf based on the values of a grid
% 1a) construct k polynomials
pol = zeros(N,k);
for i = 2:N
for l = 1:k
pol(i,l) = ((((grid(i,1) + grid(i-1,1))/2) - mu) ./ sig)^(l);
end
end
% 1b) compute the numerators for the general pdf
num = zeros(N,k);
for i = 1:N
for l = 1:k
num(i,1) = num(i,1) + pol(i,l) * theta(l,1);
end
end
num(:,1) = ((1 + num(:,1)).^2) .* ((1/sig) .* ...
(normpdf(((((grid(i,1) + grid(i-1,1))/2) - mu) ./sig),mu,sig) ./ ...
(1 - normcdf((mu / sig),mu,sig))));
% 1c) construct the denominator for the pdf by numerical integration
denominator = ones(N,1);
for i = 3:N
j = 3;
contribution = zeros(N,1);
while grid(j,1) <= grid(i,1) && j < N
sum = 1; % 1 + polynomials
for l = 1:k
sum = sum + (((0.5 * (grid(j-1,1) + 0.5 * (grid(j,1) + ...
(grid(j-2,1)))) - mu) / sig) ^ l) * theta(l,1);
end
contribution(j,1) = (sum ^ 2) * ((1/sig) * ...
(normpdf(((0.5 * (grid(j-1,1) + 0.5 * (grid(j,1) + ...
(grid(j-2,1)))) - mu) ...
/sig),mu,sig) / (1 - normcdf((mu / sig),mu,sig)))) * ...
(grid(j,1) - grid(j-2,1) * 0.5);
j = j + 1;
end
for j = 1:N
denominator(i,1) = denominator(i,1) + contribution(j,1);
end
end
% 1d) find the generic pdf
fhat = zeros(N,1);
for i =1:N
if denominator(i,1) == 0
fhat(i,1) = 0;
disp('error 0 denominator in fhat');
else
fhat(i,1) = num(i,1) ./ denominator(i,1);
end
end
Fhatx = zeros(M,1);
Fhaty = zeros(M,1);
for i = 1:M
j = 2;
while grid(j,1) <= X(i,1)
Fhatx(i,1) = Fhatx(i,1) + (fhat(j,1) * (grid(j,1) - grid(j-1,1)));
j = j + 1;
end
end
for i = 1:M
j = 2;
while grid(j,1) <= Y(i,1)
Fhaty(i,1) = Fhaty(i,1) + (fhat(j,1) * (grid(j,1) - grid(j-1,1)));
j = j + 1;
end
end
% solve the MLE problem
L = zeros(M,1);
Lsum = 0;
loglikelihood = 0;
for i = 1:M
L(i,1) = log((2 * (1 - Fhaty(i,1)) * fhaty(i,1)) / ((1 - Fhatx(i,1))^2));
end
for i = 1:M
Lsum = Lsum + L(i,1);
end
loglikelihood = (-(M^(-1)) * Lsum);
end
Could you please help me? Thank you so much!
댓글 수: 0
채택된 답변
Sruthi Ayloo
2015년 8월 5일
Hi Mic,
From the set-options doc, the default value of “MaxFunEvals” for “fminserach” is 200*length(initialValue).
You could try increasing the value to the default value (or higher) and see if your programs runs successfully.
You could also refer to the following link for more information on debugging when the solver fails:
Hope this helps.
댓글 수: 0
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Systems of Nonlinear Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!