Интеграция эмпирического CDF

У меня есть эмпирическое распределение . Я рассчитываю это следующим образом$G(x)$

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

Я обозначаю , т. - это pdf, а - это cdf.h $h(x) = dG/dx$$h$$G$

Теперь я хочу решить уравнение для верхнего предела интегрирования (скажем, ), чтобы ожидаемое значение составляло .х $a$$x$$k$

То есть, интегрируя от до , я должен иметь . Я хочу решить для .b ∫ x h ( x ) d x = k $0$$b$$\int xh(x)dx = k$$b$

Интегрируя по частям, я могу переписать уравнение как

0 $bG(b) - \int_0^b G(x)dx = k$ , где интеграл от до ------- (1)$0$$b$

Я думаю, что могу вычислить интеграл следующим образом

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

Но когда я пытаюсь использовать эту функцию с

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

где fun is eq (1), я получаю следующую ошибку

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

Я думаю, что проблема заключается в том, что моя функция intgrlоценивается по числовому значению, в то время uniroot.Allкак проходит интервалc(0,1000)

Как я должен решить для в этой ситуации в R?$b$

Пусть отсортированные данные будут . Чтобы понять эмпирический CDF , рассмотрим одно из значений в --let, называемом -, и предположим, что некоторое число в меньше и в равно . Выберите интервал в котором из всех возможных значений данных отображается только . Тогда по определению в этом интервале имеет постоянное значение для чисел, меньших$x_1 \le x_2 \le \cdots \le x_n$$G$$x_i$$\gamma$$k$$x_i$$\gamma$$t \ge 1$$x_i$$\gamma$$[\alpha, \beta]$$\gamma$$G$$k/n$$\gamma$и переходит к постоянному значению для чисел, превышающих .$(k+t)/n$$\gamma$

Рассмотрим вклад в из интервала . Хотя не является функцией - это точечная мера размера в - интеграл определяется посредством интегрирования частями, чтобы преобразовать его в интеграл честности в доброту. Давайте сделаем это за интервал :$\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.$$

Новое подынтегральное выражение, хотя оно и разрывно в , является интегрируемым. Его значение легко найти, разбив область интегрирования на части, предшествующие и следующие за скачком в :$\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).$$

Подставляя это в вышеизложенное и вспоминая дает$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}.$$

Другими словами, этот интеграл умножает местоположение (вдоль оси ) каждого прыжка на размер этого прыжка. Размер прыжка$X$

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

с одним членом для каждого из значений данных, равным . Добавление вкладов от всех таких скачков показывает, что$\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,$$

Вы сузите до интервала . Вы можете сделать не лучше, чем с помощью ECDF. (Подбирая некоторое непрерывное распределение к ECDF, вы можете интерполировать, чтобы найти точное значение , но его точность будет зависеть от точности подбора.)[ x j - 1 , x j ) $b$$[x_{j-1}, x_j)$$b$

---

Rвыполняет вычисление частичной суммы с помощью cumsumи находит, где оно пересекает любое указанное значение, используя whichсемейство поисков, как в:

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

Выходные данные в этом примере данных, извлеченных из экспоненциального распределения:

Верхний предел лежит между 0,39 и 0,57

Истинное значение, решающее составляет . Его близость к сообщенным результатам позволяет предположить, что этот код является точным и правильным. (Моделирование с гораздо большими наборами данных продолжает поддерживать этот вывод).$0.1 = \int_0^b x \exp(-x)dx,$$0.531812$

Вот график эмпирического CDF для этих данных с оценочными значениями верхнего предела, показанными в виде вертикальных пунктирных серых линий:$G$