Einbindung einer empirischen CDF

13

Ich habe eine empirische Verteilung G(x) . Ich berechne es wie folgt

    x <- seq(0, 1000, 0.1)
    g <- ecdf(var1)
    G <- g(x)

Ich bezeichne h(x)=dG/dx , dh h ist das pdf, während G das cdf ist.

Ich möchte nun eine Gleichung für die obere Integrationsgrenze (sagen wir a ) lösen , so dass der erwartete Wert von x etwas k .

Das heißt, wenn ich von nach b integriere , hätte ich x h ( x ) d x = k . Ich möchte nach b auflösen .0bxh(x)dx=kb

Nach Teilen integrierend, kann ich die Gleichung als umschreiben

, wobei das Integral von 0 bis b ist ------- (1)bG(b)0bG(x)dx=k0b

Ich denke, ich kann das Integral wie folgt berechnen

    intgrl <- function(b) {
        z <- seq(0, b, 0.01)
        G <- g(z)
        return(mean(G))
     }

Aber wenn ich versuche, diese Funktion mit zu verwenden

    library(rootSolve)
    root <- uniroot.All(fun, c(0, 1000))

Wo Spaß gleich (1) ist, erhalte ich die folgende Fehlermeldung

    Error in seq.default(0, b, by = 0.01) : 'to' must be of length 1  

Ich denke, das Problem ist, dass meine Funktion intgrlmit einem numerischen Wert ausgewertet wird, während uniroot.Alldas Intervall überschritten wirdc(0,1000)

Wie soll ich für in dieser Situation in R lösen ?b

user46768
quelle

Antworten:

13

Die sortierten Daten seien . Betrachten Sie zum Verständnis der empirischen CDF G einen der Werte des x i - -Letters als γ - und nehmen Sie an, dass eine Zahl k des x vorliegtx1x2xnGxiγk kleiner als γ sind und t 1 von x i gleich γ sind . Wählen Sie ein Intervall [ α , β ], in dem von allen möglichen Datenwerten nur γ giltxiγt1xiγ[α,β]γerscheint. Innerhalb dieses Intervalls hat dann definitionsgemäß den konstanten Wert k / n für Zahlen kleiner als γ und springt auf den konstanten Wert ( k + t ) / n für Zahlen größer als γ .Gk/nγ(k+t)/nγ

ECDF

Betrachten Sie den Beitrag zu aus dem Intervall [ α ,0bxh(x)dx . Obwohl h nicht eine Funktion ist- es ist ein Punkt Maß der Größe t / n an γ --die Integral wirddefiniertdurch Integration von Teilen, umUmrechnung in eine ehrlichen-to-Güte integral. Machen wir das über das Intervall [ α , β ] :[α,β]ht/nγ[α,β]

αβxh(x)dx=(xG(x))|αβαβG(x)dx=(βG(β)αG(α))αβG(x)dx.

Der neue Integrand ist integrierbar , obwohl er bei diskontinuierlich ist. Sein Wert kann leicht gefunden werden, indem die Integrationsdomäne in die Teile vor und nach dem Sprung in G aufgeteilt wird :γG

αβG(x)dx=αγG(α)dx+γβG(β)dx=(γα)G(α)+(βγ)G(β).

Ersetzt man dies in das Vorstehende und erinnert man sich an ergibt sich G ( β ) = ( k + t ) / nG(α)=k/n,G(β)=(k+t)/n

αβxh(x)dx=(βG(β)αG(α))((γα)G(α)+(βγ)G(β))=γtn.

Mit anderen Worten multipliziert dieses Integral die Position (entlang der Achse) jedes Sprungs mit der Größe dieses Sprungs. Die Größe des Sprungs istX

tn=1n++1n

mit einem Term für jeden der Datenwerte, der entspricht . Das Addieren der Beiträge von allen solchen Sprüngen von G zeigt dasγG

0bxh(x)dx=i:0xib(xi1n)=1nxibxi.

1/n[0,b]1/n1/mm[0,b]

kb1nxibxi=k.kj

1ni=1j1xik<1ni=1jxi,

you will have narrowed b to the interval [xj1,xj). You can do no better than that using the ECDF. (By fitting some continuous distribution to the ECDF you can interpolate to find an exact value of b, but its accuracy will depend on the accuracy of the fit.)


R performs the partial sum calculation with cumsum and finds where it crosses any specified value using the which family of searches, as in:

set.seed(17)
k <- 0.1
var1 <- round(rgamma(10, 1), 2)
x <- sort(var1)
x.partial <- cumsum(x) / length(x)
i <- which.max(x.partial > k)
cat("Upper limit lies between", x[i-1], "and", x[i])

The output in this example of data drawn iid from an Exponential distribution is

Upper limit lies between 0.39 and 0.57

The true value, solving 0.1=0bxexp(x)dx, is 0.531812. Its closeness to the reported results suggests this code is accurate and correct. (Simulations with much larger datasets continue to support this conclusion).

Here is a plot of the empirical CDF G for these data, with the estimated values of the upper limit shown as vertical dashed gray lines:

Figure of ECDF

whuber
quelle
This is a very clear and helpful answer, so thank you!
user46768