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); endend% 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); endend% 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); endend% 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); endend% 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); endendnum(:,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); endend% 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); endendFhatx = 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; endendfor 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; endend% 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)); endfor i = 1:M Lsum = Lsum + L(i,1);endloglikelihood = (-(M^(-1)) * Lsum);end
Could you please help me? Thank you so much!
Best Answer