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

Question

我想用预定义的最小值、最大值、平均值和 sd 值生成 10,000 个数字的随机分布。我已经按照此链接在 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 分布的广义形式，称为Pearson I 型分布。标准 Beta 分布定义在区间 (0,1) 上，但您可以对标准 Beta 分布变量进行线性变换以获得任意 (min, max) 之间的值。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

下面是使用 function 解决问题的代码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：感谢您的评论和链接：维基百科我实现了计算截断分布的均值和标准差的公式。现在解决方案更加优雅，如果存在，它应该非常准确地计算所需分布的均值和标准差。它的工作速度也快得多。

我实现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)
}