Frequentist Predictive Distribution für eine Cauchy-Variable

7

Ich konnte dies in der Literatur nicht finden, aber das bedeutet wahrscheinlich, dass ich an der falschen Stelle suche. Ich suche nach der frequentistischen Vorhersageverteilung für eine eindimensionale und eine n-dimensionale Cauchy-Variable, sofern sie existiert.

Das Problem bei der n-dimensionalen Version ist, dass es nichts Vergleichbares wie eine Kovariatenmatrix gibt, sondern nur einen Skalenparameter, der die Fehler hyperzirkular macht. Ich konnte sehen, dass dies die Existenz eines zentralen Wertes störte.

BEARBEITEN

Ich möchte entweder aus einer Reihe von Beobachtungen vorhersagen, die aus einer Cauchy-Verteilung mit Zentrum und Skala oder aus einer Gleichung vorhersagen wobei wie oben aus einer Cauchy-Verteilung gezogen wird. Es könnte ein Vektor oder mehrdimensional sein, aber ich versuche, die relativen Eigenschaften der Bayes'schen gegenüber der häufig auftretenden Vorhersage zu bestimmen. Meine Daten stammen entweder aus einem abgeschnittenen Cauchy oder einem Cauchy, je nachdem welcher Satz.xi+1x1xiμσ,yi+1y=mx+b,x

Ein Vorhersageintervall funktioniert, da ich das Intervall nur auf 100% setze.

Dave Harris
quelle
1
Wenn Sie Cauchy variate sagen , meinen Sie einen Regressor in der Regressionsgleichung , und suchen Sie das Vertrauen Intervall der OLS-Schätzung für ? Ich bin mir nicht sicher, ob ich die Frage verstehe. Xj,iYi=j=1JXj,iβj+εiβj
Jeremias K
1
Es gibt Arbeiten wie diese aus dem Jahr 2008 zur Schätzung mit multivariaten t-Verteilungen. Das multivariate Cauchy ist ein Sonderfall des multivariaten t. Abgesehen davon ermöglicht dies eine vollständig flexible Korrelationsstruktur. Hilft Ihnen das oder passt es nicht zu Ihrer Frage?
eric_kernfeld
1
Versuchen Sie, dies unter stats.stackexchange.com/questions/16349 für eine multivariate Verteilung mit Nullmitteln durchzuführen ?
Sextus Empiricus
1
@eric_kernfeld Ich muss es sorgfältig lesen, aber ja, das ist eine Art davon, außer dass ich wissen möchte, wie man seine prädiktive Dichte mit Frequentist-Methoden findet.
Dave Harris
1
Es scheint, dass Sie versuchen, Cauchy-Verteilungsparameter aus zu schätzen . Ist das richtig? xi
Aksakal

Antworten:

2

Die allgemeine Lösung für Ihr Problem ist die Maximum Likelihood Estimation (MLE) Ihrer Parameter . Sobald sie als , ersetzen Sie die unbekannten Parameter durch Ihr PDF, dh Sie schätzen das PDF Ihrer Zufallsvariablen als . Auf diese Weise können Sie die prädiktive Verteilung Ihrer Cauchy-Zufallsvariablen erstellen. θθ^f^(xi)=f(xi|θ^)

Für den univariaten Fall ist dieses Papier eine ausgezeichnete Ressource . Für das univariate Cauchy mit Zentrum und Skala hat man eine geschlossene Form, wenn Sie Beobachtungen haben. Wenn Sie Beobachtungen haben, existiert die MLE . Wenn Sie Beobachtungen, haben Sie zwei Gleichungen zu lösen , die durch Einstellen der erste Ableitung des Log-Likelihood auf Null leicht abgeleitet werden, finden hier für ihre genaue Form. (In ihrer Notation ist und .) Die numerische Lösung dieses Problems hat eine Implementierung in der R-Sprache, siehe hier .μσ34n>4nx0=μσ=γ

Für den multivariaten Fall müssen Sie lediglich beachten, dass die multivariate Cauchy-Verteilung einfach eine multivariate Verteilung ist, bei der der Freiheitsgradparameter auf , wie bereits in den Kommentaren ausgeführt wurde. Für die multivarate- , können Sie MLE Inferenz tun , wie hervorragend in erklärt diese Antwort , die auf dem Papier basiert, dass eric_kernfeld hingewiesen hat. Ich habe keine einsatzbereite Implementierung für diesen Algorithmus gefunden, aber wie Sie sehen werden, wenn Sie sich die bereitgestellte Antwort im Beitrag ansehen, sollte es wirklich einfach sein, sie selbst zu implementieren.t1t

