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


# we show closed-form solutions for the unbounded integrals in existing and new 
# uncertainty aware PCA
# we do for a concrete number of samples, and concrete dimensionality, feel free to change them

#number of samples 
N := 5;

#dimensionality of data points
n := 4;


#dimensionality of Covariance matrices
r := n*(n+1)/2;

#(3), (4)
T := Matrix(r,2):
d := Vector(r):
ii := 1:
for i from 1 by 1 to n do 
for j from i by 1 to n do 
T[ii,1] := i:
T[ii,2] := j:
ii := ii+1:
end do:
end do:
for i from 1 by 1 to r do 
if T[i,1]=T[i,2] then d[i]:=1: else d[i]:=sqrt(2): end if:
end do:
T;
d;

#coordinates of data points 
for i from 1 by 1 to N do 
Xj[i] := Vector(n,symbol=x[i]):
end do:

#local Mean for each data point (13)
for i from 1 by 1 to N do 
Mj[i] := Vector(n,symbol=m[i]):
end do:

#local Covariance for each data point (13)
for i from 1 by 1 to N do 
Cj[i] := Matrix(n,n,symbol=c[i]):
end do:

#set local Mean and Covariance to integer randoms and make sure that
#local Covariances are positive definite 
for i from 1 by 1 to N do 
for j from 1 by 1 to n do 
Mj[i][j] := (roll()-10)/10:
end do:
end do:
#
for i from 1 by 1 to N do 
HCj[i] := Matrix(n,n,symbol=hc[i]):
end do:
#
for ii from 1 by 1 to N do 
for i from 1 by 1 to n do 
for j from 1 by 1 to n do 
HCj[ii][i,j] := (roll()-10)/10:
end do:
end do:
end do:
#
for ii from 1 by 1 to N do 
Cj[ii] := Multiply( Transpose(HCj[ii]), HCj[ii]):
end do:

#plot all the input data
for i from 1 by 1 to N do 
print(Xj[i], Mj[i],  Cj[i]);
end do:

####

#compute input normal distributions (13)
for i from 1 by 1 to N do 
CjI[i] := MatrixInverse(Cj[i]):
end do:
#
for i from 1 by 1 to N do 
ff[i] := -(1/2)*simplify(Multiply(Multiply(Transpose(Xj[i]-Mj[i]),CjI[i]),Xj[i]-Mj[i])):
end do:
#
for i from 1 by 1 to N do 
pf[i] := (1/sqrt((2*Pi)^n*Determinant(Cj[i])))*exp( ff[i]   ):
end do:
#print all input distributions (13)
for i from 1 by 1 to N do 
print(pf[i]);
end do:





###

#(14)
Mbar := (1/N)*sum( Mj[iii] , iii=1..N );
Cbar := (1/N)*sum( Cj[iii] , iii=1..N );

#(15)
C_M := Matrix(n,n):
for i from 1 by 1 to N do 
C_M := C_M +  Multiply(Mj[i]-Mbar, Transpose(Mj[i]-Mbar)) 
end do:
C_M := C_M/N:
C_M;

#(21)
Mtilde :=  (1/N)*sum( Xj[iii] , iii=1..N );

#compute (22) by symbolic integration
M := Vector(n):
for i from 1 by 1 to n do 
j := 1:
h[j] := Mtilde[i]:
for ii from 1 by 1 to N do 
h[j+1] := h[j]*pf[ii]:
j := j+1:
end do:
for ii from 1 by 1 to N do 
for jj from 1 by 1 to n do 
h[j+1] := int( h[j] , Xj[ii][jj]=-infinity..infinity):
j := j+1:
end do:
end do:
M[i] := h[j]:
end do:
M;

#(23)
Ctilde_goertler := Matrix(n,n):
for i from 1 by 1 to N do 
Ctilde_goertler := Ctilde_goertler +  Multiply(Xj[i]-M, Transpose(Xj[i]-M)) 
end do:
Ctilde_goertler := Ctilde_goertler/N:
Ctilde_goertler;

#compute (25) by symbolic integration
C_goertler := Matrix(n,n):
for i from 1 by 1 to n do 
for k from 1 by 1 to n do 
j := 1:
h[j] := Ctilde_goertler[i,k]:
for ii from 1 by 1 to N do 
h[j+1] := h[j]*pf[ii]:
j := j+1:
end do:
for ii from 1 by 1 to N do 
for jj from 1 by 1 to n do 
h[j+1] := int( h[j] , Xj[ii][jj]=-infinity..infinity):
j := j+1:
end do:
end do:
C_goertler[i,k] := h[j]:
end do:
end do:
C_goertler;

#show (16)
simplify(M-Mbar);

#show (17)
simplify(C_goertler - C_M - Cbar);


#we now switch to the new uncertainty PCA

#(40)
Ctilde := Matrix(n,n):
for i from 1 by 1 to N do 
Ctilde := Ctilde +  Multiply(Xj[i]-Mtilde, Transpose(Xj[i]-Mtilde)) 
end do:
Ctilde := Ctilde/N:
Ctilde;

