x<-c(2,3,1,4,5,8) y<-c(3,2,2,7,6,7) n<-length(y) X<-matrix(c(rep(1,6),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 phit0<-matrix(c(1,1,1)) # 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)] sigma_q<-phit0[j+2] u<-matrix(c( 1/sigma_q*t(X)%*%y-1/sigma_q*t(X)%*%X%*%B, -n/(2*sigma_q)+1/(2*sigma_q^2)*t(y-X%*%B)%*%(y-X%*%B) ) ,ncol=1) H<-matrix(NA,ncol=(j+2),nrow=(j+2)) H[1:(j+1),1:(j+1)]<--1/sigma_q*t(X)%*%X H[j+2,j+2]<--n/(2*sigma_q^2) H[3,1:(j+1)]<-rep(0,2) # oppure anche solo 0 H[1:(j+1),3]<-rep(0,2) # oppure anche solo 0 phit1<- phit0 - solve(H)%*%u #registriamo i risultati nel data set FisherScoring: Result_iter<-c(phit1[1:(j+1)],phit1[j+2]) FisherScoring<-cbind(FisherScoring, Result_iter) #condizione logica per interrompere il ciclo: if(sum(abs(phit1 - phit0) < .0001)==3){ rownames(FisherScoring)<-c(paste("b",seq(1,j+1)),"sigma_q") colnames(FisherScoring)<-seq(0,(i-1)) print(FisherScoring) break } phit0<-phit1 } ###### Modificato (con "sbaglio") ###### Modificato (con "sbaglio") ###### Modificato (con "sbaglio") ###### Modificato (con "sbaglio") ###### Modificato (con "sbaglio") # Numero di variabili indipendenti (colonne matrice X – 1) j=1 # Numero progressivo di iterazioni i=0 # Stime iniziali dei parametri ignoti phit0<-matrix(c(1,1,1)) # 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)] sigma_q<-phit0[j+2] u<-matrix(c( 1/sigma_q*t(X)%*%y-1/sigma_q*t(X)%*%X%*%B, -n/(2*sigma_q)+1/(2*sigma_q^2)*t(y-X%*%B)%*%(y-X%*%B) ) ,ncol=1) H<-matrix(NA,ncol=(j+2),nrow=(j+2)) H[1:(j+1),1:(j+1)]<--n/c(t(y-X%*%B)%*%(y-X%*%B))*t(X)%*%X # INSERITO LO SBAGLIO NELLA MATRICE HESSIANA H[j+2,j+2]<--n/(2*sigma_q^2) H[3,1:(j+1)]<-rep(0,2) # oppure anche solo 0 H[1:(j+1),3]<-rep(0,2) # oppure anche solo 0 phit1<- phit0 - solve(H)%*%u #registriamo i risultati nel data set FisherScoring: Result_iter<-c(phit1[1:(j+1)],phit1[j+2]) FisherScoring<-cbind(FisherScoring, Result_iter) #condizione logica per interrompere il ciclo: if(sum(abs(phit1 - phit0) < .0001)==3){ rownames(FisherScoring)<-c(paste("b",seq(1,j+1)),"sigma_q") 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(2,3,1,4,5,8) y<-c(3,2,2,7,6,7) n<-length(y) j=1 #Log Likelihood function con argomento «phi» = vettore di parametri ignoti. Log_Lik<-function(phi){ B<-phi[1:(j+1)]; sigma_q<-phi[j+2]; score=-n/2*log(2*pi)-n/2*log(sigma_q)-1/(2*sigma_q)*t(y-X%*%B)%*%(y-X%*%B) (-1)*score # algoritmi di minimizzazione! } phit0<-c(1,1,1) #vettore di stime iniziali # method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent") 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") ############# GLM ############# GLM ############# GLM ############# GLM ############# GLM summary(glm(y~x,family=gaussian(link=identity))) k=1 # variabili x nel modello lineare 1.6720721*n/(n-k-1) # cambiamo il denominatore da n a (n-k-1) Equazione 24 #[1] 2.508108 sqrt(1.6720721*n/(n-k-1)) # errore standard dei residui #[1] 1.583701 -1*solve(H) #           [,1]        [,2]      [,3] #[1,]  1.0755491 -0.20787923 0.0000000       # Varianza distorta di (Intercetta) B_0 #[2,] -0.2078792  0.05422936 0.0000000       # Varianza distorta di B_1 #[3,]  0.0000000  0.00000000 0.9319417 1.0755491*n/(n-k-1) #[1] 1.613324                                # Varianza corretta di (Intercetta) B_0 sqrt(1.613324) #[1] 1.270167                                0.05422936*n/(n-k-1) #[1] 0.08134404                              # Varianza corretta di (Intercetta) B_0 sqrt(0.08134404) #[1] 0.2852088                               # Errore standard di (Intercetta) B_0 ############## Modello più complesso y<-c(3,2,2,7,6,7) # Matrice del modello X <- matrix(c(1,2,2,3, 1,3,2,7, 1,1,1,2, 1,4,5,6, 1,5,9,8, 1,8,2,9), ncol = 4, nrow = 6, byrow = TRUE ) # Stima con funzione glm(): # prendiamo solo le ultime 3 colonne di X summary(glm(y~X[,-1])) ###### Fisher Scoring modificato per renderlo modulabile a qualisiasi dimensione X j=3 # Numero di variabili indipendenti (colonne matrice X – 1) i=0 # Numero progressivo di iterazioni phit0<-matrix(rep(1,j+2)) # Stime iniziali(intercetta + 3 beta_i + sigma_q) n<-length(y) # Grandezza del campione FisherScoring<-c() # data set vuoto per i risultati: repeat{ i=i+1 # prima iterazione (la prima sarà nulla, quindi registreremo i-1) B<-phit0[1:(j+1)] sigma_q<-phit0[j+2] u<-matrix(c( 1/sigma_q*t(X)%*%y-1/sigma_q*t(X)%*%X%*%B, -n/(2*sigma_q)+1/(2*sigma_q^2)*t(y-X%*%B)%*%(y-X%*%B) ),ncol=1) H<-matrix(NA,ncol=(j+2),nrow=(j+2)) H[1:(j+1),1:(j+1)]<--1/sigma_q*t(X)%*%X H[j+2,j+2]<--n/(2*sigma_q^2) H[(j+2),1:(j+1)]<-0 H[1:(j+1),(j+2)]<-0 phit1<- phit0 - solve(H)%*%u #registriamo i risultati nel data set FisherScoring: Result_iter<-c(phit1[1:(j+1)],phit1[j+2]) FisherScoring<-cbind(FisherScoring, Result_iter) #condizione logica per interrompere il ciclo: if(sum(abs(phit1 - phit0) < .0001)==(j+2)){ rownames(FisherScoring)<-c(paste("b",seq(1,j+1)),"sigma_q") colnames(FisherScoring)<-seq(0,(i-1)) print(FisherScoring) break } phit0<-phit1 } FisherScoring k=3 # variabili x nel modello lineare 0.6927083*n/(n-k-1) # cambiamo il denominatore da n a (n-k-1) Equazione 24 #[1] 2.078125 sqrt(2.078125) # errore standard dei residui #[1] 1.44157 vcov.B<-(-1*solve(H))[1:(j+1),1:(j+1)] # minore delle var/covarianze (distorte) dei coefficienti B diag(vcov.B) # diagonale delle varianze (distorte) dei coefficienti B diag(vcov.B)*n/(n-k-1) # varianze corrette per i gdl dei coefficienti B sqrt(diag(vcov.B)*n/(n-k-1)) # errori standard dei coefficienti B