创建一个符合以下参数的假数据集:N、平均值、标准差、最小值和最大值

Create a fake dataset that fits the following parameters: N, mean, sd, min, and max

有没有办法创建一个符合以下参数的假数据集:N、平均值、标准差、最小值和最大值?

我想创建一个包含 187 个整数量表分数的样本,其平均值为 67,标准差为 17,观测值在 [30, 210] 范围内。我正在尝试演示有关统计功效的概念课,并且我想创建具有看起来像已发布结果的分布的数据。此示例中的量表分数是 30 个项目的总和,每个项目的范围从 1 到 7。我不需要构成量表分数的各个项目的数据,但这将是一个奖励。

我知道我可以使用 rnorm(),但这些值不是整数,最小值和最大值可能超过我可能的值。

scaleScore <- rnorm(187, mean = 67, sd = 17)

我也知道我可以使用 sample() 来获得保持在这个范围内的整数,但是平均值和标准偏差是不正确的。

scaleScore <- sample(30:210, 187, replace=TRUE)

@Pascal 的提示让我找到了 Runuran 包中的 urnorm()

set.seed(5)
scaleScore <- urnorm(n=187, mean=67, sd=17, lb=30, ub=210)
mean(scaleScore)
# [1] 68.51758
sd(scaleScore)
# [1] 16.38056
min(scaleScore)
# [1] 32.15726
max(scaleScore)
# [1] 107.6758

当然,均值和标准差并不精确,向量也不由整数组成。

还有其他选择吗?

我能够使用蛮力相当接近,即 method="SANN" in optim():

目标 values/constraints:

m0 <- 67
sd0 <- 17
min <- 30
max <- 210
n <- 187

设置初始值:

set.seed(101)
mm <- min:max
x0 <- sample(mm,size=n,replace=TRUE)

Objective 函数(与期望的距离 mean/sd;范围和 N 将受到限制)

objfun <- function(x) {
    (mean(x)-m0)^2+(sd(x)-sd0)^2
}

新参数集的候选分布:随机重采样一个值

candfun <- function(x) {
    x[sample(n,size=1)] <- sample(mm,size=1)
    return(x)
}
objfun(x0)  ## initial badness: 4088.621
set.seed(101)
o1 <- optim(par=x0,fn=objfun,gr=candfun,
      method="SANN",control=list(maxit=1e6))
mean(o1$par) ## 66.978
sd(o1$par) ## 17.22
plot(table(o1$par))

无模板的整数优化

因为你想要一个精确的均值、标准差、最小值和最大值,我的第一选择不会是随机数生成,因为你的样本不太可能与你的分布的均值和标准差完全匹配'从中汲取灵感。相反,我会采用整数优化方法。您可以将变量 x_i 定义为整数 i 在样本中出现的次数。您将定义决策变量 x_30x_31、...、x_210 并添加确保满足所有条件的约束:

  • 187个样本:这可以通过约束x_30 + x_31 + ... + x_210 = 187
  • 进行编码
  • 67的平均值:这可以通过约束30*x_30 + 31*x_31 + ... + 210*x_210 = 187 * 67
  • 进行编码
  • 变量的逻辑约束:变量必须取 non-negative 个整数值
  • "Looks Like Real Data"这显然是一个ill-defined概念,但我们可以要求相邻数字的频率相差不超过1 . 这是形式为 x_30 - x_31 <= 1x_30 - x_31 >= -1 的线性约束,依此类推,适用于每个连续的对。我们还可以要求每个频率不超过某个任意定义的上限(我将使用 10)。

最后,我们希望标准差尽可能接近17,也就是说我们希望方差尽可能接近17^2 = 289。我们可以定义一个变量y来是我们匹配这个方差的接近程度的上限,我们可以最小化 y:

