【问题标题】:Generating a normal distribution from known percentiles从已知百分位数生成正态分布
【发布时间】:2013-11-15 08:30:28
【问题描述】:

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

一位用户对类似问题 (link here) 给出了非常全面的答案,但是当我尝试使用现有数据对其进行测试时,差异太大了。

我是怎么做到的:

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

例如,PercRank = 2.1 表示 2.1% 的数据的值/分数

我遵循了link中的方法。这是我的代码:

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

所有这些都来自评论 (linked),到目前为止一切都很好。然后我尝试检查生成的正态分布与我的拟合值的效果:

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。

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

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

(我第一次在这里发帖,如果我违反了规则,我深表歉意 - 如果需要,我可以更清楚地说明。)

【问题讨论】:

    标签: r statistics probability


    【解决方案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
    

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

    【讨论】:

    • 成功了!比以前的方法好多了。非常感谢!
    • 如果我有一个 PercRank 值的数据框,我重复执行此操作的最简单方法是什么?
    • 嗯,您可能可以使用lapplyapply
    • 我想通了!我没有使用lapply/apply,但我不得不对其进行一些重组以使其正常工作。再次感谢您!
    猜你喜欢
    • 1970-01-01
    • 2021-05-05
    • 2020-06-18
    • 1970-01-01
    • 2020-11-13
    • 1970-01-01
    • 1970-01-01
    • 2016-05-02
    • 1970-01-01
    相关资源
    最近更新 更多