Estimating the parameter of a complex equation using maximum likelihood estimation (MLE) method
6 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (3)
Torsten
el 6 de Oct. de 2023
Editada: Torsten
el 6 de Oct. de 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 comentario
Torsten
el 6 de Oct. de 2023
Editada: Torsten
el 6 de Oct. de 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
el 6 de Oct. de 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 comentarios
Mathieu NOE
el 11 de Dic. de 2023
hello again @Kashif Naukhez
do you mind accepting my answer (if it has fullfiled your expectations ) ? tx
Walter Roberson
el 14 de Dic. de 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 comentarios
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!