Estimating the parameter of a complex equation using maximum likelihood estimation (MLE) method
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (3)
Torsten
am 6 Okt. 2023
Bearbeitet: Torsten
am 6 Okt. 2023
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
1 Kommentar
Torsten
am 6 Okt. 2023
Bearbeitet: Torsten
am 6 Okt. 2023
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
am 6 Okt. 2023
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1503884/image.png)
with y in log scale it's even better to see how the model matches the data
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1503889/image.png)
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 Kommentare
Mathieu NOE
am 11 Dez. 2023
hello again @Kashif Naukhez
do you mind accepting my answer (if it has fullfiled your expectations ) ? tx
Walter Roberson
am 14 Dez. 2023
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)
0 Kommentare
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!