MLE / Probabilidade de intervalo lognormalmente distribuído

8

Eu tenho um conjunto variável de respostas que são expressas como um intervalo, como a amostra abaixo.

> head(left)
[1]  860  516  430 1118  860  602
> head(right)
[1]  946  602  516 1204  946  688

onde esquerda é o limite inferior e direita é o limite superior da resposta. Quero estimar os parâmetros de acordo com a distribuição lognormal.

Por um tempo, quando eu estava tentando calcular as probabilidades diretamente, estava lutando com o fato de que, como os dois limites são distribuídos por diferentes conjuntos de parâmetros, eu estava obtendo alguns valores negativos, como abaixo:

> Pr_high=plnorm(wta_high,meanlog_high,sdlog_high)
> Pr_low=plnorm(wta_low, meanlog_low,sdlog_low)
> Pr=Pr_high-Pr_low
> 
> head(Pr)
[1] -0.0079951419  0.0001207749  0.0008002343 -0.0009705125 -0.0079951419 -0.0022395514

Eu realmente não conseguia descobrir como resolvê-lo e decidi usar o ponto médio do intervalo, o que é um bom compromisso até encontrar a função mledist que extrai a probabilidade de log de uma resposta de intervalo. Este é o resumo que recebo:

> mledist(int, distr="lnorm")
$estimate
meanlog     sdlog 
6.9092257 0.3120138 

$convergence
[1] 0

$loglik
[1] -152.1236

$hessian
         meanlog       sdlog
meanlog 570.760358    7.183723
sdlog     7.183723 1112.098031

$optim.function
[1] "optim"

$fix.arg
NULL

Warning messages:
1: In plnorm(q = c(946L, 602L, 516L, 1204L, 946L, 688L, 1376L, 1376L,  :
NaNs produced
2: In plnorm(q = c(860L, 516L, 430L, 1118L, 860L, 602L, 1290L, 1290L,  :
NaNs produced

Os valores dos parâmetros parecem fazer sentido e a probabilidade de logaritmo é maior que qualquer outro método que eu usei (distribuição de ponto médio ou distribuição de qualquer um dos limites).

Há uma mensagem de aviso que eu não entendo, então alguém poderia me dizer se estou fazendo a coisa certa e o que essa mensagem significa?

Agradecemos a ajuda!

Elio Druml
fonte
Sua pergunta equivale a "Como uso uma função R específica e o que essa mensagem de aviso significa?". Essa é uma pergunta para StackOverflow, em vez de CrossValidated. Além disso, quando você se refere a uma função de um pacote, deve mencionar de que pacote é . Neste caso, presumo que você queira dizer a função do pacote fitdistrplus.
Glen_b -Reinstala Monica
Bem-vindo ao site, @ElioDruml. Não sei dizer se sua pergunta principal é sobre como estimar esses parâmetros ou qual é o significado da mensagem de aviso. O primeiro seria uma boa pergunta para o CV, mas o último é realmente uma questão para o Stack Overflow (consulte nossas Perguntas frequentes ). Você pode esclarecer qual é a sua principal pergunta? Você prefere que seu Q fique aqui ou seja migrado para o SO? (Se este for o caso, sinalize sua pergunta e nós a migraremos para você , no entanto, não faça postagens cruzadas .)
gung - Reinstate Monica

Respostas:

9

Parece que você pode não estar computando a probabilidade corretamente.

x

  1. Fθ

  2. ab>abax

PrFθ(axb)=Fθ(b)Fθ(a).

RaleftbrightFθb>aba

#
# Lognormal log-likelihood for interval data.
#
lambda <- function(mu, sigma, left, right) {
  sum(log(pnorm(log(right), mu, sigma) - pnorm(log(left), mu, sigma)))
}

μσ

#
# Create an initial estimate of lognormal parameters for interval data.
#
lambda.init <- function(left, right) {
  mid <- log(left * right)/2
  c(mean(mid), sd(mid))
}

Vamos gerar alguns dados aleatórios distribuídos aleatoriamente no log e agrupá-los em intervalos:

set.seed(17)
n <- 12                     # Number of data
z <- exp(rnorm(n, 6, .5))   # Mean = 6, SD = 0.5
left <- 100 * floor(z/100)  # Bin into multiples of 100
right <- left + 100

O ajuste pode ser realizado por um otimizador multivariado de uso geral. (Como este é um minimizador por padrão, deve ser aplicado ao negativo da probabilidade de log).

fit <- optim(lambda.init(left,right), 
             fn=function(theta) -lambda(theta[1], theta[2], left, right))
fit$par

6.1188785 0.3957045

μ6.126σ0.400.512

#
# ECDF of the data.
#
F <- function(x) (1 + mean((abs(x - left) - abs(x - right)) / (right - left)))/2

y <- sapply(x <- seq(min(left) * 0.8, max(right) / 0.8, 1), F)
plot(x, y, type="l", lwd=2, lty=2, ylab="Cumulative probability")
curve(pnorm(log(x), fit$par[1], fit$par[2]), from=min(x), to=max(x), col="Red", lwd=2, 
  add=TRUE)

Parcelas

Como os desvios verticais são consistentemente pequenos e variam para cima e para baixo, parece um bom ajuste.

whuber
fonte
Muito obrigado pela sua contribuição @whuber. Recriei seu exemplo e tudo faz sentido. No entanto, não foi possível recriar com meus próprios dados n = 56 dos quais a cabeça é deixada <- c (860, 516, 430, 1118, 860, 602) e direita <- c (946, 602, 516 1204, 946, 688). Recebo a seguinte mensagem de aviso: "1: Na pnorm (log (direita), mu, sigma): NaNs produzido 2: Na pnorm (log (esquerda), mu, sigma): NaNs produzido" ao ajustar com o otimizador para extrair o estimativas. Isso me traz de volta ao meu problema anterior de ter probabilidades negativas quando calculadas. as probabilidades passo a passo e subtrair.
Elio Druml
Essas são as mesmas mensagens de aviso fornecidas pela função mledist do pacote fitdistrplus. No entanto, como você pode ver acima, ele fornece uma saída para as estimativas de mle que parecem relativamente boas. Devo confiar e / ou qual é o problema aqui? Obrigado pelo feedback.
Elio Druml
Por que você não publica seus dados, Elio, para que possamos diagnosticar o problema? Mesmo assim, não tenho certeza se esses são erros críticos. Você pode estar enfrentando os mesmos problemas relatados por outro usuário ao minimizar numericamente uma função no Mathematica ; a mesma explicação pode ser aplicada no seu caso.
whuber