机器算法验证 - 从矩和分位数近似密度，然后从中采样 - 吾爱随笔录

情况

我需要将 R 代码发送给第三方来为我运行估计（我将无法直接处理数据）。我想在将代码发送给它们之前模拟数据以测试一些估计器。

数据提供者为我提供了以下所有变量的汇总统计数据：前四个矩（均值、方差、偏度、峰度）、四个最大/最小值（最小值/最大值）以及 1%、5%、10% , 25%, 50%, 75%, 90%, 95%, 99% 百分位数。

我相信输出是在 Stata 中获得的，我在下面粘贴一个示例。该示例并非来自实际数据（实际数据有接近 600 万个观测值）。

问题

模拟数据的最佳方法是什么？最初我只是要选择一个分布并从中采样（例如二项式/多项式/（对数）正态/指数/截断正态），但根据提供的信息，我认为可以做得更好，至少对于那些不是二项式。

输入示例

      Percentiles      Smallest
 1%          163             99
 5%          216            111
10%          248            113       Obs               3,170
25%          322            114       Sum of Wgt.       3,170

50%          494                      Mean           1262.359
                        Largest       Std. Dev.      3093.165
75%          984          41584
90%         2450          54413       Variance        9567670
95%         5181          58477       Skewness       10.59025
99%        10826          59349       Kurtosis       157.7004

我目前在做什么

library(rriskDistributions)我已经尝试使用（例如get.lnorm.par()）中的函数来拟合我认为可能合适的特定分布的参数。这有时效果很好，但通常效果不佳。

目前我正在使用样条拟合 CDF，使用样条函数导数获得 PDF，然后从中采样。

一般来说，这些方法都不能很好地工作。我希望有一种通用的方法，并提供一个很好的近似值，而无需我手动观察分布并调查拟合。我知道，鉴于我掌握的数据有限，这可能会要求很多。

## function for spline interpolation
splsample <- function(p, v,
                      size = 1000000,
                      vmin = min(v), vmax = max(v),
                      gridsize = min(3*(vmax-vmin), 1000),
                      step = NULL, plot = FALSE, ...) {
  s <- splinefun(v, p, ...)
  if(is.null(step)) {
    grid <- seq(from = vmin, to = vmax, length.out = gridsize)
  } else {
    grid <- seq(from = vmin, to = vmax, by = step)
  }
  pr <- s(grid, deriv = 1)
  pr[pr < 0] <- 0
  if (plot == TRUE) {
    plot(grid, pr)
  }
  bs <- sample(grid, p = pr, size = size, replace = TRUE)
  return(bs)
}

## input
percentiles <- c(0.01, 0.05, 0.10, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99)
values <- c(163, 216, 248, 322, 494, 984, 2450, 5181, 10826)

## spline approximation of pdf
x <- splsample(percentiles, values, plot = TRUE)
summary(x)
mean(x)
var(x)

## alternative: fitting a truncated normal
library("rriskDistributions")
library("msm")
dpar <- get.tnorm.par(p = percentiles, q = values)
x <- rtnorm(10000, mean = dpar["mean"], sd = dpar["sd"],
            lower = dpar["lower"], upper = dpar["upper"])
x[x < 0] <- 0
summary(x)
mean(x)
var(x)