估计尾部密度单调下降的 PDF

Question

tldr：我正在从模拟数据中对 PDF 进行数值估计，我需要密度在 'main' 密度区域之外单调递减（因为 x-> 无穷大).我所得到的密度接近于零，但不会单调减少。

详细问题

我正在估计一个模拟的最大似然模型，这需要我在某个（观察到的）值 x 。目标是最大化这些密度的对数似然，这要求它们没有虚假的局部最大值。

因为我没有分析似然函数，所以我通过从某个已知分布函数中提取随机分量来数值模拟随机变量，并对它应用一些非线性变换。我将此模拟的结果保存在名为 simulated_stats.

的数据集中

然后我使用 density() 来近似 PDF 并使用 approxfun() 来评估 x:

处的 PDF

#some example simulation
Simulated_stats_ <- runif(n=500, 10,15)+ rnorm(n=500,mean = 15,sd = 3)
#approximation for x
approxfun(density(simulated_stats))(x)

这在模拟 simulated_stats 范围内效果很好，见图： Example PDF。问题是我需要能够评估远离模拟数据范围的 PDF。

所以在上图中，我需要在 x=50:

处评估 PDF

approxfun(density(simulated_stats))(50)
> [1] NA

所以我在密度函数中使用了 from 和 to 参数，它们正确地近似于 0 尾部，这样

approxfun(
 density(Simulated_stats, from = 0, to = max(Simulated_stats)*10)
)(50)
[1] 1.924343e-18

这很好，在一种情况下 - 我需要密度在距离 x 范围越远的地方变为零。也就是说，如果我在 x=51 处进行评估，结果必须严格更小。（否则，我的估计器可能会找到远离 'true' 区域的局部最大值，因为似然函数在远离 'main' 密度质量（即外推区域）的地方不是单调的。

为了测试这一点，我以固定的时间间隔评估了近似的 PDF，记录了日志并绘制了图表。结果令人沮丧：远离主密度质量的概率'jumps'上下。总是非常接近于零，但不是单调递减。

    a <- sapply(X = seq(from = 0, to = 100, by = 0.5), FUN = function(x){approxfun(
      density(Simulated_stats_,from = 0, to = max(Simulated_stats_)*10)
      )(x)})
    aa <- cbind( seq(from = 0, to = 100, by = 0.5), a)
    plot(aa[,1],log(aa[,2]))

结果： Non-monotonic log density far from density mass

我的问题

发生这种情况是因为 density() 中的核估计还是 approxfun() 中的不准确？（或其他？）

我可以使用哪些替代方法来提供远离模拟密度质量的单调下降 PDF？

或者 - 我如何手动更改近似的 PDF 以单调下降我离密度质量越远？我很乐意坚持一些趋于零的线性趋势...

谢谢！

Answer 1

一种可能性是使用 beta 回归模型估计 CDF；然后可以使用该模型的导数的数值估计来估计任何点的 pdf。这是我的想法的一个例子。我不确定它是否对你有帮助。

导入库

library(mgcv)
library(data.table)
library(ggplot2)

生成您的数据

set.seed(123)
Simulated_stats_ <- runif(n=5000, 10,15)+ rnorm(n=500,mean = 15,sd = 3)

使用 gam beta 回归模型估计 CDF 的函数

get_mod <- function(ss,p = seq(0.02, 0.98, 0.02)) {
  qp = quantile(ss, probs=p)
  betamod = mgcv::gam(p~s(qp, bs="cs"), family=mgcv::betar())
  return(betamod)
}

betamod <- get_mod(Simulated_stats_)

给定估计 CDF

val

est_pdf <- function(val, betamod, tol=0.001) {
  xvals  = c(val,val+tol)
  yvals = predict(betamod,newdata=data.frame(qp = xvals), type="response")
  as.numeric((yvals[1] - yvals[2])/(xvals[1] - xvals[2]))
}

让我们检查是否单调递增到 Simulated_stats

test_x = seq(0,min(Simulated_stats_), length.out=1000)
pdf = sapply(test_x, est_pdf, betamod=betamod)
all(pdf == cummax(pdf))

[1] TRUE

让我们检查是否在 Simulated_stats

test_x = seq(max(Simulated_stats_), 60, length.out=1000)
pdf = sapply(test_x, est_pdf, betamod=betamod)
all(pdf == cummin(pdf))

[1] TRUE

其他想法 3/5/22

正如评论中所讨论的那样，使用 betamod 进行预测可能会减慢估算器的速度。虽然这可以在很大程度上通过直接编写自己的预测函数来解决，但还有另一种可能的捷径。

在 X 的范围内从 betamod 生成估计值，包括极值

k <- sapply(seq(0,max(Simulated_stats_)*10, length.out=5000), est_pdf, betamod=betamod)

使用您最初使用的上述方法，即对密度进行线性插值，而不是对密度结果进行此操作，而是对 k 进行操作（即对来自 beta 的上述估计进行操作）型号)

lin_int = approxfun(x=seq(0,max(Simulated_stats_)*10, length.out=5000),y=k)

您可以在估算器中使用lin_int() 函数进行预测，速度很快。请注意，对于给定的 x

c(est_pdf(38,betamod), lin_int(38))
[1] 0.001245894 0.001245968

而且速度很快

microbenchmark::microbenchmark(
  list = alist("betamod" = est_pdf(38, betamod),"lin_int" = lint(38)),times=100
)

Unit: microseconds
    expr    min      lq     mean  median      uq    max neval
 betamod 1157.0 1170.20 1223.304 1188.25 1211.05 2799.8   100
 lin_int    1.7    2.25    3.503    4.35    4.50   10.5   100

最后，让我们检查一下您之前绘制的相同图，但使用 lin_int() 而不是 approxfun(density(....))

a <- sapply(X = seq(from = 0, to = 100, by = 0.5), lin_int)
aa <- cbind( seq(from = 0, to = 100, by = 0.5), a)
plot(aa[,1],log(aa[,2]))

估计尾部密度单调下降的 PDF

Estimating PDF with monotonically declining density at tails

r

probability

kernel-density

density-plot

log-likelihood

其他想法 3/5/22