Estimating the parameter of a complex equation using maximum likelihood estimation (MLE) method

I have the equation of the form:

y = A * (1 - (1 - q1) * beta1 * x).^ (1 / (1 - q1)) + (1 - A) * (1 - lambda / beta2 + (lambda / beta2) * exp((q2 - 1) * beta2 * x)).^ (1 / (1 - q2))

I need to determine q1, q2, A, lambda, beta1 and beta2 using the MLE method with known y and x. (attached in excel file)

Can anyone help me in this?

A = readmatrix("MLE.xlsx");
[Asortedx,idx] = sort(A(:,1));
Asortedy = A(idx,2);
plot(Asortedx,Asortedy)

Answer 1

[Torsten]

And this is a probability distribution for arbitrary values of q1, q2, A, lambda, beta1 and beta2 ? I can't believe.

If you were concerned with distribution fitting (for which mle would be the suitable tool), your Excel file would only contain one instead of two columns.

For the difference between curve and distribution fitting and the available MATLAB tools to use, you should have a look at

https://uk.mathworks.com/help/stats/curve-fitting-and-distribution-fitting.html

[Torsten]

You don't have a distribution, you have a function y(x) for which you want to fit parameters such that sum((y-y_measured)^2) is minimal.

Otherwise, you only had one column of data and a probability distribution function that integrates to 1.

Didn't you read the MATLAB page I linked to ?

Answer 2

[Mathieu NOE]

as i am lazzy, I tried first a simpler model with two decaying exponentials (one was insufficient)

y = a.*exp(-b.*x) + c.*exp(-d.*x)

now from my 4 identified parameters (a_sol, b_sol, c_sol, d_sol) you may be able to compute your own params

a_sol = 0.8631

b_sol = 29.5154

c_sol = 0.1946

d_sol = 10.1859

with y in log scale it's even better to see how the model matches the data

data = readmatrix("MLE.xlsx");
x =  data(:,1);
y =  data(:,2);

% remove duplicates 
[x,k,~] = unique(x);
y = y(k);

% resample and interpolate the data
xn = linspace(min(x),max(x),100);
yn = interp1(x,y,xn);

%% 4 parameters fminsearch optimization
f = @(a,b,c,d,x) a.*exp(-b.*x) + c.*exp(-d.*x);
obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),xn)-yn);
D_0 = [1, 20, 0.25, 10]; %initial guess

sol = fminsearch(obj_fun, D_0);
a_sol = sol(1);
b_sol = sol(2);
c_sol = sol(3);
d_sol = sol(4);


yfit = f(a_sol, b_sol, c_sol, d_sol, xn);

R2 = my_R2_coeff(yn,yfit); % correlation coefficient


figure(1);
plot(x, y, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
plot(xn, yfit, '-');
hold off
title(['Exp Fit / R² = ' num2str(R2) ], 'FontSize', 15)
eqn = " y = "+a_sol+ " * exp(-" +b_sol+" *x) + " +c_sol+ " * exp(-" +d_sol+" *x) ";
legend("data",eqn)


figure(2);
semilogy(x, y, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
semilogy(xn, yfit, '-');
hold off
title(['Exp Fit / R² = ' num2str(R2) ], 'FontSize', 15)
eqn = " y = "+a_sol+ " * exp(-" +b_sol+" *x) + " +c_sol+ " * exp(-" +d_sol+" *x) ";
legend("data",eqn)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function R2 = my_R2_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation (R squared) is
    R2 = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

[Kashif Naukhez]

Thanks @Mathieu NOE..But the only problem I am facing is that my answer changes everytime i enter a different initial guess..So i want to obtain the parameters without making any initial guess..Is this possible ?Thanks in advance.

[Torsten]

You use 6 parameters to fit a simple exponential decay that could be described by a function of 2 parameters. In mathematics, this is called "overfitting". The consequence is that you can compensate changes in one parameter by adapting another parameter and still getting the same quality of fit. So whatever initial guess you give to the parameters, there will be a set of parameters nearby that fits your function well. And I think this is the main reason that the gained parameters depend so strongly on the initial values you give to them.

[Mathieu NOE]

hello again @Kashif Naukhez

do you mind accepting my answer (if it has fullfiled your expectations ) ? tx

[Kashif Naukhez]

Mathieu NOE, actually it works very well for the function you have assumed but for my function its not converging.

Answer 3

[Walter Roberson]

The fval (residue) varies a lot, and the parameters vary a lot.

The ga() call really should add in lower bounds and upper bounds, and really should increase the population size

A = readmatrix("MLE.xlsx");
[x, idx] = sort(A(:,1));
y = A(idx,2);
%y = A * (1 - (1 - q1) * beta1 * x).^ (1 / (1 - q1)) + (1 - A) * (1 - lambda / beta2 + (lambda / beta2) * exp((q2 - 1) * beta2 * x)).^ (1 / (1 - q2))
%        1    2      1            0   1    2      10

%P(1) - A
%P(2) - q1
%P(3) - beta1
%P(4) - lambda
%P(5) - beta2
%P(6) - q2
eqn = @(P) sum((P(1) .* (1 - (1-P(2)) .* P(3) .* x).^(1./(1-P(2))) + (1-P(1)) .* (1 - P(4)./P(5) + (P(4)./P(5)) .* exp((P(6) - 1) .* P(5) .* x)) .^ (1 / (1 - P(6))) - y).^2);

numparms = 6;

[bestP, fval] = ga(eqn, numparms)

ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

bestP = 1×6

0.0822 1.3620 6.6886 4.9671 27.0463 26.9959

fval = 0.0405