Integrando un CDF empírico

Tengo una distribución empírica . Lo calculo de la siguiente manera$G(x)$

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

Denote , es decir, es el pdf mientras que es el cdf.h $h(x) = dG/dx$$h$$G$

Ahora quiero resolver una ecuación para el límite superior de integración (por ejemplo, ), de modo que el valor esperado de sea ​​algo .x $a$$x$$k$

Es decir, integrando de a , debería tener . Quiero resolver para .b ∫ x h ( x ) d x = k $0$$b$$\int xh(x)dx = k$$b$

Integrando por partes, puedo reescribir la ecuación como

$bG(b) - \int_0^b G(x)dx = k$ , donde la integral es de a ------- (1)$0$$b$

Creo que puedo calcular la integral de la siguiente manera

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

Pero cuando trato de usar esta función con

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

donde diversión es eq (1), obtengo el siguiente error

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

Creo que el problema es que mi función intgrlse evalúa en un valor numérico, mientras uniroot.Allpasa el intervaloc(0,1000)

¿Cómo debo resolver para en esta situación en R?$b$

Deje que los datos ordenados sean . Para comprender el CDF empírico , considere uno de los valores de vamos a llamarlo y suponga que algún número de es menor que y de es igual a . Elija un intervalo en el que, de todos los valores de datos posibles, solo aparezca . Entonces, por definición, dentro de este intervalo tiene el valor constante para números menores que$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$y salta al valor constante para números mayores que .$(k+t)/n$$\gamma$

Considere la contribución a desde el intervalo . Aunque no es una función, es una medida puntual del tamaño en integral se define mediante la integración por partes para convertirla en una integral honesta a la bondad. Hagamos esto durante el intervalo :$\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.$$

El nuevo integrando, aunque es discontinuo en , es integrable. Su valor se encuentra fácilmente al romper el dominio de integración en las partes que preceden y siguen al salto en :$\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).$$

Sustituyendo esto en lo anterior y recordando los rendimientos$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}.$$

En otras palabras, esta integral multiplica la ubicación (a lo largo del eje ) de cada salto por el tamaño de ese salto. El tamaño del salto es$X$

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

con un término para cada uno de los valores de datos que es igual a . Agregar las contribuciones de todos esos saltos de muestra que$\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.$$

Podríamos llamar a esto una "media parcial", ya que es igual a veces una suma parcial. (Tenga en cuenta que es no una expectativa Puede estar relacionado con la expectativa de una versión de la distribución subyacente que se ha truncado al intervalo. : debe reemplazar el del factor de , donde es el número de valores de datos dentro de .)$1/n$$[0,b]$$1/n$$1/m$$m$$[0,b]$

Dado , desea encontrar para el cualDebido a que las sumas parciales son un conjunto finito de valores, por lo general no hay una solución: tendrá que conformarse con la mejor aproximación, que se puede encontrar por horquillado entre dos medios parciales, si es posible. Es decir, al encontrar tal que$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,$$

habrás reducido al intervalo . No puedes hacer nada mejor que eso usando el ECDF. (Al ajustar una distribución continua al ECDF, puede interpolar para encontrar un valor exacto de , pero su precisión dependerá de la precisión del ajuste).[ x j - 1 , x j ) $b$$[x_{j-1}, x_j)$$b$

---

Rrealiza el cálculo de la suma parcial con cumsumy encuentra dónde cruza cualquier valor especificado utilizando la whichfamilia de búsquedas, como en:

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

El resultado en este ejemplo de datos extraídos de una distribución exponencial es

El límite superior se encuentra entre 0.39 y 0.57

El valor verdadero, resolviendo es . Su cercanía a los resultados informados sugiere que este código es preciso y correcto. (Las simulaciones con conjuntos de datos mucho más grandes siguen respaldando esta conclusión).$0.1 = \int_0^b x \exp(-x)dx,$$0.531812$

Aquí hay una gráfica del CDF empírico para estos datos, con los valores estimados del límite superior mostrados como líneas grises discontinuas verticales:$G$