clear all
set seed 12345
set obs 30

g x  = runiform()
g u  = rnormal()
g y  = .5 + x + u > 0
g um = 1
g p  = .

mkmat um x, mat(X)
mkmat y   , mat(Y)

mat beta = inv(X'*X)*X'*Y

forv k = 1/20 {
 sca a = beta[1,1]
 sca b = beta[2,1]
 replace p = exp(a+b*x)/(1+exp(a+b*x))
 mkmat p, mat(P)
 mat grad = X'*(Y-P)
 mat I = J(30,1,1)
 mat W = diag(hadamard(P, I-P))
 mat hessiano = -1*X'*W*X
 mat delta = inv(hessiano)*grad
 mat dd = delta'*delta
  if sqrt(dd[1,1]) < 1e-5 continue, break
 mat beta = beta - delta
 sca k = `k'
}

* Resultado final
sca a = beta[1,1]
sca b = beta[2,1]
replace p = exp(a+b*x)/(1+exp(a+b*x))
g li = ln(y*p + (1-y)*(1-p))
qui su li, meanonly
sca lnL = r(sum)
mkmat p, mat(P)
mat hessiano = -1*X'*diag(hadamard(P, I-P))*X

mat list beta
di "lnL = " lnL
di "dp:  " sqrt(-hessiano[1,1]) "  " sqrt(-hessiano[2,2])
di "z:   " beta[1,1]/sqrt(-hessiano[1,1]) "  " beta[2,1]/sqrt(-hessiano[2,2])

logit y x    // verificação

qui su y, meanonly
sca ybar = r(mean)
sca lnL0 = 30*(ybar*log(ybar) + (1-ybar)*log(1-ybar))
sca pseudoR2 = 1 - lnL/lnL0

di "lnL0     = " lnL0
di "pseudoR2 = " pseudoR2