# 2015.10.08

library(MASS)

ds=menarche; head(ds); tail(ds)

######  GLM

probit.fit = 
   glm(cbind(Menarche, Total - Menarche) ~ Age,
         family=binomial(link = probit), 
         data = menarche )
probit.fit
        
# 50:50 point (median and mean) probit = Z = 0

# solve z = a + b*age = 0

# age = -a/b 

age.50 = -probit.fit$coefficients[1]/
          probit.fit$coefficients[2]   ; age.50

sd = 1/probit.fit$coefficients[2] ; sd


#####  logit

logit.fit = 
   glm(cbind(Menarche, Total - Menarche) ~ Age,
         family=binomial(link = logit), 
         data = menarche )
logit.fit
         
# 50:50 point (median and mean) logit = Z = 0

# solve logit = a + b*age = 0 

age.50 = -logit.fit$coefficients[1]/
          logit.fit$coefficients[2]   ; age.50

sd = ( pi/sqrt(3) ) /logit.fit$coefficients[2] ; sd

##############


####    own    (normal for now, but east to replace 
########        by gamma or logNormal)

loglik=function(par) {
	mu=par[1]
	sigma=par[2]
	p = pnorm( (menarche$Age-mu)/sigma )
	loglik= sum( menarche$Menarche * log(p) +
	(menarche$Total - menarche$Menarche) * log(1-p)
	)
	return(loglik)
}

par=c(12,2)
fit=optim(par,loglik,control=list(fnscale=-1),hessian=TRUE)

varcov = solve(-fit$hessian)

sqrt(diag(varcov))