nhl_htwt <- read.csv("http://www.stat.ufl.edu/~winner/data/nhl_ht_wt.csv") attach(nhl_htwt); names(nhl_htwt) set.seed(1234) height_eps <- runif(length(Height)) Height <- Height + (height_eps - 0.5) htwt <- data.frame(Height, Weight) detach(nhl_htwt); attach(nhl_htwt) X <- cbind(Height, Weight) N.pop <- length(Height) one_N <- rep(1,N.pop) I_N <- diag(N.pop) J_N <- matrix(rep(1,N.pop^2),ncol=N.pop) (mu <- (1/N.pop) * t(X) %*% one_N) (Sigma <- (1/N.pop) * (t(X) %*% (I_N - (1/N.pop)*J_N) %*% X)) ### Plot Bivariate data with ellipse containing a proportion (.80 here) of the data points install.packages("car") library(car) xlimhi <- 1.2*ceiling(max(abs(Height))) ylimhi <- 1.2*ceiling(max(abs(Weight))) xlimlo <- 0.8*floor(min(abs(Height))) ylimlo <- 0.8*floor(min(abs(Weight))) dataEllipse(Height,Weight, levels=0.80,xlim=c(xlimlo,xlimhi), ylim=c(ylimlo,ylimhi),pch=16,cex=0.7) ############# ### Generate randoms samples n.samp <- 25 Num.samp <- 100000 one_n <- rep(1,n.samp) I_n <- diag(n.samp) J_n <- matrix(rep(1,n.samp^2),ncol=n.samp) Xbar <- matrix(rep(0,2*Num.samp), ncol=2) S <- matrix(rep(0,3*Num.samp), ncol=3) for (i in 1:Num.samp) { X.samp.rows <- sample(1:N.pop,n.samp,replace=F) X.samp <- X[X.samp.rows,] Xbar.samp <- (1/n.samp) * t(X.samp) %*% one_n S.samp <- (1/(n.samp-1)) * (t(X.samp) %*% (I_n - (1/n.samp)*J_n) %*% X.samp) Xbar[i,] <- t(Xbar.samp) S[i,1] <- S.samp[1,1]; S[i,2] <- S.samp[1,2]; S[i,3] <- S.samp[2,2] } colMeans(Xbar) colMeans(S) cov(Xbar) Sigma/n.samp ((Sigma/n.samp)*((N.pop-n.samp)/(N.pop-1))) # V(X-bar w/ FPCF)