#7.9 Fat and Protein Content Example

x<-rbind(
c(3.30, 3.62),
c(3.24, 3.63),
c(3.25, 3.58),
c(3.24, 3.57),
c(3.33, 3.54),
c(3.27, 3.60),
c(3.23, 3.62),
c(3.26, 3.66),
c(3.26, 3.43),
c(3.31, 3.57),
c(3.27, 3.57),
c(3.23, 3.68),
c(3.27, 3.61),
c(3.27, 3.57),
c(3.26, 3.54),
c(3.29, 3.59),
c(3.29, 3.57),
c(3.26, 3.53),
c(3.23, 3.58),
c(3.24, 3.47),
c(3.31, 3.55),
c(3.20, 3.62),
c(3.29, 3.57),
c(3.21, 3.60))

#descriptive statistics
n<-dim(x)[1]
n
mu<-apply(x,2,mean)
mu
S<-sqrt(diag(cov(x)*(n-1)/n))
S
r<-(cov(x)[2,1]*(n-1)/n)/(S[1]*S[2])
r

#transform
y<-cbind(
(x[,1]-mu[1])/S[1],
((x[,2]-mu[2])/S[2]-r*(x[,1]-mu[1])/S[1])/sqrt(1-r^2))

#calc Kozoil's statistics
m<-array(0,dim=c(4+1,4+1))
for (i in 0:(nrow(m)-1))
for (j in 0:(ncol(m)-1))
{
m[i+1,j+1]<-sum((y[,1]^i)*(y[,2]^j))/n
}
U3sq<-n*((m[2+1,1+1]^2+m[1+1,2+1]^2)/2+(m[3+1,0+1]^2+m[0+1,3+1]^2)/6)
U4sq<-n*((m[2+1,2+1]-1)^2/4+(m[3+1,1+1]^2+m[1+1,3+1]^2)/6+((m[0+1,4+1]-3)^2+(m[4+1,0+1]-3)^2)/24)
U3sq+U4sq

#Hermite-Chebyshev polynomials
g<-function(n,y)
{
if (n==0) return(1)
if (n==1) return(y)
if (n==2) return((y^2-1)/sqrt(2))
if (n==3) return((y^3-y)/sqrt(6))
if (n==4) return((y^4-6*y^2+3)/sqrt(24))
if (n==5) return((y^5-10*y^3+15*y)/sqrt(120))
if (n==6) return((y^6-15*y^4+45*y^2-15)/sqrt(720))
if (n>6) print("generate more polynomials")
}

#Vrs's
Vrs<-array(0,dim=c(4+1,4+1))
for (i in 0:(nrow(Vrs)-1))
for (j in 0:(ncol(Vrs)-1))
{
Vrs[i+1,j+1]<-sum(g(i,y[,1])*g(j,y[,2]))/sqrt(n)
}

#U3sq and (not squared) components
U3sq
c(Vrs[3+1,0+1],Vrs[0+1,3+1],Vrs[2+1,1+1],Vrs[1+1,2+1])

#U4sq and (not squared) components
U4sq
c(Vrs[4+1,0+1],Vrs[0+1,4+1],Vrs[3+1,1+1],Vrs[1+1,3+1],Vrs[2+1,2+1])