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

8 次查看（过去 30 天）
Kashif Naukhez
Kashif Naukhez 2023-10-6
回答： Walter Roberson 2023-12-14
Ran in:
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)

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

回答（3 个）

Torsten
Torsten 2023-10-6
编辑：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 个评论
Torsten
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 ?

请先登录，再进行评论。

Mathieu NOE
Mathieu NOE 2023-10-6
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更早的评论
Kashif Naukhez
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.

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2023-12-14
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Surrogate Optimization 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by