Verosimiglianza

Il grafico della funzione di verosimiglianza si può ottenere con curve. Per esempio la verosimiglianza, nel caso di un campione da una Bernoulli si rappresenta con le istruzioni seguenti:

likber = function(p) { 
           n = length(y) 
           r = sum(y)
           p^r*(1-p)^(n-r)
         }
y = c(1,0,0,1,0,0,0,0,0,0)
curve(likber(x), from = 0, to = 1, xlab = 'p', ylab = 'Verosimiglianza', lwd = 2, col = "red")

Notate che la verosimiglianza viene calcolata usando il vettore dei dati y definito fuori della funzione. La log-verosimiglianza si ottiene facilmente con

curve(log(likber(x)), from = 0, to = 1, xlab = 'p', ylab = 'Log verosimiglianza', col = "blue", lwd = 2)

Confronto di due verosimiglianze normalizzate

Definiamo la verosimiglianza normalizzata \(L(p)/L(\hat p)\) per il modello Bernoulliano. Questa volta i dati sono rappresentati dal numero totale di successi r e dal numero di prove n.

lik.rel = function(p) {
   ph = r/n
   p^r*(1-p)^(n-r) /( ph^r*(1-ph)^(n-r))
}

Quindi confrontiamo tre casi in cui la stima di massima verosimiglianza è sempre \(\hat p = 0.2\), ma la proporzione di successi è rispettivamente \(2/10\), \(8/40\), \(20/100\). Rappresentiamo le verosimiglianze relative sullo stesso grafico.

n = 10; r = 2
curve(lik.rel(x), from = 0, to = 1, xlab = 'p', ylab = 'Verosimiglianza', col = "red", lwd = 2)
n = 40; r = 8
curve(lik.rel(x), add = TRUE, col = "orange", lwd = 2)
n = 100; r = 20
curve(lik.rel(x), add = TRUE, col = "purple", lwd = 2)
legend("topright", c("n = 10", "n = 40", "n=100"), col = c("red", "orange", "purple"), lty = c(1,1,1), lwd = c(2,2,2))

Come si vede la funzione di verosimiglianza diventa sempre più concentrata sulla stima perché l’informzione contenuta nel campione aumenta.

Esempio: modello cattura-ricattura

Supponiamo di voler stimare la dimensione \(N\) di una popolazione di tassi in una zona. I tassi aumentano la diffusione della tubercolosi bovina e pretanto è importante avere un’idea della dimensione della popolazione. La tecnica per stimare \(N\) è catturare un campione di tassi, marcarli e poi rilasciarli. Quindi si estrae un campione casuale semplice senza ripetizione di tassi e si contano quanti di questi sono marcati. Supponiamo che i tassi marcati siano \(S = 25\). Supponiamo poi di estrarre un campione di 60 tassi e di trovarne 5 marcati. Qual è la probabilità di trovare 5 successi su 60 estrazioni senza ripetizione? Si usa l’ipergeometrica e si trova \[ P(5 \text{ successi} ) = \frac{{25 \choose 5}{N-25 \choose 55}}{{N \choose 60}} \] Questa probabilità è una funzione di \(N\) ed è la verosimiglianza \(L(N)\). Se la calcoliamo per vari valori si \(N\) (diciamo tra 150 e 1110) e ne facciamo il grafico otteniamo

N = 150:1100
L = choose(25,5)*choose(N-25, 55)/choose(N, 60)
plot(N, L, type = "l", xlab = "N", ylab = "L(N)", col = "purple", lwd = 2)
abline(v = 299, lty = 2)

Il massimo è in corrispondenza di \(N = 299\):

N[which(L == max(L))]
[1] 299

Approssimazione al secondo ordine della verosimiglianza

La log-verosimiglianza normalizzata è \(l(\theta) - l(\hat \theta)\). Sviluppando in serie di Taylor al secondo ordine \[ l(\theta) \simeq l(\hat \theta) + l'(\hat \theta) (\theta - \hat\theta) + \frac{1}{2} l''(\hat \theta) (\theta - \hat\theta)^2 \] Quindi se poniamo \(J(\theta)= -l''(\hat \theta)\) si ha \[ l(\theta) - l(\hat \theta) \simeq - \frac{1}{2} J(\hat \theta) (\theta - \hat\theta)^2. \] La quantità \(J(\hat \theta)\) è l’informazione osservata.

Esempi

Due successi su 10 prove.

s = 2; n = 10; th = .2; j = n/(th*(1-th))
curve(s * log(x) + (n-s) * log(1-x) - s * log(th) - (n-s) * log(1-th) , 0.01, .99, xlab = expression(theta), ylab = "log-lik normalizzata", lwd = 2, col = "blue") 
curve(- j *  (x - th)^2/2, add = TRUE, col = "red", lwd = 2) 
abline(v = s/n, lty = 2)
abline(h = 0, lty = 2)
text(0.6, -25, paste("J = ",j), col = "red")

\(n = 100\), numero di successi \(s = 20\)

s = 20; n = 100; th = .2; j = n/(th*(1-th))
curve(s * log(x) + (n-s) * log(1-x) - s * log(th) - (n-s) * log(1-th) , 0.01, .99, xlab = expression(theta), ylab = "log-lik normalizzata", lwd = 2, col = "blue", ylim = c(-30, 0)) 
curve(- j *  (x - th)^2/2, add = TRUE, col = "red", lwd = 2) 
abline(v = s/n, lty = 2)
abline(h = 0, lty = 2)
text(0.8, -20, paste("J = ",j), col = "red")

\(n = 1000\), numero di successi \(s = 200\)

s = 200; n = 1000; th = .2; j = n/(th*(1-th))
curve(s * log(x) + (n-s) * log(1-x) - s * log(th) - (n-s) * log(1-th) , 0.01, .99, xlab = expression(theta), ylab = "log-lik normalizzata", lwd = 2, col = "blue", ylim = c(-30, 0)) 
curve(- j *  (x - th)^2/2, add = TRUE, col = "red", lwd = 2) 
abline(v = s/n, lty = 2)
abline(h = 0, lty = 2)
text(0.8, -20, paste("J = ",j), col = "red")

Applicazione alla regressione logistica

Per adattare il modello logistico col metodo di massima verosimiglianza si usa la funzione glm. La sigla significa generalized linear model cioè modello lineare generalizzato. Il modello logistico è un modello lineare generalizzato perché estende il modello di regressione lineare in due modi: primo, la distribuzione della risposta può essere non normale; secondo, il modello è lineare, ma per una trasformazione del valor medio, come per esempio la trasformazione logit.

La funzione glm richiede la definizione della formula del modello (che è ancora di tipo lineare) e la specificazione della distribuzione degli errori (Binomiale). Quindi bisogna specificare la funzione che trasforma il predittore lineare nelle probabilità. In modo automatico il programma considera un modello lineare per il logit.

Per adattare il modello ai dati sui difetti del processo di produzione, se Y `e l’indicatore dei difetti e PUR `e l’indice di purezza, si calcola:

purezza = read.table("http://labdisia.disia.unifi.it/gmm/downloads/files/purezza1.txt", header = TRUE)
attach(purezza)
purezza
   PUR Y
1  7.2 0
2  6.3 1
3  8.5 1
4  7.1 0
5  8.2 1
6  4.6 1
7  8.5 0
8  6.9 1
9  8.0 0
10 8.0 1
11 9.1 0
12 6.5 0
13 4.9 1
14 5.3 1
15 7.1 0
16 8.4 1
17 8.5 0
18 6.6 1
19 9.1 0
20 7.1 1
21 7.5 0
22 8.3 0
ml = glm(Y ~ PUR, family = binomial)
summary(ml)

Call:
glm(formula = Y ~ PUR, family = binomial)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.525  -0.944  -0.117   0.972   1.597  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)    6.433      3.510    1.83    0.067 .
PUR           -0.868      0.464   -1.87    0.061 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 30.498  on 21  degrees of freedom
Residual deviance: 25.783  on 20  degrees of freedom
AIC: 29.78

Number of Fisher Scoring iterations: 4

