MLE / probabilité d'intervalle lognormalement distribué

8

J'ai un ensemble variable de réponses qui sont exprimées sous forme d'intervalle tel que l'exemple ci-dessous.

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

où gauche est la limite inférieure et droite est la limite supérieure de la réponse. Je veux estimer les paramètres en fonction de la distribution log-normale.

Pendant un moment, lorsque j'essayais de calculer directement les probabilités, je me débattais avec le fait que, comme les deux bornes sont réparties le long d'un ensemble de paramètres différent, j'obtenais des valeurs négatives comme ci-dessous:

> 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

Je ne pouvais pas vraiment comprendre comment le résoudre et j'ai décidé d'utiliser le point médian de l'intervalle à la place, ce qui est un bon compromis jusqu'à ce que je trouve la fonction mledist qui extrait la probabilité de log d'une réponse d'intervalle, voici le résumé que j'obtiens:

> 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

Les valeurs des paramètres semblent avoir un sens et la probabilité de log est plus grande que toute autre méthode que j'ai utilisée (distribution à mi-chemin ou distribution de l'une ou l'autre des bornes).

Il y a un message d'avertissement que je ne comprends pas. Quelqu'un pourrait-il me dire si je fais la bonne chose et que signifie ce message?

Appréciez l'aide!

Elio Druml
la source
Votre question revient à "Comment utiliser une fonction R particulière et que signifie ce message d'avertissement?". C'est une question pour StackOverflow plutôt que CrossValidated. De plus, lorsque vous faites référence à une fonction d'un package, vous devez mentionner de quel package elle provient . Dans ce cas, je suppose que vous voulez dire la fonction du package fitdistrplus.
Glen_b -Reinstate Monica
Bienvenue sur le site, @ElioDruml. Je ne peux pas dire si votre question principale est de savoir comment estimer ces paramètres ou quelle est la signification du message d'avertissement. Le premier serait une bonne question pour CV, mais le second est vraiment une question pour Stack Overflow (voir notre FAQ ). Pouvez-vous préciser quelle est votre question principale? Préférez-vous votre séjour Q ici, ou être migré vers SO? (Si ce dernier, marquez votre Q & nous le migrerons pour vous, veuillez ne pas faire de cross-post, cependant .)
gung - Réintégrer Monica

Réponses:

9

Il semble que vous ne calculiez pas correctement la probabilité.

Quand tout ce que vous savez sur une valeur est quex

  1. Il est obtenu indépendamment d'une distribution etFθ

  2. Il se situe entre et inclus (où et sont indépendants de ),ab>abax

alors (par définition) sa probabilité est La probabilité d'un ensemble d'observations indépendantes est donc le produit de telles expressions, une par observation. La log-vraisemblance, comme d'habitude, sera la somme des logarithmes de ces expressions.

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

À titre d'exemple, voici une Rimplémentation où les valeurs de sont dans le vecteur , les valeurs de dans le vecteur et est Lognormal. (Il ne s'agit pas d'une solution polyvalente; en particulier, elle suppose que et pour toutes les données.)aleftbrightFθ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)))
}

Pour trouver la probabilité logarithmique maximale, nous avons besoin d'un ensemble raisonnable de valeurs de départ pour la moyenne logarithmique et l'écart type logarithmique . Cette estimation remplace chaque intervalle par la moyenne géométrique de ses paramètres:μσ

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

Générons des données distribuées de façon aléatoire et lognormale et regroupons-les en intervalles:

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

L'ajustement peut être effectué par un optimiseur multivarié à usage général. (Celui-ci est un minimiseur par défaut, il doit donc être appliqué au négatif de la log-vraisemblance.)

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

6.1188785 0.3957045

L'estimation de est de , non loin de la valeur prévue de , et l'estimation de est de , non loin de la valeur prévue de : pas mal pour seulement valeurs. Pour voir à quel point l'ajustement est bon, traçons la fonction de distribution cumulative empirique et la fonction de distribution ajustée. Pour construire l'ECDF, je viens d'interpoler linéairement à travers chaque intervalle:μ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)

Parcelles

Étant donné que les écarts verticaux sont constamment faibles et varient à la fois de haut en bas, cela ressemble à un bon ajustement.

whuber
la source
Merci beaucoup pour votre contribution @whuber. J'ai recréé votre exemple et tout cela a du sens. Cependant, je n'ai pas pu recréer sur mes propres données de n = 56 dont la tête est gauche <- c (860, 516, 430, 1118, 860, 602) et droite <- c (946, 602, 516 , 1204, 946, 688). J'obtiens ce message d'avertissement: "1: Dans pnorm (log (droite), mu, sigma): NaNs produit 2: Dans pnorm (log (gauche), mu, sigma): NaNs produit" lors de l'ajustement avec l'optimiseur pour extraire le mle des estimations. Cela me ramène à mon problème précédent d'avoir des probabilités négatives lors du calcul. les probabilités étape par étape et soustraction.
Elio Druml
Ce sont les mêmes messages d'avertissement donnés par la fonction mledist du package fitdistrplus. Cependant, comme vous pouvez le voir ci-dessus, cela me donne une sortie pour les estimations mle qui semblent relativement bonnes. Dois-je lui faire confiance et / ou quel est le problème ici? Merci pour les commentaires.
Elio Druml
Pourquoi ne publiez-vous pas vos données, Elio, afin que nous puissions diagnostiquer le problème? Même ainsi, je ne suis pas sûr que ce soient des erreurs critiques. Vous pouvez rencontrer les mêmes problèmes signalés par un autre utilisateur lors de la réduction numérique d'une fonction dans Mathematica ; la même explication peut s'appliquer dans votre cas.
whuber