###################################################
# Generiamo i dati per la regressione lineare #####

  X<-cbind(1,runif(100))

# Valore dei coefficienti ignoti [2,3] e della varianza dell'errore [1]
  phi.true<-c(2,3,1)
  y<-X%*%phi.true[1:2] + rnorm(100,0,phi.true[3]) 
###################################################

fisher.scrng = function(y,X,tolerance=1e-12,max.iter=1000) {
# primo livello, preparazione del vettore e della log verosimiglianza iniziale, che verranno via via aggiornati
 n = nrow(X)
 j = ncol(X)-1
phi0 = runif(ncol(X)+1) # stima iniziale random-uniforme (t=0). Usate nel primo ciclo while, poi aggiornate
iterations = 0
cambiamento = TRUE # C'è differenza nella log Lik tra la stima (t) e (t+1)
result<-c()

		while(cambiamento & (iterations < max.iter)) { # Finchè c'è cambiamento nella log Lik...
		iterations <- iterations +1
		cambiamento <- FALSE #

	      # Stime iniziali e iterazione (t-esima):
		beta<-phi0 [1:(j+1)]
	      sigq<-phi0 [(j+1)+1]
		eps<-y-X%*%beta
		logLik_Old<- -n/2*log(2*pi)-n/2*log(sigq)-1/(2*sigq)*(t(eps)%*%eps)

		# Calcolo delle derivate prime e della matrice Hessiana, secondo i valori MLE
		u<-c(
			c(1/sigq*t(X)%*%y-1/sigq*t(X)%*%X%*%beta), # derivata prima betas
 			-n/(2*sigq) + (1/(2*sigq^2)) * crossprod(y - X %*% beta)# derivata prima varianza errore
		)  
 
#Hessian.Exp <- matrix(0,ncol=j+2,nrow = j+2)
#for(z in 1:n)
#{
# Estimate derivative of likelihood for each observation
#estVec <- c((1/sigq) * 
#t((y[z] - X[z,]%*%beta)) %*% X[z,],
#-(n/(2*sigq)) + (1/(2*sigq^2)) *
#crossprod(y[z] - X[z,] %*% beta))
# Add them up as suggested to get an estimate of the Hessian.
#Hessian.Exp <- Hessian.Exp - estVec%o%estVec
#    }

		Hessian.Exp<-matrix(0, nrow=j+2,ncol=j+2)
	 	Hessian.Exp[1:(j+1),1:(j+1)]= -1/(sigq)*t(X)%*%X
 		Hessian.Exp[(j+2),(j+2)]= -n/(2*sigq^2)

		# Stime aggiornate (t+1) e funzione di Log Lik:
		new.phi = phi0  - solve(Hessian.Exp)%*%u
	      beta_new<-new.phi [1:(j+1)]
	      sigq_new<-new.phi [(j+1)+1]
		eps_new<-y-X%*%beta_new
		logLik_new = -n/2*log(2*pi)-n/2*log(sigq_new)-1/(2*sigq_new)*(t(eps_new)%*%eps_new)
		#cambiamento.relativo = abs( (logLik_new - logLik_Old)/logLik_Old )
            cambiamento.relativo = logLik_new - logLik_Old
		cambiamento= (cambiamento.relativo > tolerance)
		phi0 = new.phi
            result=rbind(result,c(phi0,logLik_new,cambiamento.relativo))

		}
if (cambiamento) {# nel caso non abbia ancora trovato un cambiamento < tolerance, e si necessiti di aumentare max.iter
warning("Fisher scoring ha terminato prima della convergenza")
}
return(list(minimum=new.phi,value=logLik_new,deriv=u,deriv2=Hessian.Exp,
iterations=iterations,converged=!cambiamento,result=result))
}

##########################################################
# Applicazione pratica
fisher.scrng(y,X)
##########################################################

Ottimizzazione in R, con derivate ed hessiane numeriche
 logLik<-function(phi,y,X){
 n = nrow(X)
 j = ncol(X)-1
 beta<-phi [1:(j+1)]
 sigq<-phi [(j+1)+1]
 eps<-y-X%*%beta
 log_Lik<- -n/2*log(2*pi)-n/2*log(sigq)-1/(2*sigq)*(t(eps)%*%eps)
 return(-log_Lik)
}

optim(c(1,1,1),logLik,method="BFGS",hessian=T,y=y,X=X,control = list(trace = TRUE, REPORT = 1))


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

 X<-cbind(1,c(2,3,1,4,5,8))
 y=c(3,2,2,7,6,7)

 optim(c(1,1,1),logLik,method="BFGS",hessian=T,y=y,X=X,control = list(trace = TRUE, REPORT = 1))
 fisher.scrng(y,X)

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


glm(y~X[,2],gaussian(link = "identity"))