logistic <- function(u,y,deg,niter){
#
#
# This function performs Newton-Raphson iterations to find the maximum
# likelihood estimates in a logistic regression model.  It applies to a
# situation with a single covariate, and models theta as a polynomial
# of degree deg.  The covariate values are u and the responses y.  The
# desired number of Newton-Raphson iterations is niter.  
#
# The function returns all the iterates in a deg+1 by niter+1 matrix.
# The design matrix is constructed using the function poly, which
# produces values of orthogonal polynomials.
#
#
n=length(u)
X=poly(u,deg)
ymod=y
ymod[ymod==0]=0.01
ymod[ymod==1]=0.99
z=log(ymod/(1-ymod))
beta=matrix(lsfit(X,z)$coef,deg+1,1)
X=cbind(rep(1/sqrt(n),len=n),X)
iterates=matrix(0,deg+1,niter+1)
iterates[,1]=beta[,1]
for(i in 1:niter){
theta=beta[1,1]*X[,1]
for(j in 1:deg){
theta=theta+beta[(j+1),1]*X[,(j+1)]
}
bpp=exp(theta)/(1+exp(theta))^2
V=diag(bpp)
mu=exp(theta)/(1+exp(theta))
S=matrix(y-mu,n,1)
new=solve(t(X) %*% V %*% X)
new=new %*% t(X) %*% S
beta = beta + new
iterates[,i+1]=beta[,1]
}
iterates
}