Integracja empirycznego CDF

Mam rozkład empiryczny . Obliczam to w następujący sposób$G(x)$

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

Mam na myśli , tzn. to pdf, a to cdf.h $h(x) = dG/dx$$h$$G$

Chcę teraz rozwiązać równanie dla górnej granicy całkowania (powiedzmy ), tak że oczekiwana wartość wynosi jakieś .x $a$$x$$k$

To znaczy, całkując od do , powinienem mieć . Chcę rozwiązać dla .b ∫ x h ( x ) d x = k $0$$b$$\int xh(x)dx = k$$b$

Całkując przez części, mogę przepisać równanie jako

0 $bG(b) - \int_0^b G(x)dx = k$ , gdzie całka wynosi od do ------- (1)$0$$b$

Myślę, że mogę obliczyć całkę w następujący sposób

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

Ale kiedy próbuję użyć tej funkcji z

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

gdzie fun to eq (1), pojawia się następujący błąd

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

Myślę, że problem polega na tym, że moja funkcja intgrljest oceniana na wartość liczbową, podczas gdy uniroot.Allmija przedziałc(0,1000)

Jak mam rozwiązać problem w tej sytuacji w R?$b$

Niech posortowane dane będą . Aby zrozumieć empiryczną CDF , rozważ jedną z wartości - nazwij ją załóż , że pewna liczba z jest mniejsza niż a z jest równa . Wybierz przedział w którym ze wszystkich możliwych wartości danych pojawia się tylko . Następnie, z definicji, w tym przedziale ma stałą wartość dla liczb mniejszych niż G x i γ k x i γ t ≥ 1 x i$x_1 \le x_2 \le \cdots \le x_n$$G$$x_i$$\gamma$$k$$x_i$$\gamma$$t \ge 1$$x_i$$\gamma$γ G k / n γ ( k + t ) / n $[\alpha, \beta]$$\gamma$$G$$k/n$$\gamma$i skacze do stałej wartości dla liczb większych niż .$(k+t)/n$$\gamma$

Rozważ wkład do z przedziału . Chociaż nie jest funkcją - jest punktową miarą wielkości w - całka jest definiowana za pomocą całkowania przez części, aby przekształcić ją w całkę uczciwą na dobroć. Zróbmy to w przedziale :[ α , β ] h t / n γ [ α , β $\int_0^b x h(x) dx$$[\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.$$

Nowy integrand, chociaż jest nieciągły w , jest całkowalny. Jego wartość można łatwo znaleźć, dzieląc domenę integracji na części poprzedzające i następujące po skoku w :$\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).$$

Podstawiając to do powyższego i przywołując wydajności$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}.$$

Innymi słowy, ta całka zwielokrotnia lokalizację (wzdłuż osi ) każdego skoku przez jego wielkość. Rozmiar skoku to$X$

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

z jednym terminem dla każdej wartości danych równej . Dodanie wkładów ze wszystkich takich skoków pokazuje to$\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.$$

Możemy to nazwać „średnią cząstkową”, widząc, że jest równa razy suma częściowa. (Należy pamiętać, że to nie oczekiwanie może być związany z oczekiwaniem wersji rozkładu bazowego, który został obcięty do przedziału. : należy zastąpić współczynnik przez gdzie to liczba wartości danych w .)[ 0 , b $1/n$$[0,b]$$1/n$$1/m$$m$$[0,b]$

Biorąc pod uwagę , chcesz znaleźć dla któregoPonieważ sumy cząstkowe są skończonym zestawem wartości, zwykle nie ma rozwiązania: trzeba się zadowolić najlepszym przybliżeniem, które można znaleźć , jeśli to możliwe, nawiasując pomiędzy dwoma średnimi cząstkowymi. To znaczy, po znalezieniu takiego$k$$b$k$\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,$$

zawęzisz do przedziału . Nie możesz zrobić nic lepszego niż korzystając z ECDF. (Dopasowując ciągły rozkład do ECDF, można interpolować w celu znalezienia dokładnej wartości , ale jej dokładność będzie zależeć od dokładności dopasowania.)[ x j - 1 , x j ) $b$$[x_{j-1}, x_j)$$b$

---

Rwykonuje obliczenie sumy częściowej za pomocą cumsumi wyszukuje, gdzie przecina dowolną określoną wartość, używając whichrodziny wyszukiwań, jak w:

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])

Dane wyjściowe w tym przykładzie danych pobranych z rozkładu wykładniczego to

Górna granica wynosi od 0,39 do 0,57

Prawdziwa wartość rozwiązująca wynosi . Jego bliskość do zgłaszanych wyników sugeruje, że ten kod jest dokładny i poprawny. (Symulacje ze znacznie większymi zestawami danych nadal potwierdzają ten wniosek).$0.1 = \int_0^b x \exp(-x)dx,$$0.531812$

Oto wykres empirycznego CDF dla tych danych, z szacowanymi wartościami górnej granicy pokazanymi jako pionowe przerywane szare linie:$G$