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

Question

0 개 추천

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