The following image is seen on wikipedia when searching for Hotelling's T-squared distribution
This is apparently the pdf of the Hotelling T-squared distribution at different parameters. However, I am a bit suspicious of it, so I tried generate my own pdfs. Since,
$$T_{p,n-1}^2=\frac{p(n-1)}{n-p}F_{p,n-p}$$
where p =#dimensions, n=sample size and m= n-1
I generated some pdfs for the F-distribution with similar parameters, which can be seen in the following image.
The code for these density functions are
#PDF for the F-distribution with various degrees of freedom
labels <- c("p=2 m=5", "p=3 m=5", "p=4 m= 5", "p=4 m=50", "p=4 m=1000l")
colors <- c("red", "blue", "darkgreen", "gold", "cyan")
x <- seq(0,5, length=10000)
plot(x,df(x=x,df1=2,df=4), type='l', col ='red', xlab=' ', ylab= 'p(x)', main = 'Probability density function for the F-distribution', lwd=3, cex.main=0.9)
curve(df(x,df1=3,df2=3), from =0, to=5, col='blue',lwd=3, add=TRUE)
curve(df(x,df1=4,df2=2), from =0, to=5, col='darkgreen', add=TRUE, lwd=3)
curve(df(x,df1=4,df2=47), from =0, to=5, col='gold', add=TRUE, lwd=3)
curve(df(x,df1=4,df2=4997), from =0, to=5, col='cyan', add=TRUE, lwd=3)
axis(side = 1, lwd = 2)
axis(side = 2, lwd = 2)
legend("topright", inset=.05, title="Distributions",
labels, lwd=2, lty=c(1, 1, 1, 1, 2), col=colors)
Now, by fooling around a bit I was able to replicate the wikipedia image with the following code
#PDF for hotelling T-squared with various degrees of freedom
#Generated by transforming the F distribution
labels <- c("p=2 m=5", "p=3 m=5", "p=4 m= 5", "p=4 m=50", "p=4 m=5000")
colors <- c("black", "red", "green", "darkblue", "cyan")
x <- seq(0,15, length=10000)
plot(x,df(x=(1/2.5)*x,df1=2,df=4), type='l', col ='black', xlab=' ', ylab= 'p(x)', main = 'Probability density function for the \nHotelling T-squared distribution', lwd=3, cex.main=0.8)
curve(df((1/5)*x,df1=3,df2=3), from =0, to=15, col='red', add=TRUE, lwd=3)
curve(df((1/10)*x,df1=4,df2=2), from =0, to=15, col='green', add=TRUE, lwd=3)
curve(df((46/200)*x,df1=4,df2=47), from =0, to=15, col='darkblue', add=TRUE, lwd=3)
curve(df((4996/20000)*x,df1=4,df2=4997), from =0, to=15, col='cyan', add=TRUE, lwd=3)
axis(side = 1, lwd = 2)
axis(side = 2, lwd = 2)
legend("topright", inset=.05, title="Distributions",
labels, lwd=2, lty=c(1, 1, 1, 1, 1), col=colors, cex=0.75)
Getting the following image,
Which compared to the wikipedia is basically the same image. The thing I found strange was that examining the code, it appears that they have used the formula
$$T_{p,n-1}^2=\frac{n-p}{p(n-1)}F_{p,n-p)}$$
Instead of using the original formula. I believe the code should be
labels <- c("p=2 m=5", "p=3 m=5", "p=4 m= 5", "p=4 m=50", "p=4 m=5000")
colors <- c("black", "red", "green", "darkblue", "cyan")
x <- seq(0,2, length=10000)
plot(x,df(x=(2.5)*x,df1=2,df=4), type='l', col ='black', xlab=' ', ylab= 'p(x)', main = 'Probability density function for the \nHotelling T-squared distribution', lwd=3, cex.main=0.8)
curve(df((5)*x,df1=3,df2=3), from =0, to=15, col='red', add=TRUE, lwd=3)
curve(df((10)*x,df1=4,df2=2), from =0, to=15, col='green', add=TRUE, lwd=3)
curve(df((200/46)*x,df1=4,df2=47), from =0, to=15, col='darkblue', add=TRUE, lwd=3)
curve(df((20000/4996)*x,df1=4,df2=4997), from =0, to=15, col='cyan', add=TRUE, lwd=3)
axis(side = 1, lwd = 2)
axis(side = 2, lwd = 2)
legend("topright", inset=.05, title="Distributions",
labels, lwd=2, lty=c(1, 1, 1, 1, 1), col=colors, cex=0.75)
Producing the following image
This did not seem right as well. Thus, I tried generate the pdfs from hotelling T2 tests using the following methodology.
- Conduct a lot of hotelling T2 test
- Collect all the T2 statistics in a vector
- Determine their distributions
- Compare their distributions to $\frac{p(n-1)}{n-p}F_{p,n-p}$
I started off with the parameters p=4, m=50, n=51. Using the following code
library(ICSNP)
library(MVN)
T_Squared<-c()
for (i in 1:10000) {
x<-rnorm(51)
y<-rnorm(51)
z<-rnorm(51)
w<-rnorm(51)
x_bind<-cbind(x,y,z,w) #Create a four-sample vector
T<-HotellingsT2(x_bind,mu=c(0,0,0,0)) #Computing the Hotelling T-squated test
T2<-T$statistic #Save the T2 statistic and putting it into a vector
T_Squared[i]<-T2
}
hist(T_Squared, breaks='Scott', freq=FALSE, xlim=c(0,15), ylim=c(0,1), main = "Distribution of T2", xlab = ' ')
lines(density(T_Squared)) #Very similar to Hotelling T2 distriubtion
x<- seq(0,5, length=10000)
plot(x,df(x=x,df1=4,df2=47), type='l', col ='red', xlab=' ', ylab= 'p(x)', main = 'Probability density function for the F-distribution and T-squared', lwd=3, cex.main=0.9, ylim=c(0,1))
lines(density(T_Squared), col='blue')
legend("topright", inset=.05, title="Distributions",
labels, lwd=2, lty=c(1, 1), col=colors, cex=0.75)
As can be seen in the following image, the hotelling T-squared statistics have an identical distribution to the F-distribution.
Which is confirmed by conducting a KS-test
f<-rf(10000, df1=4, df2=47)
ks.test(T_Squared, f)
Two-sample Kolmogorov-Smirnov test
data: T_Squared and f
D = 0.0135, p-value = 0.3219
alternative hypothesis: two-sided
Perhaps the sample size is too large, after all Hotell T-squared test is most effective with smaller samples. I tried with parameters n=21, m=20, p=3 and got the following result
labels <- c("F_{3,18}", "T_{3,20}")
colors <- c("red", "blue")
T_Squared<-c()
for (i in 1:10000) {
x<-rnorm(21)
y<-rnorm(21)
z<-rnorm(21)
x_bind<-cbind(x,y,z) #Create a four-sample vector
T<-HotellingsT2(x_bind,mu=c(0,0,0)) #Computing the Hotelling T-squated test
T2<-T$statistic #Save the T2 statistic and putting it into a vector
T_Squared[i]<-T2
}
hist(T_Squared, breaks='Scott', freq=FALSE, xlim=c(0,15), ylim=c(0,1), main = "Distribution of T2", xlab = ' ')
lines(density(T_Squared)) #Very similar to Hotelling T2 distriubtion
x<- seq(0,5, length=10000)
plot(x,df(x=x,df1=3,df2=18), type='l', col ='red', xlab=' ', ylab= 'p(x)', main = 'Probability density function for the F-distribution and T-squared', lwd=3, cex.main=0.9, ylim=c(0,1))
lines(density(T_Squared), col='blue')
legend("topright", inset=.05, title="Distributions",
labels, lwd=2, lty=c(1, 1), col=colors, cex=0.75)
f<-rf(10000, df1=3, df2=18)
ks.test(T_Squared, f)
With the resulting p-value
Two-sample Kolmogorov-Smirnov test
data: T_Squared and f
D = 0.013, p-value = 0.3667
alternative hypothesis: two-sided
Last try now with n=6, m=5, p=2
labels <- c("F_{2,4}", "T_{2,6}")
colors <- c("red", "blue")
T_Squared<-c()
for (i in 1:10000) {
x<-rnorm(6)
y<-rnorm(6)
x_bind<-cbind(x,y) #Create a four-sample vector
T<-HotellingsT2(x_bind,mu=c(0,0)) #Computing the Hotelling T-squated test
T2<-T$statistic #Save the T2 statistic and putting it into a vector
T_Squared[i]<-T2
}
hist(T_Squared, breaks='Scott', freq=FALSE, xlim=c(0,15), ylim=c(0,1), main = "Distribution of T2", xlab = ' ')
lines(density(T_Squared)) #Very similar to Hotelling T2 distriubtion
x<- seq(0,5, length=10000)
plot(x,df(x=x,df1=2,df2=4), type='l', col ='red', xlab=' ', ylab= 'p(x)', main = 'Probability density function for the F-distribution and T-squared', lwd=3, cex.main=0.9, ylim=c(0,1))
lines(density(T_Squared), col='blue')
legend("topright", inset=.05, title="Distributions",
labels, lwd=2, lty=c(1, 1), col=colors, cex=0.75)
f<-rf(10000, df1=2, df2=4)
ks.test(T_Squared, f)
Two-sample Kolmogorov-Smirnov test
data: T_Squared and f
D = 0.0124, p-value = 0.4255
alternative hypothesis: two-sided
Same result as the last times. Why is this the case? I expected the hotelling T-sqaured distribution to have similar pdfs as the F-distribution, but not identically the same. Is my code completely off? Is wikipedia right? Can you generate the pdf from Hotelling T-sqaured distribution from condcuting hotelling T2 tests?
If anyone has any insight on this I would really appreciate some clarification
