2

如果我已经知道特定的百分位数，我正在尝试了解如何生成正态分布。

一位用户对类似问题（此处的链接）给出了非常全面的答案，但是当我尝试使用现有数据对其进行测试时，差异太大了。

我是怎么做到的：

x <- c(5,8,11)
PercRank <- c(2.1, 51.1, 98.8)

例如 PercRank = 2.1 表示 2.1% 的数据的值/分数 <= 5（x 的第一个值）。同样，PercRank = 51.1 表示 51.1% 的数据的值/分数 <= 8。

我按照这个链接中的方法。这是我的代码：

cum.p <- c(2.1, 51.1, 98.8)/100
prob <- c( cum.p[1], diff(cum.p), .01)
x <- c(5,8,11)

freq <- 1000 # final output size that we want

# Extreme values beyond x (to sample)
init <- -(abs(min(x)) + 1) 
fin  <- abs(max(x)) + 1

ival <- c(init, x, fin) # generate the sequence to take pairs from
len <- 100 # sequence of each pair

s <- sapply(2:length(ival), function(i) {
  seq(ival[i-1], ival[i], length.out=len)
})
# sample from s, total of 10000 values with probabilities calculated above
out <- sample(s, freq, prob=rep(prob, each=len), replace = T)

quantile(out, cum.p) 
# 2% 51.1% 98.8% 
# 5     8    11 

c(mean(out), sd(out))
# [1] 7.834401 2.214227

所有这些都来自评论（链接），到目前为止一切都很好。然后我尝试检查生成的正态分布与我的拟合值的配合情况：

data.frame(sort(rnorm(1000, mean=mean(out), sd=sd(out))))
...
# 988                                          13.000904
# 989                                          13.028881
# 990                                          13.076649
...
# 1000                                         14.567080

我很担心，因为第 988 个值（例如，1000 个样本的 98.8%）是13.000904，而我为 98.8% 的百分位数拟合的值是11.0。

我多次重新生成分布，方差始终大于需要的值。

我难住了。如果有人能告诉我一种使方差更准确的方法，我将不胜感激。或者，这是不可避免的吗？

（我第一次在这里发帖，如果我违反了规则，我深表歉意——如果需要，我可以说得更清楚。）

4

1 回答 1

1

为什么不将其视为优化问题？

x <- c(5,8,11)
PercRank <- c(2.1, 51.1, 98.8)

fun <- function(par, pq) {
  sum((log(pq[,1]/100)-pnorm(pq[,2], mean=par[1], sd=par[2], log.p=TRUE))^2)
}

par.estimates <- optim(c(0,1), fn=fun, pq=cbind(PercRank, x))

pnorm(11, par.estimates[[1]][1], par.estimates[[1]][2])
#[1] 0.9816948

结果似乎合理，但与 q=11 的预期值存在一些差异。但是，我怀疑这是您的数据的问题（例如，由于四舍五入），因为以下方法效果很好：

PercRank <- pnorm(x, 8, 2)*100
par.estimates <- optim(c(0,1), fn=fun, pq=cbind(PercRank, x))
par.estimates[[1]]
#[1] 7.999774 1.999953

当然，对于这个特定问题，可能会有更好的优化器。

于 2013-11-15T09:05:55.260 回答