Unterschied zur Bayes'schen Vorhersage : In der Bayes'schen Einstellung würden Sie den Parametern und einen Vorrang und Ihre Unsicherheit darüber als Zufallsvariable modellieren. Auf diese Weise erhalten Sie für beide Parameter hintere Verteilungen, die die relative Sicherheit angeben, die Sie angesichts Ihrer Daten über sie haben. Wenn Sie das hintere , erhalten Sie Ihre Vorhersageverteilung als , integriert Ihre Unsicherheit. Im Gegensatz dazu erhalten Sie mit der MLE-Einstellung Punktschätzungen für undμσq(μ,σ|x1,,xn)f(x|μ,σ)q(μ,σ|x1,,xn)dμdσμσdass Sie sich in die Funktionsform Ihres PDFs einfügen. Entsprechend könnte man sagen, dass MLE zu einem Posterior mit der Punktmasse am Tupel und einer Wahrscheinlichkeit von bei jedem anderen Wert führt. Daher ignorieren Sie in diesem Fall alle Parameterunsicherheiten und verlassen sich auf die Tatsache, dass asymptotisch , was bedeutet, dass (gleichmäßig über) ).1(μ^,σ^)0θ^θf^(x)f(x)x

Nun, es sei denn, für den exotischen Fall, in dem ist und Ihrer Beobachtungen den Wert annehmen, nimmt die andere Hälfte den Wert , was mit der Wahrscheinlichkeit Null geschieht, weil die Cauchy-Verteilung kontinuierlich ist.nn/2x1x2

Jeremias K.
quelle
Jeremias. Glauben Sie, dass es Möglichkeiten gibt, die Unsicherheit über in das Vorhersageintervall einzubeziehen? Und wie konstruieren wir eine prädiktive Verteilung aus dem PDF einer mehrdimensionalen Cauchy-Verteilung? θ^
Sextus Empiricus
Wenn Sie davon ausgehen, dass es sich bei um zufällige Ziehungen aus einer Cauchy-Zufallsvariablen mit unbekannten Parametern handelt, erhalten Sie durch direktes Einfügen der geschätzten Parameter in die Funktionsform die prädiktive Verteilung der nächsten Ziehungen von . xixi
Jeremias K
Wenn Sie die Parameterunsicherheit berücksichtigen möchten, müssen Sie den Bayes'schen Weg gehen. Beachten Sie, dass sich der Parameter posterior der Bayes'schen Inferenz als Nebenprodukt asymptotisch über den Satz von Bernstein Mises auf die MLE konzentriert.
Jeremias K
1

Man könnte eine Monte-Carlo-Methode verwenden, um empirische Schätzungen für Beziehungen zwischen und dem Vorhersageintervall für .x1....xixi+n

Motivation: Wenn wir das Vorhersageintervall basierend auf den Quartilen / CDF einer Verteilung schätzen, die sich aus Schätzungen der maximalen Wahrscheinlichkeit (oder anderen Arten von Parameterschätzungen) ergibt, unterschätzen wir die Größe des Intervalls. In der Praxis fällt der Punkt tatsächlich häufiger als vorhergesagt aus dem Bereich heraus.xi+n

Die folgende Abbildung zeigt, um wie viel wir die Größe des Intervalls unterschätzen, indem wir ausdrücken, wie oft eine neue Messung außerhalb des Vorhersagebereichs liegt, basierend auf Parameterschätzungen. (basierend auf Berechnungen mit 2000 Wiederholungen für die Vorhersage)xi

Wenn wir beispielsweise ein Vorhersageintervall von 99% verwenden (wodurch 1% Fehler erwartet werden), erhalten wir fünfmal mehr Fehler, wenn die Stichprobengröße 3 betrug.

Diese Art von Berechnungen kann verwendet werden, um empirische Beziehungen herzustellen, wie wir den Bereich korrigieren können, und die Berechnungen zeigen, dass für große die Differenz kleiner wird (und irgendwann kann man sie für irrelevant halten).n

Differenz zwischen MLE-Schätzung und effektivem Konfidenzintervall

set.seed(1)

# likelihood calculation
like<-function(par, x){
  scale = abs(par[2])
  pos   = par[1]
  n <- length(x)
  like <- -n*log(scale*pi) - sum(log(1+((x-pos)/scale)^2))
  -like
}

# obtain effective predictive failure rate rate
tryf <- function(pos, scale, perc, n) {

  # random distribution
  draw <- rcauchy(n, pos, scale)

  # estimating distribution parameters based on median and interquartile range
  first_est <- c(median(draw), 0.5*IQR(draw))

  # estimating distribution parameters based on likelihood
  out <- optim(par=first_est, like, method='CG', x=draw)
  # making scale parameter positive (we used an absolute valuer in the optim function)
  out$par[2] <- abs(out$par[2])

  # calculate predictive interval
  ql <- qcauchy(perc/2, out$par[1], out$par[2])
  qh <- qcauchy(1-perc/2, out$par[1], out$par[2])

  # calculate effective percentage outside predicted predictive interval
  pl <- pcauchy(ql, pos, scale)
  ph <- pcauchy(qh, pos, scale)
  error <- pl+1-ph
  error
}

# obtain mean of predictive interval in 2000 runs
meanf <- function(pos,scale,perc,n) {
  trueval <- sapply(1:2000,FUN <- function(x) tryf(pos,scale,perc,n))
  mean(trueval)
}


#################### generate image

# x-axis chosen desired interval percentage
percentages <- 0.2/1.2^c(0:30)

# desired sample sizes n
ns <- c(3,4,5,6,7,8,9,10,20,30)

