How to convert into double?

3 Ansichten (letzte 30 Tage)
Max
Max am 5 Feb. 2016
Kommentiert: Walter Roberson am 6 Feb. 2016
Hello,
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];
n=length(tfail);
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)+...
(heaviside(t-3000)-heaviside(t-14000))*(390))))/C;
ogL=0;
for i=1:n
tfail(i);
I(i) = int(y(t),t,0,tfail(i));
y_new(i)=subs(y,t,tfail(i));
logL =logL+log((beta*y_new(i)*(I(i))^(beta-1))*exp(-((I(i))^beta)));
end
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});
v=[beta;B;C];
H=hessian(logL,v);
H_negatv=-1*H;
now I would like to calculate the inverse of H_negatv by using:
H_inverse=inv(H_negatv);
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:
w=subs(h,[beta,B,C],[beta_hat,B_hat,C_hat]);
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?
  7 Kommentare
Max
Max am 6 Feb. 2016
Bearbeitet: Walter Roberson am 6 Feb. 2016
Dear John,
first of all, you´re right concerning the matlab command: 1\A. That doesn´t calculate the inverse matrix.
But instead of doing Ainv = inv(A), I also can do: A_inv = A^(-1):
A = [ 1 2 3; 2 3 4; 1 2 5]
A =
1 2 3
2 3 4
1 2 5
inv(A)
ans =
-3.5000 2.0000 0.5000
3.0000 -1.0000 -1.0000
-0.5000 0 0.5000
A_inv
ans =
-3.5000 2.0000 0.5000
3.0000 -1.0000 -1.0000
-0.5000 0 0.5000
My problem is that I can´t do it in MATLAB because of the "out of memory"!!!
Walter Roberson
Walter Roberson am 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.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Walter Roberson
Walter Roberson am 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)
else
Hess = @(M,V) maple('Student[VectorCalculus][Hessian]', M, maple('convert', V, 'list'));
end
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);
end
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 Kommentar
Walter Roberson
Walter Roberson am 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.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by