# Simple Linear Regression. Simulate N points from 5 number summary
SD = function(x) {
	n = length(x)
	sqrt((n-1)/n)*sd(x)
}

su = function(x) (x-mean(x))/SD(x)

getxy = function(N = 928, mx = 68.31, my = 68.09, sx = 1.79, sy = 2.52, r = 0.46) {
	x = rnorm(N,mx,sx)
	y = rnorm(N,my + (x-mx)/sx*r*sy,sqrt(1-r^2)*sy)
	plot(x,y,main=paste(N,"heights of Sons vs. Fathers"), xlab="Father's height [in]",
         ylab="Sons's height [in]")
	list(x=x,y=y)
}

xy = getxy()
x = xy$x
y = xy$y

ptaves = function(mx = mean(x),my = mean(y)) {
	abline(v=mx,col="red")
	abline(h=my,col="red")
}

sdline = function(mx=mean(x),my=mean(y),sx=sd(x),sy=sd(y),r=cor(x,y)) {
	b = (sy/sx)*sign(r)
	a = my - b*mx
	abline(a,b,col="red")
	list(a=a,b=b)
}

regline = function(mx=mean(x),my=mean(y),sx=sd(x),sy=sd(y),r=cor(x,y)) {
	a = my - r*sy/sx*mx
	b = r*sy/sx
	abline(a,b,col="blue")
	list(a=a,b=b)
}

corr = function(x,y) { mean(su(x)*su(y)) }


xys = function(N = 10, mx = 0, my = 0, sx = 1, sy = 1, r = 0.5,...) {
	x = rnorm(N,mx,sx)
	y = rnorm(N,my + (x-mx)/sx*r*sy,sqrt(1-r^2)*sy)
	plot(x,y,...)
	list(x=x,y=y)
}


add.ellipse = function(XY=xy,...) {
    x <- XY$x
    y <- XY$y
    require("ellipse")
    Z <- ellipse(cor(x,y))
    Elli <- cbind(mean(x)+SD(x)*Z[,1],mean(y)+SD(y)*Z[,2])
    points(Elli,type="l",...)
}