# computations
y <- matrix(rep(percentages, length(ns)), length(percentages))
for (i in which(ns>0)) {
  y[,i] <- sapply(percentages, FUN <- function(x) meanf(0,1,x,ns[i]))
}

# plotting
plot(NULL,
     xlim=c(0.0008,1), ylim=c(0,10),
     log="x",
     xlab="aimed error rate",
     ylab="effective error rate / aimed error rate",
     yaxt="n",xaxt="n",axes=FALSE)
axis(1,las=2,tck=-0.0,cex.axis=1,labels=rep("",2),at=c(0.0008,1),pos=0.0008)
axis(1,las=2,tck=-0.005,cex.axis=1,at=c(0.001*c(1:9),0.01*c(1:9),0.1*c(1:9)),labels=rep("",27),mgp=c(1.5,1,0),pos=0.0008)
axis(1,las=2,tck=-0.01,cex.axis=1,labels=c(0.001,0.01,0.1,1), at=c(0.001,0.01,0.1,1),mgp=c(1.5,1,0),pos=0.000)
#axis(2,las=1,tck=-0.0,cex.axis=1,labels=rep("",2),at=c(0.0008,1),pos=0.0008)
#axis(2,las=1,tck=-0.005,cex.axis=1,at=c(0.001*c(1:9),0.01*c(1:9),0.1*c(1:9)),labels=rep("",27),mgp=c(1.5,1,0),pos=0.0008)
#axis(2,las=1,tck=-0.01,cex.axis=1,labels=c(0.001,0.01,0.1,1), at=c(0.001,0.01,0.1,1),mgp=c(1.5,1,0),pos=0.0008)
axis(2,las=2,tck=-0.01,cex.axis=1,labels=0:15, at=0:15,mgp=c(1.5,1,0),pos=0.0008)


colours <- hsv(c(1:10)/20,1,1-c(1:10)/15)
for (i in which(ns>0)) {
  points(percentages,y[,i]/percentages,pch=21,cex=0.5,col=colours[i],bg=colours[i])
}

legend(x=0.4,y=4.5,pch=21,legend=ns,col=colours,pt.bg=colours,title="sample size")

title("difference between confidence interval and effective confidence interval")


plot(ns,y[31,]/percentages[31],log="")
Sextus Empiricus
quelle
Was sagt uns das Diagramm außer der Verwendung einer kleinen Stichprobengröße, was zu einer schlechten Schätzung Ihrer Parameter bei Verwendung von mle führt ? Ich kann nicht erkennen, wie es mit mle ungültig wird, da die Fehlerraten selbst bei einer sehr kleinen Stichprobengröße von 30 hervorragend aussehen. Ich bin mir auch nicht sicher, ob ich die von Ihnen vorgeschlagene Alternative verstehe. Würde es Ihnen etwas ausmachen, die Berechnungsmethoden zu erweitern? Sie erwähnen am Anfang Ihrer Antwort?
Jeremias K
1
@JeremiasK In praktischen Anwendungen mit kleinen Stichprobengrößen könnte man diese Berechnungen als empirisch ermittelte Korrekturfaktoren verwenden.
Sextus Empiricus
Das macht Sinn! Ich glaube nicht, dass Sie es in der Post erwähnen, vielleicht sollten Sie es bearbeiten, damit die Leute die Kommentare nicht durchlesen müssen
Jeremias K
@MartijnWeterings bisher machen Sie am sinnvollsten. Der Drehpunktn(μ^μ)σ^$ folgt dem Standardnormal, sobald die Stichprobengröße ungefähr 100 erreicht hat, aber ich habe festgestellt, dass ich nicht in der Lage bin, dies abzuwickeln, da ich anstelle einer Variablen eine Funktion für die Minimierung auswähle und dies zuvor noch nicht getan habe.
Dave Harris
@ DaveHarris Ich glaube, meine Methode unterscheidet sich nicht so sehr von der von Jeremia, außer dass ich einen Ausdruck (und nur durch einen experimentellen mathematischen Ansatz) für den unterschätzten Bereich mache, der aufgrund der Verteilung auftritt f(x,x^0,γ^) ist eine überstreute Version von f(x,x0,γ).
Sextus Empiricus
0

Es scheint, dass Sie lediglich die Parameter der Cauchy-Verteilung aus dem Datensatz abschätzen müssen xi. Hier ist, was Stephens vorschlägt, es ist nicht MLE, und der Autor behauptet, diese Methode sei konsistent und stabiler als MLE, obwohl Sie berücksichtigen müssen, dass dies im letzten Jahrhundert geschrieben wurde.

Geben Sie hier die Bildbeschreibung ein

wobei Cauchy wie folgt parametrisiert wird: Geben Sie hier die Bildbeschreibung ein

Sobald Sie die Verteilung haben, wird Ihre Punktprognose sein α^. Beachten Sie, dass Sie, da es keine Momente gibt, nicht zeigen können, dass Ihre Prognose im üblichen Sinne optimal ist, z. B. um die erwarteten Quadratkosten zu minimieren.

Aksakal
quelle