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

Question

Kashif Naukhez 2023-10-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2030009-estimating-the-parameter-of-a-complex-equation-using-maximum-likelihood-estimation-mle-method

回答： Walter Roberson 2023-12-14

MLE.xlsx

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)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-10-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2030009-estimating-the-parameter-of-a-complex-equation-using-maximum-likelihood-estimation-mle-method#answer_1327824

编辑：Torsten 2023-10-6

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2023-10-6

编辑：Torsten 2023-10-6

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 2023-10-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2030009-estimating-the-parameter-of-a-complex-equation-using-maximum-likelihood-estimation-mle-method#answer_1327869

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mathieu NOE 2023-12-11

hello again @Kashif Naukhez

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

Kashif Naukhez 2023-12-14

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

请先登录，再进行评论。

Answer 3

Walter Roberson 2023-12-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2030009-estimating-the-parameter-of-a-complex-equation-using-maximum-likelihood-estimation-mle-method#answer_1371224

在 MATLAB Online 中打开

MLE.xlsx

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