如何计算R中内核估计的Kullback-leibler散度答案

【问题标题】：How to calculate Kullback-leiber divergence of Kernel estimation in R如何计算R中内核估计的Kullback-leibler散度
【发布时间】：2020-05-12 22:13:03
【问题描述】：

我使用核估计来获得非参数概率密度函数。然后，我想使用 Kullback-leiber 散度比较两个连续变量的内核分布之间的尾部“距离”。我试过以下代码：

kl_l <- function(x,y) {
    integrand <- function(x,y) {

            f.x <- fitted(density(x, bw="nrd0"))
            f.y <- fitted(density(y, bw="nrd0"))

            return((log(f.x)-log(f.y))*f.x) 
    }
    return(integrate(integrand, lower=-Inf,upper=quantile(density(x,  bw="nrd0"),0.25))$value)
    #the Kullback-leiber equation
}

当我为 a, b = 19 个连续变量运行 kl_l(a,b) 时，它会返回警告

 Error in density(y, bw = "nrd0") : argument "y" is missing, with no default

有什么方法可以计算出来吗？

（如果有人想查看实际方程式：https://www.bankofengland.co.uk/-/media/boe/files/working-paper/2019/attention-to-the-tails-global-financial-conditions-and-exchange-rate-risks.pdf第 13 页。）

【问题讨论】：

你可以看看LaplacesDemon包中的KLD函数
感谢您的建议。但是，我只想研究分布的尾部，这意味着我需要一个函数来指定积分的边界（例如，从分布的第 25 个分位数的值到 -inf 的积分）。检查方程式：bankofengland.co.uk/-/media/boe/files/working-paper/2019/… 第 13 页
R 中还有其他可用的这种分歧的实现，但也许您已经检查过它们。在上面的代码中，我认为问题在于integrate 与 1 变量的函数一起使用，或者至少与积分在第一个变量上但第二个变量取常数值的函数一起使用。在任何情况下，您都没有传递 y 参数。
非常感谢您的帮助。我去看看。

标签： r

【解决方案1】：

简而言之，我认为您只需将f.x 和f.y 移到被积函数之外（并可能将fitted 替换为approxfun）：

kl_l <- function(x, y) {
    f.x <- approxfun(density(x, bw = "nrd0"))
    f.y <- approxfun(density(y, bw = "nrd0"))
    integrand <- function(z) {
        return((log(f.x(z)) - log(f.y(z))) * f.x(z)) 
    }
    return(integrate(integrand, lower = -Inf, upper = quantile(density(x, bw="nrd0"), 0.25))$value)
    #the Kullback-leiber equation
}

扩大一点：

查看您引用的论文，您似乎需要首先创建两个拟合分布f 和g。因此，如果您的变量 a 包含在全球金融条件下增加 1 个标准差的观察值，并且 b 包含在平均全球金融条件下的观察值，您可以像示例中那样创建两个函数：

f <- approxfun(density(a))
g <- approxfun(density(b))

然后定义被积函数：

integrand <- function(x) log(f(x) / g(x)) * f(x)

上限：

upper <- quantile(density(b, bw = "nrd0"), 0.25)

最后在指定范围内对x 进行积分。注意在数值计算中x的每个值都必须同时进入f和g；在你的函数kl_l 中，x 和 y 分别进入被积函数，我认为这是不正确的；在任何情况下，integrate 只会对第一个变量进行操作。

integrate(integrand, lower = -Inf, upper = upper)$value

要检查的一件事是，approxfun 返回的值超出了密度中指定的范围，这可能会扰乱您的操作，因此您需要针对这些值进行调整（如果您希望密度为例如，归零）。

【讨论】：