url <- "http://stats191.stanford.edu/data/flu.table"
flu.table <- read.table(url, header =T)
g <- function(x) {
	return(log(x/(1 - x)))
}

g.inv <- function(y) {
	return(exp(y)/(1 + exp(y)))
}

g.prime<-function(x) {
	return(1/(x * (1 - x)))
}

V<-function(mu) {
	return(mu * (1 - mu))
}

niter <- 10
mu <- rep(mean(flu.table$Shot), length(flu.table$Shot))

for (i in 1:niter) {
	Z <- g(mu) + g.prime(mu) * (flu.table$Shot - mu)
      W <- 1/(g.prime(mu)^2 * V(mu))
      Z.lm <- lm(Z ~ Age + Health.Aware, weights = W, data = flu.table)
      eta <- predict(Z.lm)
      mu <- g.inv(eta)
	beta <- Z.lm$coef
}

irls<-function(glink,ginv,gdot,Vfunc,X,Y,nits=10){
	mu <- rep(mean(Y),length(Y))
	for (i in 1:nits) {
		Z <- glink(mu) + gdot(mu) * (Y - mu)
      	W <- 1/(gdot(mu)^2 * Vfunc(mu))
      	Z.lm <- lm(Y ~ X, weights = W)
      	eta <- predict(Z.lm)
      	mu <- ginv(eta)
		beta <- Z.lm$coef
	}
	return(beta)
}

irls(g,g.inv,g.prime,V,test.data$X1,test.data$Y,20)
glm(Y~X1,data=test.data)