Einbindung einer empirischen CDF

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 .$0$$b$$\int xh(x)dx = k$$b$

Nach Teilen integrierend, kann ich die Gleichung als umschreiben

, wobei das Integral von 0 bis b ist ------- (1)$bG(b) - \int_0^b G(x)dx = k$$0$$b$

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$

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 vorliegt$x_1 \le x_2 \le \cdots \le x_n$$G$$x_i$$\gamma$$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$x_i$$\gamma$$t \ge 1$$x_i$$\gamma$$[\alpha, \beta]$$\gamma$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 γ .$G$$k/n$$\gamma$$(k+t)/n$$\gamma$

Betrachten Sie den Beitrag zu aus dem Intervall [ α $\int_0^b x h(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 [ α , β ] :$[\alpha,\beta]$$h$$t/n$$\gamma$$[\alpha,\beta]$

$$\int_\alpha^\beta x h(x) dx = \left(x G(x)\right)\vert_\alpha^\beta - \int_\alpha^\beta G(x) dx = \left(\beta G(\beta) - \alpha G(\alpha)\right) -\int_\alpha^\beta 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 :$\gamma$$G$

$$\int_\alpha^\beta G(x)dx = \int_\alpha^\gamma G(\alpha) dx + \int_\gamma^\beta G(\beta) dx = (\gamma-\alpha)G(\alpha) + (\beta-\gamma)G(\beta).$$

Ersetzt man dies in das Vorstehende und erinnert man sich an ergibt sich G ( β ) = ( k + t ) /$G(\alpha)=k/n, G(\beta)=(k+t)/n$

$$\int_\alpha^\beta x h(x) dx = \left(\beta G(\beta) - \alpha G(\alpha)\right) - \left((\gamma-\alpha)G(\alpha) + (\beta-\gamma)G(\beta)\right) = \gamma\frac{t}{n}.$$

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

$$\frac{t}{n} = \frac{1}{n} + \cdots + \frac{1}{n}$$

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

$$\int_0^b x h(x) dx = \sum_{i:\, 0 \le x_i \le b} \left(x_i\frac{1}{n}\right) = \frac{1}{n}\sum_{x_i\le b}x_i.$$

$1/n$$[0,b]$$1/n$$1/m$$m$$[0,b]$

$k$$b$$\frac{1}{n}\sum_{x_i\le b}x_i = k.$$k$$j$

$$\frac{1}{n}\sum_{i=1}^{j-1} x_i \le k \lt \frac{1}{n}\sum_{i=1}^j x_i,$$

you will have narrowed $b$ to the interval $[x_{j-1}, x_j)$. 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 = \int_0^b x \exp(-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: