R Estimating Parameter der Binomialverteilung

Question

I Schätzung ist versucht, Parameter n und p von Binomialverteilung von Maximum Likelihood in R.

ich die Funktion optim von Statistik-Paket bin mit, aber es gibt eine Error.

Das ist mein Code:

xi = rbinom(100, 20, 0.5) # Sample 
n = length(xi) # Sample size 

# Log-Likelihood 
lnlike <- function(theta){ 
log(prod(choose(theta[1],xi))) + sum(xi*log(theta[2])) + 
(n*theta[1] - sum(xi))*log(1-theta[2]) 
} 

# Optimizing 
optim(theta <- c(10,.3), lnlike, hessian=TRUE) 

Fehler bei Optim (Theta < - c (10, 0,3), lnlike, Hessische = TRUE): Funktion bei Anfangsparameter werden nicht

ausgewertet

Wer hat das getan? Welche Funktion benutzt?

Answer 1

tl; dr Sie gehen eine Wahrscheinlichkeit von Null (und damit ein negativen unendlichen Log-Likelihood) zu erhalten, wenn die Reaktionsvariable größer ist als die binomische N (was der theoretische Maximum Wert von die Antwort). Bei den meisten praktischen Problemen wird N als bekannt vorausgesetzt und nur die Wahrscheinlichkeit geschätzt. Wenn Sie N schätzen wollen, müssen Sie (1) festlegen, dass es> = der größte Wert in der Stichprobe ist; (2) Etwas Spezielles zu tun, um über einen Parameter zu optimieren, der diskret sein muss (dies ist ein fortgeschrittenes/kniffliges Problem).

Der erste Teil dieser Antwort zeigt Debugging-Strategien zur Identifizierung des Problems, der zweite zeigt eine Strategie zur simultanen Optimierung von N und p (durch Brute Force über einen vernünftigen Bereich von N).

Setup:

set.seed(101) 
n <- 100 
xi <- rbinom(n, size=20, prob=0.5) # Sample 

Log-Likelihood-Funktion:

lnlike <- function(theta){ 
    log(prod(choose(theta[1],xi))) + sum(xi*log(theta[2])) + 
     (n*theta[1] - sum(xi))*log(1-theta[2]) 
} 

Lassen Sie uns diese brechen.

theta <- c(10,0.3) ## starting values 
lnlike(c(10,0.3)) ## -Inf 

OK, die Log-Likelihood ist -Inf an dem Startwert. Kein Wunder, dass optim() damit nicht arbeiten kann.

Lassen Sie uns die Bedingungen durcharbeiten.

OK, wir sind bereits in der ersten Amtszeit in Schwierigkeiten.

prod(choose(theta[1],xi)) ## 0 

Das Produkt ist Null ... warum?

choose(theta[1],xi) 
## [1] 120 210 10 0 0 10 120 210 0 0 45 210 1 0 

Viele Nullen. Warum? Was sind die Werte von xi, die problematisch sind?

## [1] 7 6 9 12 11 9 7 6 

Aha! Wir sind OK für 7, 6, 9 ... aber in Schwierigkeiten mit 12.

badvals <- (choose(theta[1],xi)==0) 
all(badvals==(xi>10)) ## TRUE 

Wenn Sie das wirklich tun möchten, können Sie es durch Brute-Force-Enumeration über vernünftige Werte von n tun können .. .

## likelihood function 
llik2 <- function(p,n) { 
    -sum(dbinom(xi,prob=p,size=n,log=TRUE)) 
} 
## possible N values (15 to 50) 
nvec <- max(xi):50 
Lvec <- numeric(length(nvec)) 
for (i in 1:length(nvec)) { 
    ## optim() wants method="Brent"/lower/upper for 1-D optimization 
    Lvec[i] <- optim(par=0.5,fn=llik2,n=nvec[i],method="Brent", 
        lower=0.001,upper=0.999)$val 
} 
nvec[which.min(Lvec)] ## 20 
par(las=1,bty="l") 
plot(nvec,Lvec,type="b") 

Answer 2

Warum Sie in Schwierigkeiten geraten?

Wenn Sie lnlike(c(10, 0.3)) tun, erhalten Sie -Inf. Deshalb lautet die Fehlermeldung lnlike statt optim.

Oft ist n bekannt, und nur muss geschätzt werden. In dieser Situation ist entweder der Momentschätzer oder der Maximum-Likelihood-Schätzer in geschlossener Form, und es ist keine numerische Optimierung erforderlich. Es ist wirklich seltsam, n zu schätzen.

Wenn Sie schätzen möchten, müssen Sie sich bewusst sein, dass es eingeschränkt ist. Prüfen

range(xi) ## 5 15 

Sie Beobachtungen haben Bereich [5, 15] daher kommt es, dass n >= 15 erforderlich. Wie können Sie einen Anfangswert von 10 übergeben? Die Suchrichtung für n, sollte von einem großen Anfangswert sein, und dann allmählich nach unten suchen, bis es max(xi) erreicht. Sie könnten also 30 als Anfangswert für n versuchen.

Zusätzlich müssen Sie lnlike nicht in der aktuellen Weise definieren. Tun Sie dies:

lnlike <- function(theta, x) -sum(dbinom(x, size = theta[1], prob = theta[2], log = TRUE)) 

• optim wird häufig zur Minimierung verwendet (obwohl es Maximierung tun können). Ich habe ein Minuszeichen in die Funktion gesetzt, um eine negative Log-Wahrscheinlichkeit zu erhalten. Auf diese Weise minimieren Sie lnlike w.r.t. theta.
• Sie sollten auch Ihre Beobachtungen xi als zusätzliches Argument an lnlike weiterleiten, anstatt es aus der globalen Umgebung zu nehmen.

Naive try mit optim:

In meinem Kommentar, ich schon gesagt, dass ich nicht glaube, optim mit n arbeiten zu schätzen, weil n ganze Zahlen sein müssen, während optim für kontinuierliche Variablen verwendet wird, . Diese Fehler und Warnungen sollen Sie überzeugen.

optim(c(30,.3), fn = lnlike, x = xi, hessian = TRUE) 

Error in optim(c(30, 0.3), fn = lnlike, x = xi, hessian = TRUE) : 
    non-finite finite-difference value [1] 
In addition: There were 15 or more warnings (use warnings() to see the 
first 15 

> warnings() 
Warning messages: 
1: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 
2: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 
3: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 
4: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 
5: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 

Lösung?

Ben hat Ihnen einen Weg gegeben. Anstatt optimn zu schätzen, führen wir manuell eine Rastersuche nach n durch. Für jeden Kandidaten n führen wir eine univariate Optimierung w.r.t. p. (Ups, in der Tat, es gibt keine Notwendigkeit, numerische Optimierung hier zu tun.) Auf diese Weise erhalten Sie eine Profilwahrscheinlichkeit von n. Dann finden wir n auf dem Raster, um diese Profilwahrscheinlichkeit zu minimieren.

Ben hat Ihnen alle Details zur Verfügung gestellt, und ich werde das nicht wiederholen. Schöne (und schnelle) Arbeit, Ben!