Error on Symbolic calculation : "Empty sym : 0-by-1"

I need the help of people who have excellent talents in calculating using Matlab.

The equations to be implemented through Matlab coding is as follows.

The results to be obtained through this coding is written in "Answer) ~" above.

And, the data related to this problem(or eqution) is below.

I have tried with the following coding but kept getting the same results as the title. I would like to get helpful advice on this part.

I would appreciate it if you could give me a guide on the solution or error cause.

v1 = [393 393 393 393 393 393 393 393 393 393 ; 3850 4340 4760 5320 5740 6160 6580 7140 7980 8960];
v2 = [408 408 408 408 408 408 408 408 408 408 ; 3300 3720 4080 4560 4920 5280 5640 6120 6840 7680];
v3 = [423 423 423 423 423 423 423 423 423 423 ; 2750 3100 3400 3800 4100 4400 4700 5100 5700 6400];

syms B;                     % parameter Beta
syms BB;					% parameter B
syms C;

%------------- equation 1 -------------%
a1_v1 = 0;					
a1_v2 = 0;
a1_v3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(v1(i, :)))
  N = numel(v1(i,j));			
  a1_v1 = N*log(v1(i+1,j)/C/exp(BB/v1(i,j))) + a1_v1;	
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(v2(i, :)))
  N = numel(v2(i,j));			
  a1_v2 = N*log(v2(i+1,j)/C/exp(BB/v2(i,j))) + a1_v2;	
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(v3(i, :)))
  N = numel(v3(i,j));			
  a1_v3 = N*log(v3(i+1,j)/C/exp(BB/v3(i,j))) + a1_v3;	
 end

a2_v1 = 0;
a2_v2 = 0;
a2_v3 = 0;

i = 1;					% "Stress Lv. 1" in the third term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(v1(i, :)))
  N = numel(v1(i,j));			
  a2_v1 = N*(v1(i+1,j)/C/exp(BB/v1(i,j)))^B*log(v1(i+1,j)/C/exp(BB/v1(i,j))) + a2_v1;	
 end

i = 1;					% "Stress Lv. 2" in the third term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(v2(i, :)))
  N = numel(v2(i,j));			
  a2_v2 = N*(v2(i+1,j)/C/exp(BB/v2(i,j)))^B*log(v2(i+1,j)/C/exp(BB/v2(i,j))) + a2_v2;	
 end

i = 1;					% "Stress Lv. 3" in the third term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(v3(i, :)))
  N = numel(v3(i,j));			
  a2_v3 = N*(v3(i+1,j)/C/exp(BB/v3(i,j)))^B*log(v3(i+1,j)/C/exp(BB/v3(i,j))) + a2_v3;
 end

a1 = a1_v1 + a1_v2 + a1_v3;
a2 = a2_v1 + a2_v2 + a2_v3; 
N = numel(v1(1,:)) + numel(v2(1,:))+ numel(v3(1,:));		% the first term of ∂Λ/∂β(equation 1)

one = N/B + a1 - a2;			

%------------- equation 2 -------------%
b1_v1 = 0;					
b1_v2 = 0;
b1_v3 = 0;

i = 1;					% "Stress Lv. 1" in the first term of ∂Λ/∂B(equation 2)
 for j = 1:(numel(v1(i, :)))
  N = numel(v1(i,j));
  b1_v1 = N/v1(i,j) + b1_v1;	
 end

i = 1;					% "Stress Lv. 2" in the first term of ∂Λ/∂B(equation 2)
 for j = 1:(numel(v2(i, :)))
  N = numel(v2(i,j));
  b1_v2 = N/v2(i,j) + b1_v2;	
 end

i = 1;					% "Stress Lv. 3" in the first term of ∂Λ/∂B(equation 2)
 for j = 1:(numel(v3(i, :)))
  N = numel(v3(i,j));
  b1_v3 = N/v3(i,j) + b1_v3;	
 end

