Partea III

Estimatorul de verosimilitate maximă

Slide-uri

Mai jos regăsiți slide-urile pentru partea a 3-a:

Seturi de date

Mai jos regăsiți o parte din seturile de date folosite pentru exemplificare:

Nr.	Setul de date	Descărcare
1	Setul de date Sea
2	Setul de date EColi

Codul R folosit în exemple

Exemplul 1 - Nivelul maxim anual al mării în Veneția

Mai jos regăsiți codul R folosit în exemplul referitor la comparația dintre modelul de regresie clasic și cel în care termenul eroare este repartizat dintr-o repartiție extremă (Gumbel):

# Exemplul - nivelul marii in Venetia
#------------------------------------

sea <- read.csv("Labs 2024-2025/dataIn/sea.csv")

sea <- sea %>%
     rename(x = X, 
            y = x)

sea <- sea %>%
     mutate(an = 1931:1981)


# Modelul clasic de regresie - OLS
# formula: y ~ x1 + x2 + x3 -> y = beta_0 + beta_1 x1 + beta_2 x2 + beta_3 x3

model1 <- lm(y ~ x, 
             data = sea)

coefficients(model1)


# Modelul extremal 

dextval <- function(x, mu, sigma){
     y <- exp(-(x - mu)/sigma)
     
     return(y*exp(-y)/sigma)
}


# logaritmul functiei de verosimilitate
loglik_extval <- function(theta, x, y){
     # theta = (g0, g1, sigma)
     # x - variabila explicativa 
     # y - variabila raspuns 
     
     loglik <- sum(log(dextval(y, theta[1] + theta[2]*x, theta[3])))
     
     return(-loglik)
}


# MLE prin metode numerice
theta_extval <- nlm(loglik_extval, 
                    c(coefficients(model1), sigma(model1)), 
                    x = sea$x, 
                    y = sea$y)

theta_extval$estimate

# ilustrarea grafica a ajustarii modelului pe date
plot(sea$x, 
     sea$y, 
     xlab = "Anul (1 = 1931)", 
     ylab = "Nivelul maxim anul al marii",
     bty = "n")

abline(coefficients(model1), 
       lwd = 2, 
       col = "blue")

abline(theta_extval$estimate[-3], 
       lwd = 2, 
       col = "red")

legend("topleft", 
       legend = c("OLS", "MLE extrem"),
       lwd = c(2, 2),
       col = c("blue", "red"),
       bty = "n")

# graficul cuantila - cuantila

# functia cuantila
qextval <- function(x, mu, sigma){
     mu - sigma*log(-log(x))
} # functia cuantila pentru rep Gumbel


theta_mle <- theta_extval$estimate

res_extval <- sea$y - theta_mle[1] - theta_mle[2]*sea$x # valorile reziduale 

# q-q plot
plot(qextval(ppoints(res_extval), 
             theta_mle[1] + theta_mle[2]*sea$x, theta_mle[3]), 
     sort(res_extval), 
     main = "QQ-plot al reziduurilor pentru modelul extremal",
     cex.main = 0.7, 
     xlab = "Cuantilele modelului extremal",
     ylab = "Valorile reziduale ordonate",
     pch = 16, 
     bty = "n")

Exemplul 2 - Modelul Deming

Mai jos regăsiți codul R folosit în exemplul referitor la modelul de regresie Deming:

# Modelul Deming 
#-----------------