#(51) in matrix form
C := Matrix(n,n):
for i from 1 by 1 to n do 
for k from 1 by 1 to n do 
j := 1:
h[j] := Ctilde[i,k]:
for ii from 1 by 1 to N do 
h[j+1] := h[j]*pf[ii]:
j := j+1:
end do:
for ii from 1 by 1 to N do 
for jj from 1 by 1 to n do 
h[j+1] := int( h[j] , Xj[ii][jj]=-infinity..infinity):
j := j+1:
end do:
end do:
C[i,k] := h[j]:
end do:
end do:
C;

#(42)
VCtilde := Vector(r):
for i from 1 by 1 to r do 
VCtilde[i]:= d[i]*Ctilde[ T[i,1] , T[i,2]]:
end do:
VCtilde;


#(50)
M := Vector(n):
for i from 1 by 1 to n do 
j := 1:
h[j] := Mtilde[i]:
for ii from 1 by 1 to N do 
h[j+1] := h[j]*pf[ii]:
j := j+1:
end do:
for ii from 1 by 1 to N do 
for jj from 1 by 1 to n do 
h[j+1] := int( h[j] , Xj[ii][jj]=-infinity..infinity):
j := j+1:
end do:
end do:
M[i] := h[j]:
end do:
M;

#(51)
M_C := Vector(r):
for i from 1 by 1 to r do 
j := 1:
h[j] := VCtilde[i]:
for ii from 1 by 1 to N do 
h[j+1] := h[j]*pf[ii]:
j := j+1:
end do:
for ii from 1 by 1 to N do 
for jj from 1 by 1 to n do 
h[j+1] := int( h[j] , Xj[ii][jj]=-infinity..infinity):
j := j+1:
end do:
end do:
M_C[i] := h[j]:
end do:
M_C;

#(52)
MMtilde := Multiply( Mtilde-M, Transpose(Mtilde-M));

#(53)
C_C_tilde := Multiply( VCtilde - M_C, Transpose(VCtilde - M_C));

#(54)
MM := Matrix(n,n):
for i from 1 by 1 to n do 
for k from 1 by 1 to n do 
j := 1:
h[j] := MMtilde[i,k]:
for ii from 1 by 1 to N do 
h[j+1] := h[j]*pf[ii]:
j := j+1:
end do:
for ii from 1 by 1 to N do 
for jj from 1 by 1 to n do 
h[j+1] := int( h[j] , Xj[ii][jj]=-infinity..infinity):
j := j+1:
end do:
end do:
MM[i,k] := h[j]:
end do:
end do:
MM;

#(55) - this may take a few minutes as many integrations are involved
C_C := Matrix(r,r):
for i from 1 by 1 to r do 
for k from 1 by 1 to r do 
j := 1:
h[j] := C_C_tilde[i,k]:
for ii from 1 by 1 to N do 
h[j+1] := h[j]*pf[ii]:
j := j+1:
end do:
for ii from 1 by 1 to N do 
for jj from 1 by 1 to n do 
h[j+1] := int( h[j] , Xj[ii][jj]=-infinity..infinity):
j := j+1:
end do:
end do:
C_C[i,k] := h[j]:
end do:
end do:
C_C;

#show (56)
simplify(M-Mbar);

#show (57)
simplify( MM - (1/N)*Cbar );

#show (58)
simplify(C - C_M - ((N-1)/N)*Cbar);


VC := Vector(r):
for i from 1 by 1 to r do 
VC[i]:= d[i]*C[ T[i,1] , T[i,2]]:
end do:
#show (59)
M_C - VC;


for i from 1 by 1 to N do 
Mjhat[i] := Mj[i]-Mbar:
end do:

#(70)
C_C_2 := Matrix(r,r):
for j from 1 by 1 to r do 
for k from 1 by 1 to r do 
s2 := 
+Cbar[T[j,1],T[k,1]]*Cbar[T[j,2],T[k,2]]
+Cbar[T[j,1],T[k,2]]*Cbar[T[j,2],T[k,1]]
+Cbar[T[j,2],T[k,2]]*Cbar[T[j,1],T[k,1]]
+Cbar[T[j,2],T[k,1]]*Cbar[T[j,1],T[k,2]]:
s3 := 0:
for i from 1 by 1 to N do 
s3 := s3 
+Cj[i][T[j,1],T[k,1]]*Mjhat[i][T[j,2]]*Mjhat[i][T[k,2]]
+Cj[i][T[j,1],T[k,2]]*Mjhat[i][T[j,2]]*Mjhat[i][T[k,1]]
+Cj[i][T[j,2],T[k,2]]*Mjhat[i][T[j,1]]*Mjhat[i][T[k,1]]
+Cj[i][T[j,2],T[k,1]]*Mjhat[i][T[j,1]]*Mjhat[i][T[k,2]]:
end do:
s4 := 0:
for i from 1 by 1 to N do 
s4 := s4 
+Cj[i][T[j,1],T[k,1]]*Cj[i][T[j,2],T[k,2]]
+Cj[i][T[j,1],T[k,2]]*Cj[i][T[j,2],T[k,1]]
+Cj[i][T[j,2],T[k,2]]*Cj[i][T[j,1],T[k,1]]
+Cj[i][T[j,2],T[k,1]]*Cj[i][T[j,1],T[k,2]]:
end do:
s1 := (1/2)*s2 + s3 + ((N-2)/(2*N))*(s4):
C_C_2[j,k] := ((d[j]*d[k])/N^2)*s1:
end do:
end do:
C_C_2;

#show (55)=(70)
C_C - C_C_2;