x <- rnorm(100)
y <- 1 + x + rnorm(100)
y <- as.integer(y<0)
probitfit <- glm(y ~ x, family = binomial(link = probit))
summary(probitfit)
index <- probitfit$coef[1] + probitfit$coef[2]*x
invmills <- dnorm(index)/pnorm(index)