Proprietăți ale estimatorilor

Exemplu de comparare a trei estimatori

Exercițiul 1 Fie $X_1,X_2,\ldots,X_n$ un eșantion de volum $n$ dintr-o populație normală de medie $\mu$ și varianță $\sigma^2$. Atunci

\[ \hat{\mu}_1 = \frac{1}{n}\sum_{i=1}^{n}X_i, \quad \hat{\mu}_2 = M_n\,(\text{mediana}\,), \quad \hat{\mu}_3 = \frac{X_{(1)} + X_{(n)}}{2} \]

sunt trei estimatori punctuali pentru $\mu$. Creați o funcție care să ilustreze cum sunt repartizați cei trei estimatori. Începeți cu $n = 10$, $\mu = 0$ și $\sigma^2 = 1$ și trasați histogramele pentru a-i compara. Ce se întâmplă dacă schimbați $n$, $\mu$ sau $\sigma^2$ ?

Soluție. Vom crea o funcție numită comparare_estimatori care va construi repartițiile celor trei estimatori:

comparare_estimatori <- function(n, mu, sigma, S){
  # Initializam
  mu1 <- numeric(S)
  mu2 <- numeric(S)
  mu3 <- numeric(S)
  
  # repetam experimentul de S ori
  for (i in 1:S){
    x <- rnorm(n, mean = mu, sd = sigma)
    
    # calculam estimatorii
    mu1[i] <- mean(x)
    mu2[i] <- median(x)
    mu3[i] <- (min(x)+max(x))/2
  }
  
 return(list(M = cbind(mu1 = mu1, mu2 = mu2, mu3 = mu3), 
             vM = cbind(var_mu1 = var(mu1), var_mu2 = var(mu2), var_mu3 = var(mu3))))
  
}

Ca exemplu să considerăm $\mu = 3$ și $\sigma^2 = 2$:

# Exemplu
mu <- 3
sigma <- 2

n <- 1000
N <- 10000

mu_hat <- comparare_estimatori(n, mu, sigma, N)

# variantele celor trei estimatori
mu_hat$vM

        var_mu1    var_mu2   var_mu3
[1,] 0.00404305 0.00646494 0.2425589

# varianta teoretica a lui \bar{X}_n
sigma^2/n

[1] 0.004

Pentru a ilustra grafic histogramele celor trei estimatori, considerăm $\mu = 0$ și $\sigma^2 = 1$ și avem:

mu <- 0
sigma <- 1

n <- 1000
S <- 10000

a <- comparare_estimatori(n, mu, sigma, S)

par(mfrow = c(1,3))
hist(a$M[,1], freq=FALSE, xlab=expression(hat(mu)[1]), 
     col="gray80", border="white", 
     main = expression(paste("Histograma lui ", hat(mu)[1])),
     ylab = "Densitatea")
abline(v=mu, col = "brown3", lty = 2)

hist(a$M[,2], freq=FALSE, xlab=expression(hat(mu)[2]), 
     col="gray80", border="white",
     main = expression(paste("Histograma lui ", hat(mu)[2])),
     ylab = "Densitatea")
abline(v=mu, col = "brown3", lty = 2)

hist(a$M[,3], freq=FALSE, xlab=expression(hat(mu)[3]), 
     col="gray80", border="white",
     main = expression(paste("Histograma lui ", hat(mu)[3])),
     ylab = "Densitatea")
abline(v=mu, col = "brown3", lty = 2)

Figura 1: Repartițile de eșantionare ale celor trei estimatori

Ilustrarea consistenței unui estimator

Exercițiul 2 Fie $X_1,X_2,\ldots,X_n$ un eșantion de volum $n$ dintr-o populație $Pois(\theta)$. Ilustrați grafic consistența estimatorului $\hat{\theta}_n = S_n^2$ trasând histograma repartiției lui $\hat{\theta}_n$ pentru $n\in\{10,25,50,100\}$. Ce observați?

Soluție. Considerăm funcția pois_est care pentru $\theta$ fixat simulează repartiția estimatorului $\hat{\theta}_n$:

pois_est1 <- function(n, theta, S){
  # initializare
  sigma1 <- numeric(S)
  
  for (i in 1:S){
    x <- rpois(n, theta)
    sigma1[i] <- var(x)
  }
  # afisam varianta estimatorului
  print(paste0("Pentru n = ", n," varianta estimatorului este ", var(sigma1)))
  return(sigma1)
}

Considerând $\theta = 3$ și $n\in\{10,25,50,100\}$ avem:

theta <- 3

par(mfrow=c(2,2))

vals <- c(10, 25, 50, 100)

for (i in vals){
     a <- pois_est1(i, theta, 50000)
     hist(a, freq=FALSE, 
          xlab=expression(hat(theta)[n]), 
          col="gray80", border="white", 
          main = paste0("n = ", i), xlim = c(0,12),
          ylab = "Densitatea")
     abline(v=theta, col = "brown3", lty = 2)
}

[1] "Pentru n = 10 varianta estimatorului este 2.2957469740778"

[1] "Pentru n = 25 varianta estimatorului este 0.870294631741341"

[1] "Pentru n = 50 varianta estimatorului este 0.420902500419613"

[1] "Pentru n = 100 varianta estimatorului este 0.213758498328129"

Figura 2: Ilustrarea unor densități cu coeficientul de asimetrie negativ, respectiv pozitiv.

Exercițiul 3 Ce se întâmplă dacă în loc de $\hat{\theta}_n$ considerăm estimatorul $\tilde{\theta}_n = \bar{X}_n$ sau estimatorul $\dot{\theta}_n = \sqrt{\bar{X}_n S_n^2}$ ?

Soluție. Considerăm acum funcția est_Poisson care va simula repartițiile celor trei estimatori și varianțele acestora:

est_Poisson <- function(n = 10, theta = 1, N = 1000){
  e <- matrix(0, nrow = N, ncol = 3)
  
  for (i in 1:N){
    x <- rpois(n, theta)
    
    e[i, 1] <- mean(x)
    e[i, 2] <- var(x)
    e[i, 3] <- sqrt(e[i, 1] * e[i, 2])
  }
  
  return(list(e_Pois = e, 
              v_Pois = c(var(e[, 1]), var(e[, 2]), var(e[, 3]))))
}

Ilustrăm grafic repartițiile de eșantionare pentru cei trei estimatori considerând $\theta = 3$ și $n\in\{10,25,50,100\}$ și observăm cum acestea se concentrează în jurul parametrului $\theta$ pe care îl estimează:

theta <- 3
n <- c(10, 25, 50, 100)

N <- 10000

par(mfrow = c(3, length(n)))

par(mai=c(0.5,0.5,0.5,0.5), bty = "n")

# par(fig = c(0, 1, 0, 0.45))

titles <- c("hat(theta)[n]", 
            "tilde(theta)[n]", 
            "dot(theta)[n]")

for (i in 1:3){
  for (j in 1:length(n)){
    e_hat <- est_Poisson(n[j], theta, N)$e_Pois
    
    hist(e_hat[, i], probability = TRUE, 
         main = paste("n =", n[j]),
         xlab = str2lang(titles[i]),
         ylab = "Densitate", 
         xlim = c(0, 12))
    
    abline(v = theta, 
           col = "red", 
           lwd = 2, 
           lty = 2)
  }
}

Figura 3: Repartițile de eșantionare ale celor trei estimatori pentru diverse valori ale lui n

Observăm în figura de mai jos convergența varianțelor celor trei estimatori către $0$:

# variantele se duc la 0 
n <- seq(10, 1000, 10)

ve <- matrix(0, nrow = length(n), ncol = 3)

for(i in 1:length(n)){
  ve[i, ] <- est_Poisson(n[i], theta, N)$v_Pois
}

matplot(ve, 
        type = "l", 
        lwd = 2, 
        col = c(myred, myblue, mygreen), 
        bty = "n")

