slight edit the last line of the code shouldve been

y = sumbgg(p,x);

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!

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.

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

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

Start Hunting!