ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 21 Jul 2021 09:58:41 +0200Assuming matrix is positive semi-definitehttps://ask.sagemath.org/question/58070/assuming-matrix-is-positive-semi-definite/ I'm trying to find moments of of transformations of the multivariate normal distribution.
Let $x,y,z\sim N(\mu, \Sigma)$, I want to find out (for example) $E[e^x], E[e^{2x}], E[e^xy]$.
I have defined the following:
x, y, z = var('x y z')
X = vector([x,y,z])
mx, my, mz = var('mx my mz')
mu = vector([mx, my, mz])
sx, sy, sz = var('sx sy sz')
sxy, sxz, syz = var('sxy sxz syz')
assume(sx > 0)
assume(sy > 0)
assume(sz > 0)
S = matrix(SR, [[sx^2, sxy, sxz], [sxy, sy^2, syz], [sxz, syz, sz^2]])
f = (1 / sqrt(S.det() * (2*pi)^3))*exp(-0.5*((X - mu) * S.inverse() * (X - mu)))
And now I'm trying some integration, for example:
(x*f).integrate(z, -Infinity, Infinity).integrate(y, -Infinity, Infinity).integrate(x, -Infinity, Infinity)
I'd expect this to yield `mx` but of course life ain't easy so I get errors related to the covariance matrix:
ValueError: Computation failed since Maxima requested additional constraints; using the 'assume' command before evaluation *may* help (example of legal syntax is 'assume(sx^2*(sy^2*sz^2-syz^2)
>0)', see `assume?` for more details)
Is sx^2*(sy^2*sz^2-syz^2)
-sxy^2*sz^2
+2*sxy*sxz*syz
-sxz^2*sy^2 positive or negative?
Even for 1D integration (i.e only over $z$, `(x*f).integrate(z, -Infinity, Infinity)`) I get errors:
ValueError: Computation failed since Maxima requested additional constraints; using the 'assume' command before evaluation *may* help (example of legal syntax is 'assume(sx^2*(sy^2*sz^2-syz^2)
>0)', see `assume?` for more details)
Is sx^2*(sy^2*sz^2-syz^2)
-sxy^2*sz^2
+2*sxy*sxz*syz
-sxz^2*sy^2 positive or negative?
I have tried assuming everything he asks but I still get the same errors. Is there any nice way of telling `SAGE` to assume the matrix `S` is positive-semidefinite?
PS I have managed to work with two 2D case, by using correlation and the "full" version of the formula (rather than the vectorized):
x = var('x')
y = var('y')
mux = var('mux')
muy = var('muy')
sx = var('sx')
sy = var('sy')
r = var('r')
assume(r<=1)
assume(r>=-1)
k = 1/(2*pi*sx*sy*sqrt(1-r^2))
fxy = k*exp(-(1/(2*(1-r^2)))*(((x-mux)/sx)^2 - 2*r*((x-mux)/sx)*((y-muy)/sy) + ((y-muy)/sy)^2))
If I run `integrate(integrate(x*fxy, y, -Infinity, Infinity), x, -Infinity, Infinity).simplify()`, the result is eventually the correct one (`mux`).guyaspWed, 21 Jul 2021 09:58:41 +0200https://ask.sagemath.org/question/58070/Help with matrices over multivariable polynomial ringhttps://ask.sagemath.org/question/33310/help-with-matrices-over-multivariable-polynomial-ring/ I want to work with matrices over a multivariable polynomial ring.
I want the matrix
[x0^2,x1^2,x2^2]<br>
[x0^4,x1^4,x2^4]<br>
[x0^8,x1^8,x2^8]
so I can take the determinate of it. I have
R = PolynomialRing(GF(2), 3, 'x')
which is a "Multivariate Polynomial Ring in x0, x1, x2 over Finite Field of size 2". I try
M = MatrixSpace(R,3,3,sparse=True)
which is the "Full MatrixSpace of 3 by 3 sparse matrices over Multivariate Polynomial Ring in x0, x1, x2 over Finite Field of size 2". I am not even sure what "sparse" is.
Then I try
A = M([x0^2,x1^2,x2^2, x0^4,x1^4,x2^4, x0^8,x1^8,x2^8])<br>
or<br>
A = M([[x0^2,x1^2,x2^2], [x0^4,x1^4,x2^4], [x0^8,x1^8,x2^8]])
And it says "name 'x0' is not defined"
I have looked for examples in the Sage documentation, but I just can get Sage to make the matrix above.
Eventually, I want to do arbitrate number of variables and arbitrary n-by-n matrices.
Thank you for your help. MrotsliahTue, 03 May 2016 20:04:31 +0200https://ask.sagemath.org/question/33310/