MLE / Wahrscheinlichkeit eines logarithmisch normal verteilten Intervalls

8

Ich habe einen variablen Satz von Antworten, die als Intervall ausgedrückt werden, wie im folgenden Beispiel.

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

Dabei ist links die Untergrenze und rechts die Obergrenze der Antwort. Ich möchte die Parameter anhand der logarithmischen Normalverteilung schätzen.

Als ich eine Weile versuchte, die Wahrscheinlichkeiten direkt zu berechnen, hatte ich Probleme damit, dass ich einige negative Werte wie die folgenden erhielt, da die beiden Grenzen auf verschiedene Parametersätze verteilt sind:

> 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

Ich konnte nicht wirklich herausfinden, wie ich es lösen sollte, und entschied mich stattdessen für die Verwendung des Mittelpunkts des Intervalls, was ein guter Kompromiss ist, bis ich eine mledist-Funktion fand, die die Loglikelihood einer Intervallantwort extrahiert. Dies ist die Zusammenfassung, die ich erhalte:

> 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

Die Parameterwerte scheinen sinnvoll zu sein und die Loglikelihood ist größer als bei jeder anderen Methode, die ich verwendet habe (Mittelpunktverteilung oder Verteilung einer der Grenzen).

Es gibt eine Warnmeldung, die ich nicht verstehe. Kann mir jemand sagen, ob ich das Richtige tue und was diese Meldung bedeutet?

Schätzen Sie die Hilfe!

Elio Druml
quelle
Ihre Frage lautet "Wie verwende ich eine bestimmte R-Funktion und was bedeutet diese Warnmeldung?". Das ist eher eine Frage für StackOverflow als für CrossValidated. Wenn Sie auf eine Funktion aus einem Paket verweisen, sollten Sie außerdem angeben, aus welchem ​​Paket sie stammt . In diesem Fall meine ich vermutlich die Funktion aus dem Paket fitdistrplus.
Glen_b -Rate State Monica
Willkommen auf der Website @ElioDruml. Ich kann nicht sagen, ob Ihre Hauptfrage darin besteht, wie diese Parameter geschätzt werden oder welche Bedeutung die Warnmeldung hat. Ersteres wäre eine gute Frage für den Lebenslauf, letzteres ist wirklich eine Frage für den Stapelüberlauf (siehe unsere FAQ ). Können Sie klarstellen, was Ihre Hauptfrage ist? Möchten Sie Ihren Q-Aufenthalt hier bevorzugen oder nach SO migriert werden? (Wenn letzteres der Fall ist, kennzeichnen Sie Ihr Q & wir werden es für Sie migrieren, aber bitte nicht überkreuzen .)
gung - Reinstate Monica

Antworten:

9

Es hört sich so an, als würden Sie die Wahrscheinlichkeit möglicherweise nicht richtig berechnen.

x

  1. Fθ

  2. ab>abax

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

Als Beispiel ist hier eine RImplementierung, bei der die Werte von im Vektor , die Werte von im Vektor und Lognormal sind. (Dies ist keine Allzwecklösung. Insbesondere wird davon ausgegangen, dass und für alle Daten gelten.)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)))
}

Um die maximale Log-Wahrscheinlichkeit zu ermitteln, benötigen wir einen angemessenen Satz von Startwerten für den Log-Mittelwert und die Log-Standardabweichung . Diese Schätzung ersetzt jedes Intervall durch das geometrische Mittel seiner Endpunkte:μσ

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

Lassen Sie uns einige zufällig logarithmisch verteilte Daten generieren und diese in Intervalle unterteilen:

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

Die Anpassung kann durch einen universellen multivariaten Optimierer durchgeführt werden. (Dieser ist standardmäßig ein Minimierer , daher muss er auf das Negativ der Log-Wahrscheinlichkeit angewendet werden.)

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

6,1188785 0,3957045

Die Schätzung von ist , nicht weit vom beabsichtigten Wert von , und die Schätzung von ist , nicht weit vom beabsichtigten Wert von : nicht schlecht für nur Werte. Um zu sehen, wie gut die Anpassung ist, zeichnen wir die empirische kumulative Verteilungsfunktion und die angepasste Verteilungsfunktion auf. Um das ECDF zu konstruieren, interpoliere ich einfach linear durch jedes Intervall:μ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)

Grundstücke

Da die vertikalen Abweichungen konstant klein sind und sowohl nach oben als auch nach unten variieren, scheint dies eine gute Anpassung zu sein.

whuber
quelle
Vielen Dank für Ihre Eingabe @whuber. Ich habe Ihr Beispiel neu erstellt und alles macht Sinn. Ich konnte jedoch meine eigenen Daten von n = 56 nicht neu erstellen, von denen der Kopf links <- c (860, 516, 430, 1118, 860, 602) und rechts <- c (946, 602, 516) ist 1204, 946, 688). Ich erhalte die folgende Warnmeldung: "1: In pnorm (log (rechts), mu, sigma): NaNs erzeugt 2: In pnorm (log (links), mu, sigma): NaNs erzeugt", wenn mit dem Optimierer zum Extrahieren der mle Schätzungen. Das bringt mich zurück zu meinem früheren Problem, negative Wahrscheinlichkeiten bei der Berechnung zu haben. die Wahrscheinlichkeiten Schritt für Schritt und subtrahieren.
Elio Druml
Dies sind die gleichen Warnmeldungen, die von der Funktion mledist aus dem Paket fitdistrplus ausgegeben werden. Wie Sie oben sehen können, gibt es mir jedoch eine Ausgabe für die mle-Schätzungen, die relativ gut aussehen. Soll ich ihm vertrauen und / oder worum geht es hier? Danke für die Rückmeldung.
Elio Druml
Warum postest du deine Daten nicht, Elio, damit wir das Problem diagnostizieren können? Trotzdem bin ich mir nicht sicher, ob dies kritische Fehler sind. Möglicherweise treten dieselben Probleme auf , die von einem anderen Benutzer gemeldet wurden, wenn Sie eine Funktion in Mathematica numerisch minimieren . Die gleiche Erklärung könnte in Ihrem Fall gelten.
whuber