R Оценка параметров биномиального распределения

Question

Я пытаюсь оценить параметры п и р из биномиального распределения методом максимального правдоподобия в R.

Я использую функцию optim из пакета статистики, но есть ошибка.

То есть мой код:

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) 

Ошибка в Optim (тета < - с (10, 0.3), lnlike, мешковина = TRUE): функция не может быть оценена при начальных параметрах

Кто-нибудь это сделал? Какая функция используется?

Answer 1

Т.Л., д-р вы собираетесь получить вероятность нуля (и, следовательно, отрицательно-бесконечное журнала правдоподобия), если переменная отклики больше, чем биномиальный N (который является теоретическим максимального значения ответ). В большинстве практических задач N принимается как известно, и оценивается вероятность. Если вы хотите оценить N, вам нужно (1) ограничить это> = наибольшее значение в выборке; (2) сделать что-то особенное для оптимизации по параметру, который должен быть дискретным (это сложная/сложная проблема).

В первой части этого ответа показаны стратегии отладки для идентификации проблемы, вторая иллюстрирует стратегию оптимизации N и p одновременно (путем грубой силы в разумном диапазоне N).

Установка:

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

Log-функция правдоподобия:

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

Давайте разберем это вниз.

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

OK, логарифмическая функция правдоподобия является -Inf в стартовое значение. Не удивительно, что optim() не может с этим работать.

Давайте рассмотрим условия.

log(prod(choose(theta[1],xi))) ## -Inf 

ОК, у нас уже есть проблемы на первом сроке.

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

Продукт равен нулю ... почему?

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

Много нулей. Зачем? Каковы значения xi, которые являются проблематичными?

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

Aha! Мы хорошо для 7, 6, 9 ... но в ладах с 12.

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

Если вы действительно хотите сделать это, вы можете сделать это с помощью грубой силы перечисления более разумных значений n .. ,

## 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

Почему вы попадаете в неприятности?

Если вы делаете lnlike(c(10, 0.3)), вы получаете -Inf. Вот почему ваше сообщение об ошибке жалуется lnlike, а не optim.

Часто известен n, и только p нуждается в оценке. В этой ситуации оценка момента времени или оценка максимального правдоподобия находятся в закрытом виде, и никакая численная оптимизация не требуется. Итак, это действительно странно оценивать n.

Если вы хотите оценить, вы должны знать, что он ограничен. Проверить

range(xi) ## 5 15 

Вы наблюдения имеют диапазон [5, 15], следовательно, требуется, чтобы n >= 15. Как вы можете передать начальное значение 10? Поисковое направление для n, должно быть от большого стартового значения, а затем постепенно искать вниз, пока оно не достигнет max(xi). Итак, вы можете попробовать 30 в качестве начального значения для n.

Кроме того, текущим способом не нужно определять lnlike. Сделайте это:

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

• optim часто используется для минимизации (хотя это может сделать максимизацию). Я положил знак минус в функцию, чтобы получить отрицательную вероятность регистрации. Таким образом, вы минимизируете lnlike w.r.t. theta.
• Вы также должны передать свои замечания xi в качестве дополнительного аргумента lnlike, а не принимать его из глобальной среды.

Наивная попытка с optim:

В моем комментарии, я уже сказал, что я не верю, что с помощью optim оценить n будет работать, потому что n должны быть целыми числа, а optim используется для непрерывных переменных , Эти ошибки и предупреждения должны вас убедить.

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 

Решение?

Бен предоставил вам способ. Вместо того, чтобы дать optim для оценки n, мы вручную проведем поиск сетки для n. Для каждого кандидата n мы выполняем одномерную оптимизацию w.r.t. p. (Ой, на самом деле, здесь нет необходимости делать численную оптимизацию.) Таким образом, вы получаете вероятность профиля n. Затем мы найдем n на сетке, чтобы минимизировать вероятность этого профиля.

Бен предоставил вам полную информацию, и я не буду повторять этого. Хорошая (и быстрая) работа, Бен!