dat <- matrix(c(
  0,  0,  0,  1,  3,  4,  15,  7,
  0,  0,  0,  1,  5, 13,  11,  0,
  0,  1,  1,  2, 25, 15,   6,  0,
  0,  1,  3,  7, 14,  7,   4,  2,
  0,  1,  7, 15, 28,  8,   2,  0,
  0,  1,  7, 18, 15,  6,   0,  0,
  0,  4, 10, 12,  8,  2,   0,  0,
  0,  5, 11,  2,  3,  0,   0,  0,
  9, 12, 10,  3,  1,  0,   0,  0),
  nrow=9,ncol=8,byrow=T)
rownames(dat) <- 71:63
colnames(dat) <- seq(16.25,19.75,0.5)
rowtots <- apply(dat,1,sum)
coltots <- apply(dat,2,sum)
grandtot <- sum(dat)
print(dat)
freq <- as.vector(dat)
sta <- as.numeric(rep(rownames(dat),ncol(dat)))
cub <- as.numeric(rep(colnames(dat),rep(nrow(dat),ncol(dat))))
stature <- rep(sta,freq)
cubit <- rep(cub,freq)
mean.stature <- sum(sta*freq)/grandtot
mean.cubit <- sum(cub*freq)/grandtot
var.stature <- sum((sta-mean.stature)^2*freq)/(grandtot-1)
var.cubit <- sum((cub-mean.cubit)^2*freq)/(grandtot-1)
cov.stature.cubit <- sum((sta-mean.stature)*(cub-mean.cubit)*freq)/(grandtot-1)
cor.stature.cubit <- cor(stature,cubit)
cat("Correlation between cubit and stature is",cor(stature,cubit),"\n")
cat("Note while column totals are correctly given as\n")
print(coltots)
cat("some row totals are wrong - right values are\n")
print(rowtots)
cat("\nNeglecting extreme values data reduces to\n")
print(dat[2:8,2:7])
freq.co <- freq
freq.co[(sta<64)|(sta>70)|(cub>19.5)] <- 0
grandtot.co <- sum(freq.co)
stature.co <- rep(sta,freq.co)
cubit.co <- rep(cub,freq.co)
mean.stature.co <- sum(sta*freq.co)/grandtot.co
mean.cubit.co <- sum(cub*freq.co)/grandtot.co
var.stature.co <- sum((sta-mean.stature.co)^2*freq.co)/(grandtot.co-1)
var.cubit.co <- sum((cub-mean.cubit.co)^2*freq.co)/(grandtot.co-1)
cov.stature.cubit.co <- 
  sum((sta-mean.stature.co)*(cub-mean.cubit.co)*freq.co)/(grandtot.co-1)
cor.stature.cubit.co <- cor(stature.co,cubit.co)
cat("Corrected correlation between cubit and stature is",
    cor(stature.co,cubit.co),"\n")
d <- as.vector(t(dat))
cub.factor <- gl(8,1,72)
stat.factor <- gl(9,8,72)
cub.length <- 15.25 + as.numeric(cub.factor)
stat.length <- 70.5 - as.numeric(stat.factor)
cub.reps <- rep(cub.length,d)
stat.reps <- rep(stat.length,d)
plot(cub.length,stat.length,type="n",
  xlab="Left cubit", ylab="Stature")
text(cub.length,stat.length,d)