How to convert into double?
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Max
el 5 de Feb. de 2016
Comentada: Walter Roberson
el 6 de Feb. de 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 comentarios
Walter Roberson
el 6 de Feb. de 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.
Respuesta aceptada
Walter Roberson
el 5 de Feb. de 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 comentario
Walter Roberson
el 6 de Feb. de 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.
Más respuestas (0)
Ver también
Categorías
Más información sobre Linear Algebra en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!