R生成有界随机样本周围特定平均值
我一直坚持这一段时间,所以我决定写一个问题。R生成有界随机样本周围特定平均值
问题:如何使用结合的下部/上部和角落找寻一个特定意味着生成(的lenght n)的随机样本。
观察:分布不需要具体(可能是正常的,测试版等)。
Aproaches认为:
- 一种形式给出是使用
rtnorm
功能(package msm)与指定范围内的正态分布产生一个随机数,但它不会把你想要的平均值。 -
我已经尝试了第二的aproach是这个功能,我在一个问题中我找不到了
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>下界+(2 *平均 - lowerBound)(上限超过了从lowerBound到mean的距离的两倍)。
特别是,我想生成一个长度为1,800的随机样本,值在50,000到250,000之间,平均值= 70,000。
您应该使用截断的正态分布,但mean
应该重新校准。如果您看rtnorm
中的mean
,则明确指出:mean
是截断前原始正态分布的均值。
如果你想可观测平均值等于期望值,只是用公式从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
之后,它向下传递到rtnorm为mean
参数。
注意:您可能想要做类似的工作与sigma
- 这是怎么回事成rtnorm是不是你打算在抽样观察,再次看到维基参考
UPDATE
好东西,就到了自己编码,尽管第一次剪切是在Python中完成的(现在正在查看R)。问题在于,对于给定的可观察平均值mu
在f(a)
,f(b)
,F(a)
和F(b)
中,其将问题转换为搜索非线性方程的根。但它是可以解决的,请检查code。请注意,它遵循几乎维基表示法。
例如,对于您的参数和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 = 24000我
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,这是你想要的。
UPDATE II
这里是[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)))
谢谢,在这种情况下,f(a)将是我想要的下限?如果F(a)= 0,反之f(b)和F(b)= 1? –
@AlfredoLozano我已经更新了wrt'a'和'b'。不,如果你看wiki的话,'\ phi(x)'是纯高斯的,'\ Phi(x)'是高斯的累积(误差函数的变化),所以F(a)不是0,F b)不是1 –
@AlfredoLozano您必须对西格玛有一定的价值 - 而且,重要的是,如果您希望它成为公式sigma或OBSERVABLE sigma。 –
你想从你的生成随机样本分布是什么?此链接可能有所帮助:http://r.789695.n4.nabble.com/how-to-generate-a-normal-distribution-with-mean-1-min-0-2-max-0-8-td3481450 .html – Chrisss
谢谢@Chrisss,我发现我并不是在寻找一个特定的发行版,尽管我所做的所有研究都是以普通版和测试版为导向的,但我相信这两者中的一个可以通过观察它们的密度函数形状。 –
顺便说一句,你想要什么西格玛?同样,你想要公式西格玛还是可观察西格玛?我在两个小时内飞行,但只要我回来,我会尝试写一些** R ** ... –