MLE / Probabilidad de intervalo distribuido lognormalmente

Tengo un conjunto variable de respuestas que se expresan como un intervalo, como la muestra a continuación.

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

donde izquierda es el límite inferior y derecha es el límite superior de la respuesta. Quiero estimar los parámetros de acuerdo con la distribución lognormal.

Durante un tiempo, cuando estaba tratando de calcular las probabilidades directamente, estaba luchando con el hecho de que, dado que los dos límites se distribuyen a lo largo de diferentes conjuntos de parámetros, obtuve algunos valores negativos como a continuación:

> 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

Realmente no pude averiguar cómo resolverlo y decidí usar el punto medio del intervalo, lo cual es un buen compromiso hasta que encontré la función mledist que extrae la probabilidad de una respuesta de intervalo, este es el resumen que obtengo:

> 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

Los valores de los parámetros parecen tener sentido y la probabilidad de logl es mayor que cualquier otro método que haya utilizado (distribución de punto medio o distribución de cualquiera de los límites).

Hay un mensaje de advertencia que no entiendo, ¿alguien podría decirme si estoy haciendo lo correcto y qué significa este mensaje?

Agradezco la ayuda!

Parece que no estás calculando la probabilidad correctamente.

Cuando todos que sabes sobre un valor$x$ es eso

1. Se obtiene independientemente de una distribución. $F_\theta$ y

2. Se encuentra entre $a$ y $b \gt a$ inclusivo (donde $b$ y $a$ son independientes de $x$),

entonces (por definición) su probabilidad es

Como ejemplo, aquí hay un R implementación donde los valores de$a$ están en el vector left , los valores de$b$en el vector righty$F_\theta$es lo normal (Esta no es una solución de propósito general; en particular, supone que$b \gt a$ y $b \ne a$ para todos los datos)

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

Para encontrar la probabilidad de registro máxima, necesitamos un conjunto razonable de valores iniciales para la media del registro $\mu$ y registrar la desviación estándar $\sigma$. Esta estimación reemplaza cada intervalo por la media geométrica de sus puntos finales:

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

Generemos algunos datos aleatorios distribuidos de forma lognormal y los agrupamos en 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

El ajuste puede realizarse mediante un optimizador multivariante de uso general. (Este es un minimizador por defecto, por lo que debe aplicarse al negativo de la probabilidad de registro).

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

6.1188785 0.3957045

La estimación de $\mu$ es $6.12$, no muy lejos del valor previsto de $6$y la estimación de $\sigma$ es $0.40$, no muy lejos del valor previsto de $0.5$: no está mal solo $12$valores. Para ver qué tan bueno es el ajuste, grafiquemos la función empírica de distribución acumulativa y la función de distribución ajustada. Para construir el ECDF, solo interpolo linealmente a través de cada intervalo:

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

Debido a que las desviaciones verticales son consistentemente pequeñas y varían tanto hacia arriba como hacia abajo, parece un buen ajuste.