#we prove the closed form solutions for the unbounded integrals in
# (73) and (76)

restart;
roll := rand(1 .. 20):
with(LinearAlgebra):
#with(Optimization):
#with(plots):
#with(CurveFitting):

#we start with the univariate case
assume(c>0);

#(71)
pf := sqrt(2)*exp(-(x - m)^2/(2*c))/(2*sqrt(Pi*c))  ;

#(72)
rf := q[0] + r[1]*(x-m) 
+ (1/(2!))*r[2]*(x-m)^2
+ (1/(3!))*r[3]*(x-m)^3
+ (1/(4!))*r[4]*(x-m)^4;

#proof of (73)
int( pf*rf , x=-infinity..infinity );

###
#now we switch to the multivariate case
restart;
with(LinearAlgebra):
roll := rand(1 .. 20):


# we write down the proof for a 5D space, any other dimensionality works similarly 
n := 5;

#coordinates of a general point in the space
X := Vector(n,symbol=x);

#Mean and Covariance in (74)
M := Vector(n,symbol=m);
C := Matrix(n,n,symbol=c,shape=symmetric );

#we set M and C to integer random numbers to accelerate computation. Computation with pure #symbolic M and C are possible but make the sheet extremely slow. For setting an random C, we #have to make sure that C is positive definite 
for i from 1 by 1 to n do 
M[i] := (roll()-10)/10:
end do:
HC := Matrix(n,n,shape=symmetric,symbol=hc):
for i from 1 by 1 to n do 
for j from 1 by 1 to n do 
HC[i,j] := (roll()-10)/10:
end do:
end do:
C := Multiply( Transpose(HC), HC):

M;
C;


#compute (74)
CI := MatrixInverse(C):
ff := -(1/2)*simplify(Multiply(Multiply(Transpose(X-M),CI),X-M)):
pf := (1/sqrt((2*Pi)^n*Determinant(C)))*exp( ff   );


#compute (75)
rf := 0:
for i1 from 0 by 1 to 4 do 
for i2 from 0 by 1 to 4-i1 do 
for i3 from 0 by 1 to 4-i1-i2 do 
for i4 from 0 by 1 to 4-i1-i2-i3 do 
for i5 from 0 by 1 to 4-i1-i2-i3-i4 do 
rf := rf + (1/i1!)*(X[1]-M[1])^i1 *(1/i2!)*(X[2]-M[2])^i2 *(1/i3!)*(X[3]-M[3])^i3
*(1/i4!)*(X[4]-M[4])^i4 *(1/i5!)*(X[5]-M[5])^i5 * r[i1,i2,i3,i4,i5]:
end do:
end do:
end do:
end do:
end do:
rf;

#compute the bounded integral (left-hand side of (76) )
h[0] := pf*rf:
for i from 1 by 1 to n do 
h[i] := int( h[i-1], X[i]=-infinity..infinity ):
end do:
int_p_r := h[n];

#compute the right-hand side of (76) 
iint_p_r := r[0, 0, 0, 0, 0]:
#
for i from 1 by 1 to n do 
for j from 1 by 1 to n do 
for k from 1 by 1 to n do 
o[k] := 0:
end do:
o[i] := o[i]+1:
o[j] := o[j]+1:
iint_p_r := iint_p_r + (1/2)*C[i,j]*r[o[1],o[2],o[3],o[4],o[5]]:
end do:
end do:
#
for i1 from 1 by 1 to n do 
for j1 from 1 by 1 to n do 
for i2 from 1 by 1 to n do 
for j2 from 1 by 1 to n do 
for k from 1 by 1 to n do 
o[k] := 0:
end do:
o[i1] := o[i1]+1:
o[j1] := o[j1]+1:
o[i2] := o[i2]+1:
o[j2] := o[j2]+1:
iint_p_r := iint_p_r + (1/8)*C[i1,j1]*C[i2,j2]*r[o[1],o[2],o[3],o[4],o[5]]:
end do:
end do:
end do:
end do:

#now proof that lefthand side and right-hand side of (76) coincide
iint_p_r - int_p_r;