How to convert into double?

Max
2016 2 5
1 个回答
回答已采纳
更新时间：2016 2 6
3 次查看（30 天）

0 个投票

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 个评论
显示 5更早的评论 隐藏 5更早的评论

John D'Errico
John D'Errico 2016-2-5
编辑：John D'Errico 2016-2-5
I assume that you meant to type logL = 0, instead of ogL = 0. That may in fact be why you have trouble. I don't know.
But NO NO NO NO! Then you write:
H_inverse=inv(H_negatv);
But that doesn´t work. So I tried out:
h = 1\H_negatv.
The above does not compute the inverse of a matrix!!!!!!!!!!!!
A = sym(magic(3))
A =
[ 8, 1, 6]
[ 3, 5, 7]
[ 4, 9, 2]
1\A
ans =
[ 8, 1, 6]
[ 3, 5, 7]
[ 4, 9, 2]
So, NO, it is NOT good so far.
Walter Roberson
Walter Roberson 2016-2-5
My tests show that inverting H_negatv is not practical on my system (8 Gb) because inverting it runs out of memory. The length of the hessian is about 4 megabytes.
Switching to using \ instead of inv() is not going to help.
As I am pretty sure I said before at some point, 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.
Max
Max 2016-2-5
Hey John,
but what´s about that:
"Note: MATLAB® computes X^(-1) and inv(X) in the same manner, and both are subject to the same limitations." see, also: http://de.mathworks.com/help/matlab/ref/inv.html
John D'Errico
John D'Errico 2016-2-5
编辑：John D'Errico 2016-2-5
READ WHAT I SAID! Think about what you see in my comment.
As I showed, 1\A does NOT compute a matrix inverse, unless you think that magic(3) is it own inverse.
In fact, when you write the expression 1\A in MATLAB, what you did was to solve the system
1*X = A
The solution to that problem is NOT the inverse of the matrix A. The solution is A.
So giving me lots of documentation that shows that inv(A) is equivalent to x^-1 is irrelevant. None of those sources claim that you will get a matrix inverse from the form 1\A.
PERIOD.
Walter Roberson
Walter Roberson 2016-2-5
You need eye(size(A))\A with square nonsingular A to have the equivalent of inv(A)
Max
Max 2016-2-6
编辑：Walter Roberson 2016-2-6
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 2016-2-6
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.

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Walter Roberson
Walter Roberson 2016-2-5

0 个投票

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 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Walter Roberson
Walter Roberson 2016-2-6
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Number Theory 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by