regresie_deming <- function(x, y, var.ratio = 1){
     # x - covariabila (masurata cu eroare)
     # y - variabila raspuns
     # var.ratio - raportul variantelor (y pe x)
     
     # calculul s-urilor 
     
     n <- length(x) # volumul esantionului 
     
     s_xx <- (n - 1) * var(x)
     s_yy <- (n - 1) * var(y)
     s_xy <- (n - 1)  * cov(x, y)
     
     # Estimatorul MLE pentru beta - panta
     
     delta <- (s_yy - var.ratio * s_xx)^2 + 4 * var.ratio * s_xy^2
     
     beta <- (s_yy - var.ratio * s_xx + sqrt(delta)) / (2*s_xy)
     
     # Estimatorul MLE pentru alpha - ordonata la origine
     
     alpha <- mean(y) - beta * mean(x)
     
     # Estimatorul MLE pentru x_i*
     res <- y - alpha - beta * x
     
     x_star <- x + beta * res / (var.ratio + beta^2)
     
     # Estimatorul MLE pentru y_i*
     
     y_star <- y - var.ratio * res / (var.ratio + beta^2)
     
     # Estimatorul pentru sigma_x^2 (folosind la numitor n - 2 nu 2n)
     
     sigma2_x <- (var.ratio * sum((x - x_star)^2) 
                  + sum((y - alpha - beta * x_star)^2)) / (var.ratio * (n-2))
     
     # Estimatorul pentru sigma_y^2 
     
     sigma2_y <- var.ratio * sigma2_x
     
     
     # rezultate
     return(list(x_star = x_star, 
                 y_star = y_star, 
                 intercept = alpha, 
                 slope = beta, 
                 sigma_x = sqrt(sigma2_x), 
                 sigma_y = sqrt(sigma2_y)))
     
     
}


# Exemplu de test 
#.................

set.seed(20230305)

b0 <- 25
b1 <- 1.9

n <- 1000

# valorile adevarate
x_star <- runif(n, 0, 5)
y_star <- b0 + b1 * x_star

sx <- 1 # sigma_x nu sigma_x^2
sy <- 2

vr <- (sy/sx)^2 # cunoscut

# ce observam 
x <- x_star + rnorm(n, 0, sx)
y <- y_star + rnorm(n, 0, sy)


# modelul clasic de regresie 
model_ols <- lm(y ~ x)

coefficients(model_ols)
sigma(model_ols)

# modelul Deming
model_deming <- regresie_deming(x, y, var.ratio = vr)

cbind(model_deming$intercept, 
      model_deming$slope)

cbind(model_deming$sigma_x, 
      model_deming$sigma_y)

# grafic 
plot(x, y, 
     pch = 16, 
     bty = "n", 
     main = "OLS vs Deming")

abline(coefficients(model_ols), 
       lwd = 3, 
       col = "blue")

abline(model_deming$intercept, model_deming$slope, 
       lwd = 3, 
       col = "red")

legend("bottomright", 
       legend = c("OLS", "Deming"), 
       col = c("blue", "red"), 
       lwd = 3, bty = "n")

Exemplul 3 - Lanț Markov

Mai jos regăsiți codul R folosit în exemplul referitor la modelul Markov:

• exemplul 1 - model Markov general

# Exemplul 1
#-------------

# Construim o functie care sa estimeze matricea de tranzitie
estimare_matrice_tranzitie <- function(secventa){
    # valorile unice ale sirului 
    stari <- unique(secventa)
    nstari <- length(stari)
    
    # Initializam matricea de frecventa
    matrice_frecventa <- matrix(0, nrow = nstari, ncol = nstari)
    rownames(matrice_frecventa) <- stari
    colnames(matrice_frecventa) <- stari
    
    # Numar de tranzitii
    for (t in 1:(length(secventa) - 1)){
        i <- which(stari == secventa[t])
        j <- which(stari == secventa[t + 1])
        
        matrice_frecventa[i, j] <- matrice_frecventa[i, j] + 1 #N_{i,j}
    }
    
    # Matricea de tranzitie
    matrice_tranzitie <- matrix(0, nrow = nstari, ncol = nstari)
    rownames(matrice_tranzitie) <- stari
    colnames(matrice_tranzitie) <- stari
    
    for(i in 1:nstari){
        suma_linie <- sum(matrice_frecventa[i, ]) # N_{i.}
        if (suma_linie > 0){
            matrice_tranzitie[i, ] <- matrice_frecventa[i, ] / suma_linie
        }
    }
    
    return(list(matrice_frecventa = matrice_frecventa, 
                matrice_tranzitie = matrice_tranzitie))
    
}


secventa <- sample(c("a", "b", "c"), 1000, replace = TRUE)

estimare_matrice_tranzitie(secventa)


# Consideram un model pentru vreme: Soare (S), Innorat (N), Ploaie (P)
set.seed(1234)

stari_vreme <- c("S", "N", "P")

# adevarata matrice de tranzitie 

matrice_tranzitie_adev <- matrix(c(
    0.65, 0.25, 0.1, # 
    0.23, 0.57, 0.2,
    0.25, 0.42, 0.33
), nrow = 3, byrow = TRUE)

