How to convert into double?

13 views (last 30 days)
I wrote the following code:
tfail = [5571.760,5573.742,5654.457,6079.693,6081.927,6172.915,6515.064,6517.515,6617.308,7095.558,7098.298,7209.831,7530.929,7533.885,7654.224,7966.300,7969.472,8098.617,8401.671,8405.059,8543.009,8982.166,8985.843,9135.533,9852.908,9857.017,10024.38,10868.774,10873.387,11061.234];
beta_hat = 4.2915822;
B_hat = 1861.6186657;
C_hat = 58.9848692;
syms t B beta C
y(t) = (exp(-B/((heaviside(t)-heaviside(t-2000))*(330)+(heaviside(t-2000)-heaviside(t-3000))*(350)+...
for i=1:n
I(i) = int(y(t),t,0,tfail(i));
logL =logL+log((beta*y_new(i)*(I(i))^(beta-1))*exp(-((I(i))^beta)));
p = int(y(t),t,0,14000);
u = beta*log(p);
du_dB = diff(u,B);
du_dbeta = diff(u,beta);
du_dC = diff(u,C);
du_dB_sub = subs(du_dB,{beta,B,C},{beta_hat,B_hat,C_hat});
du_dbeta_sub = subs(du_dbeta,{B,C},{B_hat,C_hat});
du_dC_sub = subs(du_dC,{beta,B,C},{beta_hat,B_hat,C_hat});
now I would like to calculate the inverse of H_negatv by using:
But that doesn´t work. So I tried out:
h = 1\H_negatv.
That´s good so far.
But now I would do sth. like that:
F_direct = w;
In according to calculate:
Var_B_hat_direct = double(F_direct(2,2));
But I can´t do that in MATLAB.
Does somebody have an idea how to solve that problem?
Walter Roberson
Walter Roberson on 6 Feb 2016
You are not looping substituting different values for the symbols in H_negatv, so what you should be doing is substituting in the numeric values and double() the result before doing the inverse, so that you are taking the inverse of a numeric 3 x 3 instead of a large symbolic matrix.

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 5 Feb 2016
Working notes for me.
rational = @(V) sym(V, 'r');
if ismember( exist('hessian'), [2, 3, 5, 6, 8]) %is it an executable function?
Hess = @(M,V) hessian(M,V)
Hess = @(M,V) maple('Student[VectorCalculus][Hessian]', M, maple('convert', V, 'list'));
tfail = rational([5571.760, 5573.742, 5654.457, 6079.693, 6081.927, 6172.915, 6515.064, 6517.515, 6617.308, 7095.558, 7098.298, 7209.831, 7530.929, 7533.885, 7654.224, 7966.300, 7969.472, 8098.617, 8401.671, 8405.059, 8543.009, 8982.166, 8985.843, 9135.533, 9852.908, 9857.017, 10024.38, 10868.774, 10873.387, 11061.234]);
n = length(tfail);
beta_hat = rational(4.2915822);
B_hat = rational(1861.6186657);
C_hat = rational(58.9848692);
syms t B beta C
y = (exp(-B/((heaviside(t) - heaviside(t-2000)) * (330) + (heaviside(t-2000) - heaviside(t-3000)) * (350) + (heaviside(t-3000) - heaviside(t-14000)) * (390))))/C;
Z = rational(0);
LogL = Z;
I = rational(zeros(1,n));
y_new = rational(zeros(1,n));
new_term = rational(zeros(1,n));
for i = 1:n
I(i) = simplify( int(y, t, Z, tfail(i)) );
y_new(i) = simplify( subs(y, t, tfail(i)) );
new_term(i) = log((beta * y_new(i) * (I(i))^(beta-1)) * exp(-((I(i))^beta)));
logL = logL + new_term(i);
p = int(y, t, Z, rational(14000));
u = beta * log(p);
du_dB = diff(u, B);
du_dbeta = diff(u, beta);
du_dC = diff(u, C);
du_dB_sub = subs(du_dB, {beta, B, C}, {beta_hat, B_hat, C_hat});
du_dbeta_sub = subs(du_dbeta, {B,C}, {B_hat,C_hat});
du_dC_sub = subs(du_dC, {beta,B,C}, {beta_hat,B_hat,C_hat});
v = [beta; B; C];
H = Hess(logL, v);
%the next will probably fail, running out of memory
Hs = simplify(H);
H_negatv = -1*Hs;
  1 Comment
Walter Roberson
Walter Roberson on 6 Feb 2016
The above is not intended as a solution: it is a recoding for compatibility with the Maple interface. On my system it will run out of memory attempting to simplify H. And if you skip the simplify() step and go on to taking the inv() it will run out of memory trying to take the inverse.
Your expression is simply too large to take a reasonable symbolic inverse of. But since you have specific numeric values to substitute in, you should substitute those specific numeric values into H, take double() of the result (since it will be completely numeric) and inv() that purely numeric interface.
Except, of course, we do not recommend that you take inv() for any reason other than proving that you can.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by