At least one END missing, yet I have written end
10 views (last 30 days)
Show older comments
Hi, I am trying to estimate the parameters of a distribution using the maximum likelihood estimator, and this is my loglikelihood function/objective function
I have written a code for the formula for the objective function as so (filename sumbgg)
function F = sumbgg(p,x)
sum1 = 0;
sum2 = 0;
sum3 = 0;
sum4 = 0;
sum5 = 0;
sum6 = 0;
sum7 = 0;
sum8 = 0;
sum9 = 0;
sum10 = 0;
for i=1:length(x)
sum1 = sum1 + log(1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)));
sum2 = sum2 + ((1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))^p(1))*log(1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))/log(1-(1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))^p(1));
sum3 = sum3 + exp(p(3)*x(i));
sum4 = sum4 + (exp(p(3)*x(i)-1)*exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))/(1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)));
sum5 = sum5 + ((exp(p(3)*x(i)-1)*exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))*(1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))^(p(1)-1))/(1-(1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))^p(1));
sum6 = sum6 + x(i);
sum7 = sum7 + (exp(p(3)*x(i)))*(1-p(3)*x(i));
sum8 = sum8 + (exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))*(p(3)*x(i)*exp(p(3*x(i)))+exp(p(3)*x(i))+1)/(1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)));
sum9 = sum9 + (((1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))^p(1))*exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))*(p(3)*x(i)*exp(p(3*x(i)))+exp(p(3)*x(i))+1)/(1-(1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1)))^p(1));
sum10 = sum10 + log(1-((1-exp(-(p(2)/p(3))*(exp(p(3)*x(i))-1))))^p(1));
end
F = [(length(x)/p(1))+(p(4)*sum1)-((p(5)-1)*sum2);
(length(x)/p(2))+(length(x)/p(3))-(1/p(3))*sum3+(((p(1)*p(4)-1)/p(3))*sum4)+((p(1)*(p(5)-1)/p(3))*sum5);
(-length(x)*p(2)/(p(3)^2))+ sum6 + p(2)*sum7/(p(3)^2)-p(2)*(p(4)*p(1)-1)*sum8/(p(3)^2)+p(1)*p(2)*(p(5)-1)*sum9/(p(3)^2);
-length(x)*(psi(p(4))+psi(p(4)+p(5)))+p(1)*sum1;
-length(x)*(psi(p(4))+psi(p(4)+p(5)))+p(1)*sum10];
F = double(F);
end
And I want to find the values for each parameter as below (filename calc)
p = [2 0.005 1.16 0.01 0.01];
x = xlsread('Time of First Failure of 500 MW Generators.xlsx');
i = 1;
while(max(sumbgg(p,x))>=10^4 && i<=9999)
option = optimset();
y = fsolve(@sumbgg,p,option,x);
p = y;
fprintf('Iterasi ke-%d\n',i);
i = i + 1;
end
y = sumbgg(x,a);
When I tried to run it returned this error:
Error: File: calc.m Line: 5 Column: 1
At least one END is missing. The statement beginning here does not have a matching end.
I don't understand what I have to do since I have ended the while loop, does anybody know how to fix this?
Also, I have read about optimset but I still don't understand what method it uses to optimize the functions is it perhaps Newton Raphson?
I am new to Matlab so any advice would be much appreciated, thank you!
7 Comments
Accepted Answer
Jan
on 8 Jun 2022
Edited: Jan
on 8 Jun 2022
At least one bug is here:
sum8 = sum8 + (exp(-(p(2) / p(3)) * (exp(p(3) * x(i)) - 1))) * ...
(p(3) * x(i) * exp(p(3 * x(i))) + exp(p(3) * x(i)) + 1) / ...
... % ^^^^^^^^^ This is no valid index for p
... % exp(p(3) * x(i)) is meant
(1 - exp(-(p(2) / p(3)) * (exp(p(3) * x(i)) - 1)));
This appears again:
sum9 = sum9 + (((1 - exp(-(p(2) / p(3)) * ...
(exp(p(3) * x(i)) - 1)))^p(1)) * exp(-(p(2) / p(3)) * ...
(exp(p(3) * x(i)) - 1))) * ...
(p(3) * x(i) * exp(p(3*x(i))) + exp(p(3) * x(i))+1) / ...
... % ^^^^^^^^
(1 - (1 - exp(-(p(2)/p(3)) * (exp(p(3) * x(i)) - 1))) ^ p(1));
I think, this is much easier to debug:
function F = sumbgg(p,x)
sum1 = 0;
sum2 = 0;
sum3 = 0;
sum4 = 0;
sum5 = 0;
sum6 = 0;
sum7 = 0;
sum8 = 0;
sum9 = 0;
sum10 = 0;
n = length(x);
for i = 1:n
e = exp(p(3) * x(i));
d = exp(-p(2) / p(3) * (e - 1));
c = 1 - d;
b = exp(p(3) * x(i) - 1);
a = c^p(1);
sum1 = sum1 + log(c);
sum2 = sum2 + a * log(c) / log(1 - a);
sum3 = sum3 + e;
sum4 = sum4 + b * d / c;
sum5 = sum5 + b * d / c * a / (1 - a);
sum6 = sum6 + x(i);
sum7 = sum7 + e * (1 - p(3)*x(i));
sum8 = sum8 + d * (p(3)*x(i) * e + e + 1) / c;
sum9 = sum9 + (a * d) * (p(3) * x(i) * e + e + 1) / (1 - a);
sum10 = sum10 + log(1 - a);
end
n = length(x);
p3_2 = p(3)^2;
F = [n / p(1) + p(4) * sum1 - (p(5) - 1) * sum2;
n / p(2) + n / p(3) - 1 / p(3) * sum3 + ...
(p(1) * p(4) - 1) / p(3) * sum4 + p(1) * (p(5) - 1) / p(3) * sum5;
-n * p(2) / p3_2 + sum6 + p(2) * sum7 / p3_2 - ...
p(2) * (p(4) * p(1) - 1) * sum8 / p3_2 + ...
p(1) * p(2) * (p(5) - 1) * sum9 / p3_2;
-n * (psi(p(4)) + psi(p(4) + p(5))) + p(1) * sum1;
-n * (psi(p(4)) + psi(p(4) + p(5))) + p(1) * sum10];
F = double(F);
end
Remove clutter from the code: Unnecessary parentheses (if they do not support clarity), repeated expensive calculations (exp and power are very expensive!). USe temporary variables for repeated calculations. Then the inner structure of the formula gets more obvious and this supports the debugging.
3 Comments
Jan
on 9 Jun 2022
Matlab displays the values of arrays in a scaled format. 1.0e+03 mean 1.0 * 10^3, or 1000.
Usually NaN are produced by dividing by zero. Simply let Matlab stop, if this happens and check the local variables:
dbstop if naninf
More Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!