dat <- as.matrix(read.table("galton_corr.dat"))
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")