y >= ((30-67)^2 * x_30 + (31-67)^2 * x_31 + ... + (210-67)^2 * x_210) - (289 * (187-1))
y >= -((30-67)^2 * x_30 + (31-67)^2 * x_31 + ... + (210-67)^2 * x_210) + (289 * (187-1))

这是一个非常简单的优化问题,可以使用 lpSolve:

这样的求解器来解决
library(lpSolve)
get.sample <- function(n, avg, stdev, lb, ub) {
  vals <- lb:ub
  nv <- length(vals)
  mod <- lp(direction = "min",
            objective.in = c(rep(0, nv), 1),
            const.mat = rbind(c(rep(1, nv), 0),
                              c(vals, 0),
                              c(-(vals-avg)^2, 1),
                              c((vals-avg)^2, 1),
                              cbind(diag(nv), rep(0, nv)),
                              cbind(diag(nv)-cbind(rep(0, nv), diag(nv)[,-nv]), rep(0, nv)),
                              cbind(diag(nv)-cbind(rep(0, nv), diag(nv)[,-nv]), rep(0, nv))),
            const.dir = c("=", "=", ">=", ">=", rep("<=", nv), rep("<=", nv), rep(">=", nv)),
            const.rhs = c(n, avg*n, -stdev^2 * (n-1), stdev^2 * (n-1), rep(10, nv), rep(1, nv), rep(-1, nv)),
            all.int = TRUE)
  rep(vals, head(mod$solution, -1))
}
samp <- get.sample(187, 67, 17, 30, 210)
summary(samp)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#      30      64      69      67      74     119
sd(samp)
# [1] 17
plot(table(samp))

对于您提供的参数,我们能够在返回所有整数值的同时得到准确的均值和标准差,并且在我的计算机上在 0.4 秒内完成了计算。

使用模板进行整数优化

获得类似于“真实数据”的东西的另一种方法是定义一个起始连续分布(例如,您在原始 post 中包含的 urnorm 函数的结果)并以最能实现均值和标准差目标的方式将值四舍五入为整数。这实际上只引入了两个新的 类 约束:某个值的样本数的上限是可以向上或向下舍入以达到该值的样本数,以及总和的下限两个连续的频率是落在这两个整数之间的连续样本的数量。同样,这很容易用 lpSolve 实现,并且对于 运行:

来说并不是非常低效
library(lpSolve)
get.sample2 <- function(n, avg, stdev, lb, ub, init.dist) {
  vals <- lb:ub
  nv <- length(vals)
  lims <- as.vector(table(factor(c(floor(init.dist), ceiling(init.dist)), vals)))
  floors <- as.vector(table(factor(c(floor(init.dist)), vals)))
  mod <- lp(direction = "min",
            objective.in = c(rep(0, nv), 1),
            const.mat = rbind(c(rep(1, nv), 0),
                              c(vals, 0),
                              c(-(vals-avg)^2, 1),
                              c((vals-avg)^2, 1),
                              cbind(diag(nv), rep(0, nv)),
                              cbind(diag(nv) + cbind(rep(0, nv), diag(nv)[,-nv]), rep(0, nv))),
            const.dir = c("=", "=", ">=", ">=", rep("<=", nv), rep(">=", nv)),
            const.rhs = c(n, avg*n, -stdev^2 * (n-1), stdev^2 * (n-1), lims, floors),
            all.int = TRUE)
  rep(vals, head(mod$solution, -1))
}

library(Runuran)
set.seed(5)
init.dist <- urnorm(n=187, mean=67, sd=17, lb=30, ub=210)
samp2 <- get.sample2(187, 67, 17, 30, 210, init.dist)
summary(samp2)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#      32      57      66      67      77     107
sd(samp2)
# [1] 17
plot(table(samp2))

这种方法甚至更快(不到 0.1 秒)并且仍然 returns 一个完全符合所需均值和标准差的分布。此外,给定来自连续分布的足够高质量的样本,这可用于获得不同形状的分布,这些分布采用整数值并满足所需的统计特性。