интегрируя квадрат плотности вероятности?

Question

Предположим, у меня есть

set.seed(2020)    # make the results reproducible
a <- rnorm(100, 0, 1)

Моя плотность вероятности оценивается через ядро оценки плотности (гауссовой) в R , используя R встроенный в функции density. Вопрос в том, как интегрировать квадрат оцененной плотности. Неважно, между какими значениями, допустим между -Infи +Inf. Я пробовал следующее:

f <- approxfun(density(a)$x, density(a)$y)
integrate (f*f, min(density(a)$x), max(density(a)$x))

Answer 1

Здесь есть пара проблем. Во- первых у вас есть xи yокруглить неправильный путь в approxfun. Во-вторых, нельзя умножать имена функций вместе. Вам нужно указать новую функцию, которая даст вам квадрат исходной функции:

set.seed(2020)   
a  <- rnorm(100, 0, 1)
f  <- approxfun(density(a)$x, density(a)$y)
f2 <- function(v) ifelse(is.na(f(v)), 0, f(v)^2)

integrate (f2, -Inf, Inf)
#> 0.2591153 with absolute error < 0.00011

Мы также можем построить исходную функцию плотности и функцию плотности в квадрате:

curve(f, -3, 3)
curve(f2, -3, 3, add = TRUE, col = "red")

Answer 2

Я думаю, вам следует написать целевую функцию как function(x) f(x)**2, а не f*f, например,

> integrate (function(x) f(x)**2, min(density(a)$x), max(density(a)$x))
0.2331793 with absolute error < 6.6e-06

Answer 3

Вот способ использования package caTools, function trapz. Он вычисляет интеграл по вектору xи его соответствующему изображению, yиспользуя правило трапеций .
Я также включаю функцию, trapzfоснованную на оригинале, чтобы вычислить интеграл с функцией, возвращаемойapproxfun

trapzf <- function(x, FUN) trapz(x, FUN(x))

set.seed(2020)    # make the results reproducible
a <- rnorm(100, 0, 1)

d <- density(a)
f <- approxfun(d$x, d$y)

int1 <- trapz(d$x, d$y^2)
int2 <- trapzf(d$x, function(x) f(x)^2)

int1
#[1] 0.2591226

identical(int1, int2)
#[1] TRUE