fminsearch fails to return the optimal parameters to approximate a function my MLE

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:

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:

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!

Answer 1

[Sruthi Ayloo]

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:

http://www.mathworks.com/help/optim/ug/when-the-solver-fails.html

Hope this helps.