# 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).

edit retag close merge delete
{"ft_lang_label":"__label__en","ft_lang_prob":0.6783799,"math_prob":0.9975668,"size":2313,"snap":"2022-05-2022-21","text_gpt3_token_len":790,"char_repetition_ratio":0.12949328,"word_repetition_ratio":0.19277108,"special_character_ratio":0.34716818,"punctuation_ratio":0.14203455,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9999598,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-05-27T22:07:26Z\",\"WARC-Record-ID\":\"<urn:uuid:a44778e8-61ef-47af-89b9-0e7a1d7f3214>\",\"Content-Length\":\"50423\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:93079604-bd22-4368-96c5-6542b827f0e3>\",\"WARC-Concurrent-To\":\"<urn:uuid:b34bd16d-23c1-454d-bb25-ff3ebbf0f00e>\",\"WARC-IP-Address\":\"194.254.163.53\",\"WARC-Target-URI\":\"https://ask.sagemath.org/question/58070/assuming-matrix-is-positive-semi-definite/\",\"WARC-Payload-Digest\":\"sha1:S6TWPAIV3UNIRUKRUW6SIQN5FQ2IL2E3\",\"WARC-Block-Digest\":\"sha1:GH6O67KMNOBGYIXZHERLCEQJPZQGZ775\",\"WARC-Identified-Payload-Type\":\"application/xhtml+xml\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-21/CC-MAIN-2022-21_segments_1652663006341.98_warc_CC-MAIN-20220527205437-20220527235437-00459.warc.gz\"}"}