b2_v1 = 0;					
b2_v2 = 0;
b2_v3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂B(equation 2)
 for j = 1:(numel(v1(i, :)))
  N = numel(v1(i,j));
  b2_v1 = N/v1(i,j)*(v1(i+1,j)/C/exp(BB/v1(i,j)))^B+b2_v1;	
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂B(equation 2)
 for j = 1:(numel(v2(i, :)))
  N = numel(v2(i,j));
  b2_v2 = N/v2(i,j)*(v2(i+1,j)/C/exp(BB/v2(i,j)))^B+b2_v2;	
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂B(equation 2)
 for j = 1:(numel(v3(i, :)))
  N = numel(v3(i,j));
  b2_v3 = N/v3(i,j)*(v3(i+1,j)/C/exp(BB/v3(i,j)))^B+b2_v3;	
 end

b1 = b1_v1 + b1_v2 + b1_v3;
b2 = b2_v1 + b2_v2 + b2_v3; 

two = B*(- b1 + b2);		  

%------------- equation 3 -------------%
c1 = 0;					
c2 = 0;
c3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂C(equation 3)
 for j = 1:(numel(v1(i, :)))
  N = numel(v1(i,j));			
  c1 = N*(v1(i+1,j)/C/exp(BB/v1(i,j)))^B+c1;
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂C(equation 3)
 for j = 1:(numel(v2(i, :)))
  N = numel(v2(i,j));			
  c2 = N*(v2(i+1,j)/C/exp(BB/v2(i,j)))^B+c2;
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂C(equation 3)
 for j = 1:(numel(v3(i, :)))
  N = numel(v3(i,j));			
  c3 = N*(v3(i+1,j)/C/exp(BB/v3(i,j)))^B+c3;
 end

N = numel(v1(1,:)) + numel(v2(1,:)) + numel(v3(1,:));		% the first term of ∂Λ/∂C(equation 3)
three = B/C*(- N + (c1 + c2 + c3));


eqns = [one == 0, two == 0, three == 0];
vars = [B BB C];
answer = solve(eqns, vars);
vpa(answer.B)
vpa(answer.BB)
vpa(answer.C)

[Torsten]

What is the underlying distribution ? Can you give us its pdf ?

[J]

Thank you for your interest. The distribution is a Weibull distribution.

I wonder if you mean file(*.pdf) or probability density function(pdf) when you say pdf.

If it's a pdf file, I inserted the captured file into this inquiry because it couldn't be extracted due to security issues.

There is no error in the equations, and I would appreciate if you could review whether this coding has a wrong structural arrangement or if there is a more suitable command or function.

[Torsten]

I only know Weilbull distribution with two parameters:

https://en.wikipedia.org/wiki/Weibull_distribution

Could you give a link to the one you use with three parameters or include its probability density function (pdf) as written formula ?

