x<-c(8,14,-7,6,5,6,-5,1,0,-17)
y<-c(1,1,0,0,1,0,1,0,0,0)
n<-length(y)
X<-matrix(c(rep(1,n),x),ncol=2,byrow=FALSE)
X

###### Fisher Scoring
###### Fisher Scoring
###### Fisher Scoring
###### Fisher Scoring
###### Fisher Scoring


# Numero di variabili indipendenti (colonne matrice X – 1)
j=1
# Numero progressivo di iterazioni
i=0
# Stime iniziali dei parametri ignoti b0 e b1
phit0<-matrix(c(0,0)) # le probabilità previste sono 1/2 per ogni valore di x: exp(0)/(1+exp(0)) = 1/2

# Grandezza del campione
n<-length(y)
# data set vuoto per i risultati:
FisherScoring<-c()




repeat{
           
           i=i+1 # prima iterazione (la prima sarà nulla, quindi registreremo i-1)
		
	     B<-phit0[1:(j+1)]
             p<-exp(X%*%B)/(1+exp(X%*%B))
             u<-t(X)%*%(matrix(y)-p)
             
             W<-diag( c(-p*(1-p)) )    
             H<-t(X)%*%W%*%X

             phit1<- phit0 - solve(H)%*%u
                         
             #registriamo i risultati nel data set FisherScoring:
             Result_iter<-c(phit1[1:(j+1)])
             FisherScoring<-cbind(FisherScoring, Result_iter)
             
             #condizione logica per interrompere il ciclo:
             if(all(abs(phit1 - phit0) < .00001)){
             	                                 rownames(FisherScoring)<-c(paste("b",seq(0,j)))
             	                                 colnames(FisherScoring)<-seq(0,(i-1))
             	                                 print(FisherScoring)
             	                                 break
                                                 }
	  phit0<-phit1
}



###### Metodi alternativi di ottimizzazione
###### Metodi alternativi di ottimizzazione
###### Metodi alternativi di ottimizzazione
###### Metodi alternativi di ottimizzazione
###### Metodi alternativi di ottimizzazione

#Dati:
x<-c(8,14,-7,6,5,6,-5,1,0,-17)
y<-c(1,1,0,0,1,0,1,0,0,0)
n<-length(y)
j=1 # Numero di variabili indipendenti
X<-matrix(c(rep(1,n),x),ncol=2,byrow=FALSE)
X

#Log Likelihood function con argomento «phi» = vettore di parametri ignoti.
# Versione 1: Eq. 48
Log_Lik<-function(phi){
                      B<-phi[1:(j+1)]
                      score=sum( y*X%*%B - log(1+exp(X%*%B)) )
                      (-1)*score # algoritmi di minimizzazione!
}
# Versione 2: Eq. 38
Log_Lik2<-function(phi){
                      B<-phi[1:(j+1)]
                      p<-exp(X%*%B)/(1+exp(X%*%B))
                      score=sum( y*log(p/(1-p)) + log(1-p) )
                      (-1)*score # algoritmi di minimizzazione!
}

phit0<-c(0,0) #vettore di stime iniziali 

# Versione 1: Eq. 48
optim(par=phit0, fn=Log_Lik, method="Nelder-Mead")
optim(par=phit0, fn=Log_Lik, method="BFGS")
optim(par=phit0, fn=Log_Lik, method="CG")
optim(par=phit0, fn=Log_Lik, method="L-BFGS-B")
optim(par=phit0, fn=Log_Lik, method="SANN")
# Versione 2: Eq. 38
optim(par=phit0, fn=Log_Lik2, method="Nelder-Mead")
optim(par=phit0, fn=Log_Lik2, method="BFGS")
optim(par=phit0, fn=Log_Lik2, method="CG")
optim(par=phit0, fn=Log_Lik2, method="L-BFGS-B")
optim(par=phit0, fn=Log_Lik2, method="SANN")



############# GLM 
############# GLM 
############# GLM 
############# GLM 
############# GLM 

summary(glm(y~x,family=binomial(link=logit)))

-1*solve(H)
##             [,1]        [,2]
## [1,]  0.66149431 -0.04421518        # Varianza distorta di (Intercetta) B_0
## [2,] -0.04421518  0.01333360        # Varianza distorta di B_1

vcov.B<-(-1*solve(H))                   # var/covarianze (distorte) dei coefficienti B
diag(vcov.B)                            # diagonale delle varianze (distorte) dei coefficienti B
sqrt(diag(vcov.B))                      # errori standard dei coefficienti B


############# Residui glm()
#Deviance residuals...
residuals(glm(y~x,family=binomial),type="deviance")
# manualmente
sign(y-p)*sqrt(-2*(y*X%*%B - log(1+exp(X%*%B))))

#response residuals...
residuals(glm(y~x,family=binomial),type="response")
# manualmente
y-p

# Pearson (standardized) residuals...
residuals(glm(y~x,family=binomial),type="pearson")
# manualmente
((y-p)-0)/sqrt(p*(1-p))
#############





#################### Sparse Data problem x casa
################### Sparse Data problem x casa
################# Sparse Data problem x casa
x1 = c(0.5, 2, 10, 6, 8, 1, 8, 4, 0, 7, 5, 1.1)
x2 = c(7, 9.5, 6, 0, 8.1, 9, 10, 0, 8, 5, 3, 4)
y = c(0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0)