Il programma fornisce le stime di massima verosimiglianza (notare che il metodo numerico per trovare il massimo della verosimiglianza ha richiesto 4 iterazioni), l’errore standard asintotico (Std. Error), la statistica di Wald (z value) e il relativo livello di significatività. Il programma riporta due indici chiamati Null deviance e Residual deviance, che sono le devianze del modello nullo e del modello adattato. I relativi gradi di libertà (degrees of freedom) sono rispettivamente \(n-1\) e \(n-2\). Infine, viene fornito l’indice di Akaike AIC che viene talvolta utilizzato per la scelta del modello.

La devianza con dati binari non si può interpretare come un indice assoluto di bontà di adattamento.

plot(PUR, Y, pch = 20, xlim = c(3,10))
b = coef(ml)
expit = function(x) exp(x)/(1+exp(x))
curve(expit(b[1] + b[2] * x), 3, 10, lwd  =2, col = "RoyalBlue", add = TRUE)

Valori adattati

Probabilità stimate.

fitted(ml)                     
     1      2      3      4      5      6      7      8      9     10     11 
0.5450 0.7236 0.2792 0.5665 0.3345 0.9197 0.2792 0.6085 0.3743 0.3743 0.1871 
    12     13     14     15     16     17     18     19     20     21     22 
0.6875 0.8982 0.8618 0.5665 0.2971 0.2792 0.6686 0.1871 0.5665 0.4801 0.3155

Previsione

new = data.frame(PUR = c(4, 5, 6, 7, 8))
cbind(new, logit = predict(ml, new),  prob = predict(ml, new, type = "response"))
  PUR   logit   prob
1   4  2.9594 0.9507
2   5  2.0911 0.8900
3   6  1.2227 0.7725
4   7  0.3543 0.5877
5   8 -0.5140 0.3743

Intervalli di confidenza

confint(ml)
             2.5 %   97.5 %
(Intercept)  0.495 14.71960
PUR         -1.959 -0.07635

Log-verosimiglianza

liklogi = function(b) {
lin = b[1] + b[2] * PUR
sum(Y * lin - log( 1 + exp(lin)))
}
-2 * liklogi(b)
[1] 25.78

L’ultimo valore fornisce la devianza del modello: \(-2 l(\hat \beta_0, \hat \beta_1)\).

Contour plot

k = 80
B = matrix(0, k,k)
x = seq(0,15, len = k)
y = seq(-2, 0.5, len = k)

for(i in 1:k){
  for(j in 1:k){
    B[i,j] = liklogi(c(x[i], y[j]))
  }
}
contour(x,y,B, nlevels = 40)
points(b[1], b[2], pch = 19, col = "red")

Esercizi 7

  1. Supponiamo che \((12, 3, -3)\) siano tre osservazioni iid da una Cauchy di densità \[ f(y; \theta) = \frac{1}{\pi[1 + (y - \theta)^2]} , \quad y \in \mathbb{R} \] con parametro di posizione \(\theta\). Disegnare la funzione di verosimiglianza \(L(\theta)\).

  2. In un processo industriale è stato predisposto il seguente esperimento per studiare a cosa sono imputabili certi piccoli difetti di produzione. Sono stati scelte balle di materiale grezzo e ciascuna è stata suddivisa in due parti uguali. Quindi, la prima è stata utilizzata in un processo produttivo standard, mentre l’altra è stata lavorata con un metodo alternativo (a temperature più basse). Prima della lavorazione è stato rilevato sull’intera balla un indice di pourezza del materiale. Quindi, a fine produzione, si è registrato se il prodotto derivato da ogni parte lavorata presentava o meno i piccoli difetti studiati. I risultati per 22 balle di materiale sono riportati nella tavola seguente

Purezza    Standard  Modificato Purezza Standard Modificato
 7.2          0          0        6.5       0        1  
 6.3          1          0        4.9       1        1  
 8.5          1          0        5.3       1        0  
 7.1          0          1        7.1       0        1  
 8.2          1          0        8.4       1        0  
 4.6          1          0        8.5       0        1  
 8.5          0          0        6.6       1        0  
 6.9          1          1        9.1       0        0  
 8.0          0          0        7.1       1        0  
 8.0          1          0        7.5       0        1  
 9.1          0          0        8.3       0        0

Si adatti un modello logistico lineare studiando l’effetto della purezza e l’effetto del processo.