rownames(matrice_tranzitie_adev) <- stari_vreme
colnames(matrice_tranzitie_adev) <- stari_vreme

# functie de generare Markov

simultate_Markov <- function(starea_init, mat_tranzitie, n){
    # starea_init - starea din care pleaca lantul
    # mat_tranzitie - matricea de trecere
    # n - marimea lantului
    
    stari <- rownames(mat_tranzitie)
    lant <- character(n + 1)
    
    lant[1] <- starea_init
    
    for(t in 1:n){
        # care este starea curenta
        starea_curenta <- which(stari == lant[t])
        prob_stare_curenta <- mat_tranzitie[starea_curenta, ]
        lant[t + 1] <- sample(stari, 1, prob = prob_stare_curenta)
    }
    
    return(lant)
}

simultate_Markov("S", matrice_tranzitie_adev, 10)

secventa_vreme <- simultate_Markov("S", matrice_tranzitie_adev, 10000)

# matricea de tranzitie estimata

mat <- estimare_matrice_tranzitie(secventa_vreme)

# calculam logaritmul functiei de verosimilitate
calcul_loglik_markov <- function(secventa, mat_tranzitie){
    stari <- rownames(mat_tranzitie)
    
    loglik <- 0
    
    for(t in 1:(length(secventa) - 1)){
        i <- which(stari == secventa[t])
        j <- which(stari == secventa[t + 1])
        
        loglik <- loglik + log(mat_tranzitie[i, j])
    }
    
    return(loglik)
}

# verificare
calcul_loglik_markov(secventa_vreme, matrice_tranzitie_adev)
calcul_loglik_markov(secventa_vreme, mat$matrice_tranzitie)

• exemplul 2 - modelul Markov pentru secvența ADN

# Exemplul 2 - ADN
#--------------------

library(seqinr)

ecoli <- read.fasta("Labs 2024-2025/dataIn/Ecoli.fasta")
str(ecoli)

ecoli <- ecoli[[1]]
ecoli[1:10]

length(ecoli)

tstart <- proc.time()
estimare_matrice_tranzitie(ecoli)
proc.time() - tstart


# o varianta optimizata

cuvinte <- function(marime){
    nucleo <- c("a", "c", "g", "t")
    
    if (marime == 1){
        return(nucleo)
    }else{
        s <- outer(nucleo, cuvinte(marime - 1), paste, sep = "")
        return(as.vector(s))
    }
}

cuvinte(3)

count_nucleo <- function(secventa, marime){
    nucleo <- c("a", "c", "g", "t")
    
    ind <- seq(from = 1, to = length(secventa), by = 1)
    s <- secventa
    
    s.levels <- levels(as.factor(cuvinte(marime)))
    
    if (marime >= 2){
        for (i in 2:marime){
            s <- paste(s, secventa[ind + i - 1], sep = "")
        }
    }
    
    count_frec <- table(factor(s, levels = s.levels))
    
    return(count_frec)
}

estimare_matrice_tranzitieM1 <- function(secventa){
    # valorile unice ale sirului 
    stari <- c("a", "c", "g", "t")
    nstari <- 4
    
    # Numar de tranzitii
    dinucleo <- count_nucleo(secventa, 2)
    
    # Mtricea de frecventa
    matrice_frecventa <- matrix(dinucleo, nrow = nstari, ncol = nstari, byrow = TRUE)
    rownames(matrice_frecventa) <- stari
    colnames(matrice_frecventa) <- stari

    
    # Matricea de tranzitie
    matrice_tranzitie <- matrix(0, nrow = nstari, ncol = nstari)
    rownames(matrice_tranzitie) <- stari
    colnames(matrice_tranzitie) <- stari
    
    for(i in 1:nstari){
        suma_linie <- sum(matrice_frecventa[i, ]) # N_{i.}
        if (suma_linie > 0){
            matrice_tranzitie[i, ] <- matrice_frecventa[i, ] / suma_linie
        }
    }
    
    return(list(matrice_frecventa = matrice_frecventa, 
                matrice_tranzitie = matrice_tranzitie))
    
}


tstart <- proc.time()
estimare_matrice_tranzitieM1(ecoli)
proc.time() - tstart