Partea I

Principiul estimatorului de verosimilitate maximă

Slide-uri

Mai jos regăsiți slide-urile pentru prima parte:

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 ZIP

2 Setul de date Hurricane

Codul R folosit în exemple

Exemplul 1 - Model ZIP

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

Codul R 
x <- c(rep(0, 182),
       rep(1, 41),
       rep(2, 12),
       rep(3, 2),
       rep(4, 2),
       rep(7, 1))

table(x)

# daca pp ca suntem in cadrul unui model Poisson
theta_hat <- mean(x)

expected_counts <- length(x)*dpois(0:7, lambda = theta_hat)
round(expected_counts, 2)

obs_counts <- c(182, 41, 12, 2, 2, 0, 0, 1)

barplot(rbind(obs_counts,
              expected_counts),
        names.arg = 0:7,
        beside = TRUE)

# logaritmul fctei de verosimilitate
loglik_zip <- function(theta, x){
     # theta - vectorul de parametrii (prob p si media Poisson theta de la curs)
     p <- theta[1]
     t <- theta[2]

     n <- length(x)
     n0 <- sum(x == 0)

     loglik <- n0*log(p + (1-p)*dpois(0, t)) + (n - n0)*(log(1-p) - t) +
          log(t)*sum(x[x > 0]) - log(prod(factorial(x[x>0])))

     return(-loglik)
}

?nlm

opt_zip <- nlm(loglik_zip, c(0.5, 0.5), x = x)
theta_hat_zip <- opt_zip$estimate

u <- 0:7
expected_zip <- length(x)*ifelse(u == 0,
                                 theta_hat_zip[1] + (1 - theta_hat_zip[1])*dpois(u, theta_hat_zip[2]),
                                 (1 - theta_hat_zip[1])*dpois(u, theta_hat_zip[2]))
round(expected_zip, 2)

barplot(rbind(obs_counts,
              expected_zip),
        names.arg = 0:7,
        beside = TRUE)

Exemplul 2 - Metoda capturării-recapturării

Mai jos regăsiți codul R folosit în exemplul referitor la metoda capturării-recapturării:

Codul R 
# Exemplul cu porci mistreti
M <- 25
n <- 25
k <- 7

lik_hyper <- function(N, M, n, k){
     out <- dhyper(k, M, N - M, n)

     return(out)
}

# graficul functiei de verosimilitate
N <- max(n, M + n - k) : (max(n, M + n - k) + 200)

y <- lik_hyper(N, M, n, k)

N_hat <- if(M*n/k != floor(M*n/k)){
     floor(M*n/k)
}else{
     c(M*n/k - 1, M*n/k)
}

plot(N, y,
     type = "p",
     pch = 16,
     bty = "n",
     xlim = c(N_hat[1] - 20, N_hat[1] + 20))


points(N_hat, lik_hyper(N_hat, M, n, k),
       col = "red",
       pch = 16,
       cex = 1.5)

y[which.max(y)]

n_mle <- max(n, M + n - k) + which.max(y) - 1

lik_hyper(n_mle, M, n, k)
lik_hyper(N_hat, M, n, k)

Exemplul 3 - Model Gamma

Mai jos regăsiți codul R folosit în exemplul referitor la modelul Gamma folosind setul de date Hurricane:

Codul R 
# Exemplul cu uragane
data_hurricane <- read.csv("dataIn/hurrprec.csv")
hurricane <- data_hurricane$x

hist(hurricane,
     probability = TRUE)

# functia de verosimilitate
n <- length(hurricane)

# sum(log(x_i))
S1 <- sum(log(hurricane))

# sum(x_i)
S2 <- sum(hurricane)


loglik_gamma <- function(alpha, beta){
     -n*log(gamma(alpha)) - n*alpha*log(beta) + (alpha - 1)*S1 - S2/beta
}

# grafic 3d
?persp

alpha <- seq(1, 3.5, length.out = 100)
beta <- seq(2, 6, length.out = 100)

?outer
z <- outer(alpha, beta, loglik_gamma)

persp(alpha, beta, z)

contour(alpha,
        beta,
        z,
        nlevels = 30,
        drawlabels = FALSE)


# install.packages("rgl")
library("rgl")

persp3d(alpha, beta, z,
        col = "royalblue")

# MLE

nloglik_gamma <- function(theta){
     - loglik_gamma(alpha = theta[1],
                    beta = theta[2])
}

theta0 <- c(mean(hurricane)^2/var(hurricane),
            var(hurricane)/mean(hurricane))

theta_hat <- nlm(nloglik_gamma, theta0)$estimate

hist(hurricane,
     probability = TRUE)

curve(dgamma(x, shape = theta_hat[1], scale = theta_hat[2]),
      col = "royalblue",
      lwd = 2,
      add = TRUE)