R生成围绕特定平均值的有界随机样本

Question

我已经被这个问题困扰了一段时间，所以我决定写一个问题。

问题：如何生成具有下/上界限并围绕特定均值的随机样本（长度为n）。

观察：分布不需要具体（可以是正常的、测试版等）。

考虑的方法：

• 一种方法是使用rtnorm 函数(package msm) 生成一个在指定范围内具有正态分布的随机数，但它不符合您想要的平均值。

• 我尝试过的第二种方法是这个函数，我在一个我再也找不到的问题中找到了这个函数

  rBootstrap <- function(n, mean, sd, lowerBound, upperBound){
    range <- upperBound - lowerBound
    m <- (mean-lowerBound) / range #mapping mean to 0-1 range
    s <- sd / range #mapping sd to 0-1 range
    a <- (m^2 - m^3 - m*s^2)/s^2 #calculating alpha for rbeta 
    b <- (m-2*m^2+m^3-s^2+m*s^2)/s^2 #calculating beta for rbeta
    data <- rbeta(n,a,b)  #generating data
    data <- lowerBound + data * range #remaping to given bounds
    return(data)
  }
  

  这个函数实际上给出了很好的结果，除非：upperBound > lowerBound + (2* mean - lowerBound)（上限超过了从 lowerBound 到平均值的距离的两倍）。

特别是，我想生成一个长度为 1,800 的随机样本，其值在 50,000 到 250,000 之间，平均值 = 70,000。

Answer 1

您应该使用截断的正态分布，但应该重新校准 mean。如果你看rtnorm中的mean，里面写得很清楚：mean是截断前原始正态分布的均值。

如果您希望 OBSERVABLE 均值等于所需值，只需使用来自Truncated Normal 的公式：

mu = E + sigma*(f(b) - f(a))/(F(b) - F(a))

这里 E 是您想要的平均值（在您的情况下为 70,000），f(x) 是高斯密度，F(x) 是累积函数，a 和 b 是区间边界（居中和缩放)。

a = (LB - mu)/sigma
b = (RB - mu)/sigma

计算出mu 后，将其作为mean 参数传递给rtnorm。

注意：您可能希望对 sigma 进行类似的练习 - 进入 rtnorm 的内容不是您在采样中要观察到的内容，请再次查看 wiki 参考

更新

好的，我自己去写代码，虽然现在第一次剪切是在 Python 中完成的（查看 R）。问题是，对于给定的可观察平均值mu 在f(a)、f(b)、F(a) 和F(b) 中，这将问题转化为对非线性方程根的搜索。不过是可以解决的，请查看code。请注意，它几乎遵循 wiki 符号。

例如对于您的参数和 sigma=12,000，我得到了

Found mu = 68430.372119287 for the desired mean 70000.0 and sigma 12000.0
Sampled 100000 truncated gaussians and got observed mean = 70023.15990337673

对于你的参数和 sigma=24,000，我得到了

Found mu = 52275.475000378945 for the desired mean 70000.0 and sigma 24000.0
Sampled 100000 truncated gaussians and got observed mean = 69922.16000288539

因此，mu 非常接近大 sigma 的左边界，这是预期的行为，但观察到的平均值保持在接近 70,000，这是您想要的。

更新二

这是 R 代码，也在 github 存储库中

require(rootSolve)
require(msm)

phi <- function(z) {
    dnorm(z)
}

Phi <- function(z) {
    pnorm(z)
}

Mean <- function(mu, sigma, a, b) {
    alfa <-  (a - mu) / sigma
    beta <-  (b - mu) / sigma

    Z <-  Phi(beta) - Phi(alfa)

    mu + sigma*(phi(alfa) - phi(beta))/Z
}

f <- function(mu, mean, sigma, a, b) {
    mean - Mean(mu, sigma, a, b)
}

a <-  50000.0
b <-  250000.0
mean  <- 70000.0
sigma <- 24000.0

# find mu for desired mean
q <- uniroot(f, c(a, b), mean, sigma, a, b)
mu <- q$root

print(sprintf("Found mu = %f for the desired mean %f and sigma %f", mu, mean, sigma))

# sampling test
set.seed(32345)
N = 100000
r <- rtnorm(N, mean=mu, sd=sigma, lower=a, upper=b)

print(sprintf("Sampled %d truncated gaussians and got observed mean = %f", N, mean(r)))