janwillem wrote: > Thanks for the explanation. I moved the doit()s to later in the script > where everything boils down to scalars and that works. > > I had tried jacobean but sometimes that gives me AttributeError. Here > is the real thing I was working on: > > def lpu_nonlin(Y, X, C): > """ > LPU for non-linear measurement models according to > Giovanni Mana and Francesca Pennecchi > Metrologia 44 (2007) 246-251 > Inputs: > Y measurement model e.g. Y = (X[0] / X[1]) - X[2] > X vector of input quantities > C covariance matrix > Returns symbolic equations for: > expectation using LPU with 3rd degree Taylor series > variance using LPU with 3rd degree Taylor series > expectation using LPU with 1st degree Taylor series > variance using LPU with 1st degree Taylor series > """ > #Jacobean vector > J = sympy.Matrix([Y]).jacobian(X).T > #Hessian matrix > H = sympy.hessian(Y, X) > > #Components of Mana & Pennacchi equations 2 and 3 > JTC = J.T * C > JTCJ = (J.T * C * J)[0, 0] > HC = H * C > TrHC = HC.trace() > HC2 = HC**2 > TrHC2 = HC2.trace() > dTrHC = sympy.Matrix([sympy.diff(TrHC, X)]).T > #dTrHC = sympy.Matrix([TrHC]).jacobean(X).T > > #Expectation, Mana & Pennacchi eqn 2 > exp3_y = (Y + TrHC / 2).doit() > #variance, Mana & Pennacchi eqn 3 > var3_y = (JTCJ + TrHC2/2 + (JTC * dTrHC)[0, 0]).doit() > > #according to GUM framework LPU > exp1_y = Y > var1_y = JTCJ.doit() > return exp3_y, var3_y, exp1_y, var1_y > > For getting J I already changed from diff to jacobean as you suggested > but doing the same a few lines further on with TrHC fails (the > commented version, TrHC is a scalar function of X and C). > > Why can that be?
After a short look I guess its a simple typo: you wrote 'jacobean' instead of 'jacobian' > And why is the syntax of hessian and jacobean so > different? My guess: the hessian is typically calculated for a scalar function while the jacobian is usually calculated for a vector field (i.e a n-by-1 matrix or its transpose). Therefore 'jacobian' is a method of the Matrix-class while hessian is a function that lives directly in the sympy namespace. As I said: Thats just a guess. Regards. Bastian. -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to sy...@googlegroups.com. To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
