有没有办法创建一个符合以下参数的假数据集:N、mean、sd、min 和 max?
我想创建一个包含 187 个整数尺度分数的样本,这些分数的平均值为 67,标准差为 17,观测值在 [30, 210] 范围内。我正在尝试演示有关统计功效的概念类(class),并且我想创建具有看起来像已发布结果的分布的数据。此示例中的量表分数是 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
当然,平均值和标准差并不精确,向量也不由整数组成。
还有其他选择吗?
最佳答案
无模板的整数优化
由于您希望获得精确的均值、标准差、最小值和最大值,因此我的首选不是随机数生成,因为您的样本不太可能与您的分布的均值和标准差完全匹配从绘制。相反,我会采用整数优化方法。您可以定义变量 x_i
为整数i
的次数出现在您的示例中。您将定义决策变量 x_30
, x_31
, ..., x_210
并添加确保满足所有条件的约束:
- 187 个样本:这可以通过约束
x_30 + x_31 + ... + x_210 = 187
进行编码 - 67 的平均值:这可以通过约束
30*x_30 + 31*x_31 + ... + 210*x_210 = 187 * 67
进行编码 - 变量的逻辑约束:变量必须采用非负整数值
- “看起来像真实数据” 这显然是一个定义不明确的概念,但我们可以要求相邻数字的频率之差不超过 1。这是表格
x_30 - x_31 <= 1
,x_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 秒内完成了计算。
使用模板进行整数优化
获得类似于“真实数据”的另一种方法是定义一个起始连续分布(例如,您在原始帖子中包含的 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 秒),并且仍然返回完全满足所需均值和标准差的分布。此外,给定来自连续分布的足够高质量的样本,这可用于获取具有整数值并满足所需统计特性的不同形状的分布。
关于r - 创建一个符合以下参数的假数据集 : N, mean、sd、min 和 max,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/32792824/