[John D'Errico]

@Torsten

In almost all cases I can think of, when a distribution that is normally bounded at zero is given an extra parameter, the third parameter is a translation parameter, allowing the distribution to be bounded by some other value. That is the case of this three parameter Weibull in the stats TB:

https://www.mathworks.com/help/stats/three-parameter-weibull-distribution.html

In there, you would read:

" If Xhas a two-parameter Weibull distribution, then Y=X+c has a three-parameter Weibull distribution with the added location parameter c."

[Torsten]

@John D'Errico

Do you think equation 1-3 stem from the differentiation of the log-likelihood function of this 3-parameter Weilbull distribution ? The equations look quite special - adopted to some physical reaction-kinetics application. So I think it's better if @J gives the special form to be considered here. Or can you recover /\ from its derivatives ?

Answer 1

[Torsten]

Try this code.

By the way: numel(v1(i,j)), numel(v2(i,j)) and numel(v3(i,j)) always equals 1 because v1(i,j), v2(i,j) and v3(i,j) is one single scalar value.

format long
v1 = [393 393 393 393 393 393 393 393 393 393 ; 3850 4340 4760 5320 5740 6160 6580 7140 7980 8960];
v2 = [408 408 408 408 408 408 408 408 408 408 ; 3300 3720 4080 4560 4920 5280 5640 6120 6840 7680];
v3 = [423 423 423 423 423 423 423 423 423 423 ; 2750 3100 3400 3800 4100 4400 4700 5100 5700 6400];
x0 = [4,1800,59];
sol = fsolve(@(x)fun(x,v1,v2,v3),x0)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

sol = 1×3

1.0e+03 * 0.004291582233058 1.861618685614979 0.058984866373623

function res = fun(x,v1,v2,v3)
B = x(1);
BB = x(2);
C = x(3);
%------------- equation 1 -------------%
a1_v1 = 0;					
a1_v2 = 0;
a1_v3 = 0;
i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂β(equation 1)
for j = 1:(numel(v1(i, :)))
    N = numel(v1(i,j));			
    a1_v1 = N*log(v1(i+1,j)/C/exp(BB/v1(i,j))) + a1_v1;	
end
i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂β(equation 1)
for j = 1:(numel(v2(i, :)))
    N = numel(v2(i,j));			
    a1_v2 = N*log(v2(i+1,j)/C/exp(BB/v2(i,j))) + a1_v2;	
end
i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂β(equation 1)
for j = 1:(numel(v3(i, :)))
    N = numel(v3(i,j));			
    a1_v3 = N*log(v3(i+1,j)/C/exp(BB/v3(i,j))) + a1_v3;	
end
a2_v1 = 0;
a2_v2 = 0;
a2_v3 = 0;
i = 1;					% "Stress Lv. 1" in the third term of ∂Λ/∂β(equation 1)
for j = 1:(numel(v1(i, :)))
    N = numel(v1(i,j));			
    a2_v1 = N*(v1(i+1,j)/C/exp(BB/v1(i,j)))^B*log(v1(i+1,j)/C/exp(BB/v1(i,j))) + a2_v1;	
end
i = 1;					% "Stress Lv. 2" in the third term of ∂Λ/∂β(equation 1)
for j = 1:(numel(v2(i, :)))
    N = numel(v2(i,j));			
    a2_v2 = N*(v2(i+1,j)/C/exp(BB/v2(i,j)))^B*log(v2(i+1,j)/C/exp(BB/v2(i,j))) + a2_v2;	
end
i = 1;					% "Stress Lv. 3" in the third term of ∂Λ/∂β(equation 1)
for j = 1:(numel(v3(i, :)))
    N = numel(v3(i,j));			
    a2_v3 = N*(v3(i+1,j)/C/exp(BB/v3(i,j)))^B*log(v3(i+1,j)/C/exp(BB/v3(i,j))) + a2_v3;
end
a1 = a1_v1 + a1_v2 + a1_v3;
a2 = a2_v1 + a2_v2 + a2_v3; 
N = numel(v1(1,:)) + numel(v2(1,:))+ numel(v3(1,:));		% the first term of ∂Λ/∂β(equation 1)
one = N/B + a1 - a2;			
%------------- equation 2 -------------%
b1_v1 = 0;					
b1_v2 = 0;
b1_v3 = 0;
i = 1;					% "Stress Lv. 1" in the first term of ∂Λ/∂B(equation 2)
for j = 1:(numel(v1(i, :)))
    N = numel(v1(i,j));
    b1_v1 = N/v1(i,j) + b1_v1;	
end
i = 1;					% "Stress Lv. 2" in the first term of ∂Λ/∂B(equation 2)
for j = 1:(numel(v2(i, :)))
    N = numel(v2(i,j));
    b1_v2 = N/v2(i,j) + b1_v2;	
end
i = 1;					% "Stress Lv. 3" in the first term of ∂Λ/∂B(equation 2)
for j = 1:(numel(v3(i, :)))
    N = numel(v3(i,j));
    b1_v3 = N/v3(i,j) + b1_v3;	
end
b2_v1 = 0;					
b2_v2 = 0;
b2_v3 = 0;
i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂B(equation 2)
for j = 1:(numel(v1(i, :)))
    N = numel(v1(i,j));
    b2_v1 = N/v1(i,j)*(v1(i+1,j)/C/exp(BB/v1(i,j)))^B+b2_v1;	
end
i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂B(equation 2)
for j = 1:(numel(v2(i, :)))
    N = numel(v2(i,j));
    b2_v2 = N/v2(i,j)*(v2(i+1,j)/C/exp(BB/v2(i,j)))^B+b2_v2;	
end
i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂B(equation 2)
for j = 1:(numel(v3(i, :)))
    N = numel(v3(i,j));
    b2_v3 = N/v3(i,j)*(v3(i+1,j)/C/exp(BB/v3(i,j)))^B+b2_v3;	
end
b1 = b1_v1 + b1_v2 + b1_v3;
b2 = b2_v1 + b2_v2 + b2_v3; 
two = B*(- b1 + b2);		  
%------------- equation 3 -------------%
c1 = 0;					
c2 = 0;
c3 = 0;
i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂C(equation 3)
for j = 1:(numel(v1(i, :)))
    N = numel(v1(i,j));			
    c1 = N*(v1(i+1,j)/C/exp(BB/v1(i,j)))^B+c1;
end
i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂C(equation 3)
for j = 1:(numel(v2(i, :)))
    N = numel(v2(i,j));			
    c2 = N*(v2(i+1,j)/C/exp(BB/v2(i,j)))^B+c2;
end
i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂C(equation 3)
for j = 1:(numel(v3(i, :)))
    N = numel(v3(i,j));			
    c3 = N*(v3(i+1,j)/C/exp(BB/v3(i,j)))^B+c3;
end
N = numel(v1(1,:)) + numel(v2(1,:)) + numel(v3(1,:));		% the first term of ∂Λ/∂C(equation 3)
three = B/C*(- N + (c1 + c2 + c3));
res = [one,two,three];
end

[J]

You found the answer using the nonlinear solver.

Finally, I wonder if you can provide me one more helpful instruction.

It is judged that the code you guided is possible only when we know the answer to some extent.

That is, when I enter with "x0 = [10,2000,100];", other answers came out due to the limitation on the iterations of calculation.

>> sol = 1.0e+03 *

0.000086087873255 8.006624830398081 0.095124448408611

For this reason, I want to supplement it so that I can get the accurate answer even if any initial value will be put. For example, I want to ease the maximum number of iterations. Using this type of approach, I would like to ask if you can give me a little more guidance on coding.

Thank you.

In terms of this phrase below,

"By the way: numel(v1(i,j)), numel(v2(i,j)) and numel(v3(i,j)) always equals 1 because v1(i,j), v2(i,j) and v3(i,j) is one single scalar value."

You made the right point, but if N_i are typed with the number of datat in the columns, the answer will probably be wrong.

For now, this coding method seems to be the closest way to the answer.(Structurally, it may need to be supplemented a little more.)

@Torsten @John D'Errico

[Torsten]

Your equations are nonlinear in the parameters - there is no such "foolproof" method you are searching for. Convergence to the "correct" parameters will always depend on your initial values.

If you can write down the probability density function from which you derived the maximum likelihood function /\ , you can try MATLAB's "mle" with this custom probability density function option and your data from above. Maybe using "mle" to determine maximum likelihood estimates for the parameters C, B and beta is more stable and has a greater region of attraction than simply applying a nonlinear solver to d/\ / dp = 0 as I did above.

[J]

Thanks to you, I was able to get answers to my questions. I am sorry for the delay in adopting the answer. Thank you for your clear help. @Torsten

Answer 2

[Torsten]

v1 = [393 393 393 393 393 393 393 393 393 393 ; 3850 4340 4760 5320 5740 6160 6580 7140 7980 8960];
v2 = [408 408 408 408 408 408 408 408 408 408 ; 3300 3720 4080 4560 4920 5280 5640 6120 6840 7680];
v3 = [423 423 423 423 423 423 423 423 423 423 ; 2750 3100 3400 3800 4100 4400 4700 5100 5700 6400];
v = [v1,v2,v3];
V = v(1,:);
T = v(2,:);
syms beta B C
f = log(prod(beta/C*exp(-B./V).*(T/C.*exp(-B./V)).^(beta-1).*exp(-(T/C.*exp(-B./V)).^beta)));
dfdbeta = diff(f,beta);
dfdB = diff(f,B);
dfdC = diff(f,C);
df = matlabFunction([dfdbeta,dfdB,dfdC],'Vars',{[beta,B,C]});
format long
sol = fsolve(df,[4,1800,59],optimset('TolFun',1e-12,'TolX',1e-12))

Equation solved, inaccuracy possible.

The vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

sol = 1×3

1.0e+03 * 0.004291582233070 1.861618665737442 0.058984869277753

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

df(sol)

ans = 1×3

1.0e-12 * -0.431531027276132 0.009315108149980 0.056401667616809

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>