> #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
[1] 24
> mu<-apply(x,2,mean)
> mu
[1] 3.262917 3.577917
> S<-sqrt(diag(cov(x)*(n-1)/n))
> S
[1] 0.03194777 0.05291339
> r<-(cov(x)[2,1]*(n-1)/n)/(S[1]*S[2])
> r
[1] -0.2280972
> 
> #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
[1] 11.96056
> 
> #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
[1] 5.861887
> c(Vrs[3+1,0+1],Vrs[0+1,3+1],Vrs[2+1,1+1],Vrs[1+1,2+1])
[1]  0.2046187 -2.1136136  0.5206869 -1.0399714
> 
> #U4sq and (not squared) components
> U4sq
[1] 6.098675
> 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])
[1] -0.5408338  1.7912501 -0.4006503  0.4833948 -1.4843873
>