Plotting a Weibull density function

Hello,

i've created a Maximum-Likelihood function and estimated the parameters as you can see right below:

syms A B b y0 real

%Faktorstufen A%
x11=-1;x21=-1;x31=1;x41=1;
x1=[x11;x21;x31;x41];

%Faktorstufen B%

x12=-1;x22=1;x32=-1;x42=1;
x2=[x12;x22;x32;x42];

%Ausfallzeiten%
%A low, B low%     %A low, B high%    %A high, B low%    % A high, B high%
t1=27;             t2=50;             t3=25;             t4=55;

Versuchsdaten=[t1;t2;t3;t4];

%Log.Likelihood-Funktion%
logL=0;
for i=1:length(Versuchsdaten)
    logL=logL+log(b/(exp(y0+A*x1(i)+B*x2(i))))+log(Versuchsdaten(i)/(exp(y0+A*x1(i)+B*x2(i))))^(b-1)-(Versuchsdaten(i)/exp(y0+A*x1(i)+B*x2(i)))^(b);
end

%differentielle Ableitungen%

dLdy0=diff(logL,y0)
dLdb=diff(logL,b)
%dLdt0=diff(logL,t0)
dLdA=diff(logL,A)
dLdB=diff(logL,B)
%dLdAB=diff(logL,AB)

%range=[0 Inf; 0 Inf; -1 1; -1 1];

[y0_hat, b_hat, A_hat, B_hat] = vpasolve([dLdy0==0, dLdb==0, dLdA==0, dLdB==0],[y0, b, A, B])

The result for the parameters are:

y0_hat =

2.0121168178642018928072071428139

b_hat =

0.38450126565678579374223747169347

A_hat =

0.0045872846670491337649337697413908

B_hat =

0.35116087494702178847354983047048

Then im calculating the parameter T_hat:

if A_hat<1
  x1=-1
else
  x1=1
end

if B_hat<1
      x2=-1  
else
  x2=1
end

T_hat=exp(y0_hat+A_hat*x1+B_hat*x2)

With the result:

T_hat =

5.2402471237931471278213840718739

The problem right now: When trying to plot a Weibull density function with the calculated parameters, i'm getting the following error:

x=0:0.1:130;
f=wblpdf(x,T_hat,b_hat)
plot(f);

Error using symengine
Division by zero.

Error in sym/privBinaryOp (line 973)
            Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

Error in  .^  (line 324)
      B = privBinaryOp(A, p, 'symobj::zip', '_power');

Error in wblpdf (line 50)
y = z.^(B-1) .* w .* B ./ A;

Error in zwei_parametrig (line 66)
f=wblpdf(x,T_hat,b_hat)

I don't really understand where the problem is. I took the 2 solutions for the parameters b_hat and T_hat and tried to plot a density function in an extra Matlab file. It worked.

Can someone tell me, where the problem is/could be?

Sincerly yours, Ronald

Answer 1

[Torsten]

Try

x=0.1:0.1:130;

instead of

x=0:0.1:130;

Best wishes

Torsten.

[Ronald Singer]

Hi Torsten,

it worked out, okay. Thank you.

But isn't there any other solution? The density function should start at 0, not on 0.1.

[Torsten]

If it should start at 0, then in your case (since b_hat < 1), it's +Infinity.

Best wishes

Torsten.

[Ronald Singer]

Now i understand. Thank you very much!.

Sincerly, Ronald