legend("topright", 
       legend = c(expression(hat(theta)[n]), expression(tilde(theta)[n]), expression(dot(theta)[n])), 
       col = c(myred, myblue, mygreen),
       lty = 1,
       lwd = 2,
       bty = "n", 
       seg.len = 1.5)

Figura 4: Ilustrarea convergenței varianței celor trei estimatori la 0.

Estimare prin metoda verosimilității maxime

Exemplu: EVM nu este întotdeauna media eșantionului chiar dacă $\mathbb{E}_{\theta}[\hat{\theta}_n] = \theta$

Exercițiul 4 Fie $X_1,X_2,\ldots,X_n$ un eșantion de talie $n$ dintr-o populație Laplace $L(\theta, c)$ a cărei densitate este dată de formula

\[ f_{\theta, c}(x) = \frac{1}{2c}e^{-\frac{|x-\theta|}{c}}, \quad -\infty<x<\infty \]

Ilustrați grafic densitatea și funcția de repartiție a repartiției Laplace pentru diferite valori ale parametrilor $\theta$ (de locație) și $c$ (de scală), e.g. $\theta\in\{0, 3\}$ și $c\in\{1,2,3,4\}$.
Determinați estimatorul de verosimilitate maximă $\hat{\theta}_n$ pentru $\theta$.

Soluție. Avem:

Se poate arăta cu ușurință (am văzut în ?@exr-stat-sim-minv-cont-e7) că funcția de repartiție a repartiției $\mathrm{Laplace}(\theta, c)$ este

