R - 具有预定义的最小值、最大值、平均值和 sd 值的随机分布

Question

我想生成一个随机分布，例如 10,000 个具有预定义的最小值、最大值、平均值和 sd 值的数字。我已经按照这个链接setting upper and lower limits in rnorm 来获得具有固定最小值和最大值的随机分布。但是，这样做时，平均值会发生变化。

例如，

#Function to generate values between a lower limit and an upper limit.
mysamp <- function(n, m, s, lwr, upr, nnorm) {
set.seed(1)
samp <- rnorm(nnorm, m, s)
samp <- samp[samp >= lwr & samp <= upr]
if (length(samp) >= n) {
return(sample(samp, n))
}  
stop(simpleError("Not enough values to sample from. Try increasing nnorm."))
} 
Account_Value <- mysamp(n=10000, m=1250000, s=4500000, lwr=50000, upr=5000000, nnorm=1000000)
summary(Account_Value)

# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 50060 1231000 2334000 2410000 3582000 5000000
#Note - though min and max values are good, mean value is very skewed for an obvious reason.
# sd(Account_Value) # 1397349

我不确定我们是否可以生成满足所有条件的随机正态分布。如果还有其他可以满足所有条件的随机分布，也请分享一下。

期待您的意见。

-谢谢。

Answer 1

您可以使用beta distribution 的广义形式，称为Pearson type I distribution。标准 beta 分布在区间 (0,1) 上定义，但您可以对标准 beta 分布变量进行线性变换以获得任何 (min, max) 之间的值。 this question on CrossValidated 的答案解释了如何在一定的约束条件下使用其均值和方差参数化 beta 分布。

虽然可以用所需的最小值、最大值、平均值和 sd 来制定截断正态分布和广义 beta 分布，但两种分布的形状会非常不同。这是因为截断正态分布在其支持区间的端点处具有正概率密度，而在广义 beta 分布中，密度将始终在端点处平滑地下降到零。哪种形状更可取将取决于您的预期应用。

这是 R 中的一个实现，用于生成具有均值、方差、最小值和最大值参数化的广义 beta 分布观测值。

rgbeta <- function(n, mean, var, min = 0, max = 1)
{
  dmin <- mean - min
  dmax <- max - mean

  if (dmin <= 0 || dmax <= 0)
  {
    stop(paste("mean must be between min =", min, "and max =", max)) 
  }

  if (var >= dmin * dmax)
  {
    stop(paste("var must be less than (mean - min) * (max - mean) =", dmin * dmax))
  }

  # mean and variance of the standard beta distributed variable
  mx <- (mean - min) / (max - min)
  vx <- var / (max - min)^2

  # find the corresponding alpha-beta parameterization
  a <- ((1 - mx) / vx - 1 / mx) * mx^2
  b <- a * (1 / mx - 1)

  # generate standard beta observations and transform
  x <- rbeta(n, a, b)
  y <- (max - min) * x + min

  return(y)
}

set.seed(1)

n <- 10000
y <- rgbeta(n, mean = 1, var = 4, min = -4, max = 5)

sapply(list(mean, sd, min, max), function(f) f(y))
#    [1]  0.9921269  2.0154131 -3.8653859  4.9838290

Answer 2

讨论：

嗨。这是一个非常有趣的问题。它需要相当大的努力才能得到妥善解决，而且并非总能找到解决方案。

第一件事是，当您截断分布（为其设置最小值和最大值）时，标准偏差是有限的（最大值取决于最小值和最大值）。如果你想要它的价值太大 - 你不能得到它。

第二个限制限制的意思。很明显，如果你想要低于最小值和高于最大值的平均值是行不通的，但你可能想要一些太接近极限的东西，但仍然无法满足。

第三个限制限制了这个参数的一个组合。我不确定它是如何工作的，但我很确定并非所有组合都可以满足。

但是有一些组合可能有效并且可能被发现。

解决办法：

问题是：参数是什么：mean 和 sd 具有定义限制的截断（切割）分布 a 和 b，所以最终平均值将等于 desired_mean 和标准差将等于desired_sd。

重要的是参数值：mean 和sd 在在截断之前使用。这就是为什么最终均值和偏差不同的原因。

以下是使用函数optim() 解决问题的代码。它可能不是解决此问题的最佳解决方案，但通常可以：

require(truncnorm)

eval_function <- function(mean_sd){
    mean <- mean_sd[1]
    sd <- mean_sd[2]
    sample <- rtruncnorm(n = n, a = a, b = b, mean = mean, sd = sd)
    mean_diff <-abs((desired_mean - mean(sample))/desired_mean)
    sd_diff <- abs((desired_sd - sd(sample))/desired_sd)
    mean_diff + sd_diff
}

n = 1000
a <- 1
b <- 6
desired_mean <- 3
desired_sd <- 1

set.seed(1)
o <- optim(c(desired_mean, desired_sd), eval_function)

new_n <- 10000
your_sample <- rtruncnorm(n = new_n, a = a, b = b, mean = o$par[1], sd = o$par[2])
mean(your_sample)
sd(your_sample)
min(your_sample)
max(your_sample)
eval_function(c(o$par[1], o$par[2]))

如果该问题还有其他解决方案，我非常感兴趣，所以如果您找到其他答案，请发布它们。

编辑：

@Mikko Marttila：感谢您的评论和链接：Wikipedia 我实现了计算截断分布的均值和标准差的公式。现在解决方案更加优雅，它应该非常准确地计算所需分布的均值和标准差（如果存在）。它的工作速度也快得多。

我实现了eval_function2，它应该在optim()函数中使用，而不是之前的：

eval_function2 <- function(mean_sd){
    mean <- mean_sd[1]
    sd <- mean_sd[2]

    alpha <- (a - mean)/sd
    betta <- (b - mean)/sd

    trunc_mean <- mean + sd * (dnorm(alpha, 0, 1) - dnorm(betta, 0, 1)) / 
                  (pnorm(betta, 0, 1) - pnorm(alpha, 0, 1))

    trunc_var <- (sd ^ 2) * 
                 (1 + 
                  (alpha * dnorm(alpha, 0, 1) - betta * dnorm(betta, 0, 1))/
                  (pnorm(betta, 0, 1) - pnorm(alpha, 0, 1)) -
                 (dnorm(alpha, 0, 1) - dnorm(betta, 0, 1))/
                 (pnorm(betta, 0, 1) - pnorm(alpha, 0, 1)))

    trunc_sd <- trunc_var ^ 0.5

    mean_diff <-abs((desired_mean - trunc_mean)/desired_mean)
    sd_diff <- abs((desired_sd - trunc_sd)/desired_sd)
}