\[ F_{\theta, c}(x) = \frac{1}{2} + \frac{1}{2}\operatorname{sgn}(x-\theta)\left(1-e^{-\frac{|x-\theta|}{c}}\right) = \left\{\begin{array}{ll} \frac{1}{2}e^{-\frac{|x-\theta|}{c}}, & x<\theta\\ 1-\frac{1}{2}e^{-\frac{|x-\theta|}{c}}, & x\geq\theta \end{array}\right. \]

În Figura 5 ilustrăm grafic pentru o serie de parametrii $(\theta,c)$ densitatea și funcția de repartiție pentru repartiția $\mathrm{Laplace}(\theta, c)$:

Figura 5: Ilustrarea densității și a funcției de repartiție pentru repartiția Laplace.

Pentru a determina estimatorul de verosimilitate maximă să observăm că funcția de verosimilitate este

\[ L_n(\theta|\mathbf{X}) = \prod_{i=1}^{n}\left(\frac{1}{2c}e^{-\frac{|X_i-\theta|}{c}}\right) = \frac{1}{(2c)^n}e^{-\sum_{i=1}^{n}\frac{|X_i-\theta|}{c}} \]

și acesta ia valoarea maximă pentru toate valorile lui $\theta$ care minimizează funcția de la exponent

\[ M(\theta) = \sum_{i=1}^{n}|X_i-\theta| = \sum_{i=1}^{n}|X_{(i)}-\theta|, \]

unde $x_{(i)}$ este statistica de ordine de rang $i$. Se poate vedea că funcția $M(\theta)$ este continuă și afină pe porțiuni din Figura 6 de mai jos (pentru un eșantion de volum $10$ dintr-o populație $\mathrm{Laplace}(3,1)$ - creați o funcție care vă permite să generați observații repartizate Laplace).

rLaplace <- function (n, mu = 0, b = 1) 
{
    u <- runif(n) - 0.5
    x <- mu - b * sign(u) * log(1 - 2 * abs(u))
    return(x)
}

theta0 <- 3
b <- 1

set.seed(333)
x <- rLaplace(10, mu = theta0, b = b)

M_theta <- function(x, theta){
  sapply(theta, function(t){sum(abs(x-t))})
}

theta <- seq(min(x), max(x), length.out = 500)
M <- M_theta(x, theta)

plot(theta, M, type = "l", 
     xlab = TeX("$\\theta$"),
     ylab = TeX("$M(\\theta)$"),
     bty = "n", lwd = 2,
     col = myblue)
points(x, M_theta(x, x), 
       pch = 16, 
       cex = 1.2, 
       col = myred)

Figura 6: Graficul funcției $M(\theta)$ pentru un eșantion de volum 10.

Observăm că dacă $\theta$ se află între statistica de ordine de rang $m$ și cea de rang $m+1$, i.e. $X_{(m)}\leq \theta\leq X_{(m+1)}$, atunci am avea că $X_{(i)} \leq X_{(m)} \leq \theta$ dacă $i\leq m$ și $\theta\leq X_{(m+1)}\leq X_{(i)}$ dacă $m+1\leq i\leq n$, prin urmare

\[ M(\theta) = \sum_{i=1}^{n}|X_{(i)}-\theta| = \sum_{i=1}^{m}(\theta - X_{(i)}) + \sum_{i=m+1}^{n}(X_{(i)}-\theta) \]

deci dacă $X_{(m)}< \theta< X_{(m+1)}$ atunci

\[ \frac{d}{d\theta}M(\theta) = m - (n-m) = 2m-n. \]

Astfel, $M'(\theta)<0$ (și $M(\theta)$ este descrescătoare) dacă $m<\frac{n}{2}$ și $M'(\theta)>0$ (și $M(\theta)$ este crescătoare) dacă $m>\frac{n}{2}$. Dacă $n = 2k+1$ este impar, atunci $\frac{n}{2} = k +\frac{1}{2}$ iar $M(\theta)$ este strict descrescătoare dacă $\theta<X_{(k+1)}$ și strict crescătoare dacă $\theta>X_{(k+1)}$ de unde deducem că minimul se atinge pentru $\theta = X_{(k+1)}$.

Dacă $n = 2k$ este par atunci, raționând asemănător, deducem că $M(\theta)$ este minimizată pentru orice punct din intervalul $(X_{(k)}, X_{(k+1)})$, deci orice punct din acest interval va maximiza și funcția de verosimilitate. Prin convenție alegem estimatorul de verosimilitate maximă să fie mijlocul acestui interval, i.e. $\theta = \frac{X_{(k)} + X_{(k+1)}}{2}$.

Prin urmare am găsit că estimatorul de verosimilitate maximă este mediana eșantionului

\[ \hat{\theta}_n = \left\{\begin{array}{ll} X_{\left(\frac{n+1}{2}\right)}, & \text{$n$ impar}\\ \frac{X_{\left(\frac{n}{2}\right)} + X_{\left(\frac{n}{2}+1\right)}}{2}, & \text{$n$ par}\\ \end{array}\right. \]

În Figura 7 avem ilustrat logaritmul funcției de verosimilitate pentru un eșantion de volum par (stânga) și unul de volum impar (dreapta):

theta0 <- 0
b <- 1

logLikelihoodLaplace <- function(x, theta, b){
  sapply(theta, function(t){
    -length(x)*log(2*b)-sum(abs(x-t))})
}

par(mfrow = c(1,2))

# Esantion par
n <- 10

set.seed(333)
x <- rLaplace(n, mu = theta0, b = b)

theta <- seq(-1,1, length.out = 100)

y <- logLikelihoodLaplace(x, theta, b)

plot(theta, y, type = "l",
     bty = "n", lwd = 2,
     col = myblue,
     xlab = TeX("$\\theta$"),
     ylab = TeX("$l(\\theta | x)$"), 
     main = paste0("Esantion par\n n = ", n), 
     cex.main = 0.8)
points(x, logLikelihoodLaplace(x, x, b),
       pch = 16, 
       col = myred, 
       cex = 0.7)

# Esantion impar
n <- 13

set.seed(1234)
x <- rLaplace(n, mu = theta0, b = b)

theta <- seq(-1,1, length.out = 100)

y <- logLikelihoodLaplace(x, theta, b)

plot(theta, y, type = "l",
     bty = "n", lwd = 2,
     col = myblue,
     xlab = TeX("$\\theta$"),
     ylab = TeX("$l(\\theta | x)$"), 
     main = paste0("Esantion impar\n n = ", n), 
     cex.main = 0.8)
points(x, logLikelihoodLaplace(x, x, b),
       pch = 16, 
       col = myred, 
       cex = 0.7)

Figura 7: Logaritmul funcției de verosimilitate pentru un eșantion de volum 10 (stânga) și respectiv 13 (dreapta).

Exemplu de EVM determinat prin soluții numerice

Exercițiul 5 Fie $X_1,X_2,\ldots,X_n$ un eșantion de talie $n$ dintr-o populație logistică a cărei densitate este dată de formula

\[ f_{\theta}(x) = \frac{e^{-(x-\theta)}}{\left(1+e^{-(x-\theta)}\right)^2}, \quad x\in\mathbb{R},\, \theta\in\mathbb{R} \]

Determinați estimatorul de verosimilitate maximă $\hat{\theta}_n$ pentru $\theta$.

Soluție. Densitatea de repartiție și funcția de repartiție a repartiției logistice sunt ilustrate mai jos (în R se folosesc funcțiile: rlogis, dlogis, plogis și respectiv qlogis):

Figura 8: Ilustrarea densității și a funcției de repartiție pentru repartiția logistică.

Observăm că funcția de verosimilitate este dată de

\[ L_n(\theta|\mathbf{x}) = \prod_{i=1}^{n}f_{\theta}(x_i) = \prod_{i=1}^{n}\frac{e^{-(x_i-\theta)}}{\left(1+e^{-(x_i-\theta)}\right)^2} \]

iar logaritmul funcției de verosimilitate este

\[ l_n(\theta|\mathbf{x}) = \sum_{i=1}^{n}\log{f_{\theta}(x_i)} = n\theta - n\bar{x}_n - 2\sum_{i=1}^{n}\log{\left(1+e^{-(x_i-\theta)}\right)}. \]

Pentru a găsi valoarea lui $\theta$ care maximizează logaritmul funcției de verosimilitate și prin urmare a funcției de verosimilitate trebuie să rezolvăm ecuația $l'(\theta|\mathbf{x}) = 0$, unde derivata lui $l(\theta|\mathbf{x})$ este

\[ l_n'(\theta|\mathbf{x}) = n - 2\sum_{i = 1}^{n}\frac{e^{-(x_i-\theta)}}{1+e^{-(x_i-\theta)}} \]

ceea ce conduce la ecuația

\[ \sum_{i = 1}^{n}\frac{e^{-(x_i-\theta)}}{1+e^{-(x_i-\theta)}} = \frac{n}{2} \tag{$\star$} \]

Chiar dacă această ecuație nu se simplifică, se poate arăta că această ecuația admite soluție unică. Observăm că derivata parțiala a membrului drept în ($\star$) devine

\[ \frac{\partial }{\partial \theta}\sum_{i = 1}^{n}\frac{e^{-(x_i-\theta)}}{1+e^{-(x_i-\theta)}} = \sum_{i = 1}^{n}\frac{e^{-(x_i-\theta)}}{\left(1+e^{-(x_i-\theta)}\right)^2}>0 \]

ceea ce arată că membrul stâng este o funcție strict crescătoare în $\theta$. Cum membrul stâng în ($\star$) tinde spre $0$ atunci când $\theta\to-\infty$ și spre $n$ pentru $\theta\to\infty$ deducem că ecuația ($\star$) admite soluție unică (vezi graficul de mai jos).

Figura 9: Ilustrarea ecuație de verosimilitate în cazul repartiției logistice.

Cum nu putem găsi o soluție a ecuației $l_n'(\theta|\mathbb{x}) = 0$ sub formă compactă, este necesar să apelăm la metode numerice. O astfel de metodă numerică este binecunoscuta metodă a lui Newton-Raphson. Metoda presupune să începem cu o valoare (soluție) inițială $\hat{\theta}^{(0)}$ și să alegem, plecând de la aceasta, o nouă valoare $\hat{\theta}^{(1)}$ definită prin

\[ \hat{\theta}^{(1)} = \hat{\theta}^{(0)} - \frac{l_n'\left(\hat{\theta}^{(0)}\right)}{l_n''\left(\hat{\theta}^{(0)}\right)}, \]

adică $\hat{\theta}^{(1)}$ este intersecția cu axa absciselor a tangentei în punctul $\left(\hat{\theta}^{(0)}, l_n'\left(\hat{\theta}^{(0)}\right)\right)$ la graficul funcției $l_n'(\theta)$. Ideea este de a itera procesul până când soluția converge, cu alte cuvinte pornind de la o valoare rezonabilă de start $\hat{\theta}^{(0)}$ la pasul $k+1$ avem

\[ \hat{\theta}^{(k+1)} = \hat{\theta}^{(k)} - \frac{l_n'\left(\hat{\theta}^{(k)}\right)}{l_n''\left(\hat{\theta}^{(k)}\right)} \]

și oprim procesul atunci când $k$ este suficient de mare și/sau $\left|\hat{\theta}^{(k+1)} - \hat{\theta}^{(k)}\right|$ este suficient de mic. Următorul grafic ilustrează grafic algoritmul lui Newton:

Figura 10: Ilustrarea metodei lui Newton pentru rezolvarea ecuației $l_n'(\theta|\mathbb{x}) = 0$

Obs: Singurul lucru care se schimbă atunci când trecem de la scalar la vector, este funcția $l_n(\theta)$ care acum este o funcție de $p>1$ variabile, $\theta = (\theta_1, \theta_2, \ldots, \theta_p)^{\intercal}\in\mathbb{R}^p$. În acest context $l'(\theta)$ este un vector de derivate parțiale iar $l''(\theta)$ este o matrice de derivate parțiale de ordin doi. Prin urmare iterațiile din metoda lui Newton sunt

\[ \hat{\theta}^{(k+1)} = \hat{\theta}^{(k)} - \left[l_n''\left(\hat{\theta}^{(k)}\right)\right]^{-1}l_n'\left(\hat{\theta}^{(k)}\right) \]

unde $[\cdot]^{-1}$ este pseudoinversa unei matrice.

Funcția de mai jos implementează metoda lui Newton pentru cazul multidimensional:

# Metoda lui Newton

newton <- function(f, df, x0, eps=1e-08, maxiter=1000, ...) {
  # in caz ca nu e incarcat pachetul sa putem accesa pseudoinversa
  if(!exists("ginv")) library(MASS) 
  
  x <- x0
  k <- 0
  
  repeat {
    k <- k + 1
    
    x.new <- x - as.numeric(ginv(df(x, ...)) %*% f(x, ...))
    
    if(mean(abs(x.new - x)) < eps | k >= maxiter) {
      if(k >= maxiter) warning("S-a atins numarul maxim de iteratii!")
      break
    }
    x <- x.new
  }
  out <- list(solution = x.new, value = f(x.new, ...), iter = k)
  
  return(out)
}

Să presupunem că am observat următorul eșantion de volum $20$ din repartiția logistică:

 [1]  6.996304  9.970107 12.304991 11.259549  6.326912  5.378941  4.299639
 [8]  8.484635  5.601117  7.094335  6.324731  6.868456  9.753360  8.042095
[15]  8.227830 10.977982  7.743096  7.722159  8.562884  6.968356

set.seed(112)
x <- rlogis(20, location = 7.5)

n <- length(x)
dl <- function(theta) n - 2 * sum(exp(theta - x) / (1 + exp(theta - x)))
ddl <- function(theta) {-2 * sum(exp(theta - x) / (1 + exp(theta - x))^2)}

logis.newton <- newton(dl, ddl, median(x))

și aplicând metoda lui Newton găsim estimatorul de verosimilitate maximă $\hat{\theta}_n=$ 7.7933 după numai 3 iterații (datele au fost simulate folosind $\theta = 7.5$).

Elemente de estimare punctuală în `R`

Proprietăți ale estimatorilor

Exemplu de comparare a trei estimatori

Ilustrarea consistenței unui estimator

Estimare prin metoda verosimilității maxime

Exemplu: EVM nu este întotdeauna media eșantionului chiar dacă \(\mathbb{E}_{\theta}[\hat{\theta}_n] = \theta\)

Exemplu de EVM determinat prin soluții numerice