从 R 中的数据的完全高斯拟合中获取分位数

Question

我一直在纠结R如何计算分位数和数据的正常拟合。我的数据（NDVI 值）遵循截断的正态分布（见图）

我有兴趣从数据和拟合正态分布曲线中获取最低的第 10 个百分位值 (p=0.1)。

以我的理解，因为数据被截断了，所以两者应该是完全不同的：我期望数据的分位数高于正态分布计算的分位数，但事实并非如此。据我对分位数函数的理解，数据中的分位数应该是默认的分位数函数：

q=quantile(y, p=0.1)

而正态分布的分位数是：

qx=quantile(y, p=0.1, type=9)

然而这两个结果在所有情况下都非常接近，这让我想知道 R 适合什么类型的分布来计算分位数（截断的正态分布。？）

我也尝试根据拟合正态曲线计算分位数：

fitted=fitdist(as.numeric(y), "norm", discrete = T)
fit.q=as.numeric(quantile(fitted, p=0.1)[[1]][1])

但没有区别。

所以我的问题是： R 将计算分位数的数据拟合到什么曲线，特别是对于 type=9 ？如何根据完全正态分布（包括下尾）计算分位数？

我不知道如何为此生成可重现的示例，但数据可在 https://dl.dropboxusercontent.com/u/26249349/data.csv

谢谢！

Answer 1

R 在确定分位数时使用数据的经验排序，而不是假设任何特定分布。

截断数据的第 10 个百分位数和数据的正态分布拟合恰好非常接近，尽管第 1 个百分位数有很大不同。例如：

# Load data
df = read.csv("data.csv", header=TRUE, stringsAsFactors=FALSE)

# Fit a normal distribution to the data
df.dist = fitdist(df$x, "norm", discrete = T)

现在让我们得到拟合分布和原始数据的分位数。除了第 10 个百分位之外，我还包括了第 1 个百分位。您可以看到拟合正态分布的第 10 个百分位数仅比数据低一点。但是，拟合正态分布的第一个百分位数 低很多 。

quantile(df.dist, p=c(0.01, 0.1))

Estimated quantiles for each specified probability (non-censored data)
           p=0.01    p=0.1
estimate 1632.829 2459.039

quantile(df$x, p=c(0.01, 0.1))

       1%     10% 
  2064.79 2469.90

quantile(df$x, p=c(0.01, 0.1), type=9)

        1%      10% 
  2064.177 2469.400

您还可以通过直接对数据进行排名以及获取均值和标准差等于 fitdist:

中的拟合值的正态分布的第 1 个和第 10 个百分位数来查看这一点

# 1st and 10th percentiles of data by direct ranking
df$x[order(df$x)][round(c(0.01,0.1)*5780)]
[1] 2064 2469

# 1st and 10th percentiles of fitted distribution 
qnorm(c(0.01,0.1), df.dist$estimate[1], df.dist$estimate[2])
[1] 1632.829 2459.039

让我们绘制原始数据（蓝色）和从拟合正态分布生成的假数据（红色）的直方图。重叠的区域是紫色的。

# Histogram of data (blue)
hist(df$x, xlim=c(0,8000), ylim=c(0,1600), col="#0000FF80")

# Overlay histogram of random draws from fitted normal distribution (red)
set.seed(685)
set.seed(685)
x.fit = rnorm(length(df$x), df.dist$estimate[1], df.dist$estimate[2])
hist(x.fit, add=TRUE, col="#FF000080")

或者我们可以绘制数据（蓝色）的经验累积分布函数（ecdf）和拟合正态分布（红色）的随机抽取。水平灰线标记第 10 个百分位数：

plot(ecdf(df$x), xlim=c(0,8000), col="blue")
lines(ecdf(x.fit), col="red")
abline(0.1,0, col="grey40", lwd=2, lty="11")

既然我已经完成了这些，我想知道您是否期望 fitdist 到 return 如果您的数据真的来自正态分布且未被截断。相反，fitdist return 是一个正态分布，具有手头（截断的）数据的均值和 sd，因此由 fitdist 编辑的分布 return 向右移动到我们可能 "expected" 的地方。

c(mean=mean(df$x), sd=sd(df$x))

     mean        sd 
3472.4708  790.8538

df.dist$estimate

     mean        sd 
3472.4708  790.7853

或者，另一个简单的例子：x 服从均值 ~ 0 和 sd ~ 1 的正态分布。xtrunc 删除所有小于 -1 的值，xtrunc.dist 是输出在 xtrunc 上 fitdist：

set.seed(55)
x = rnorm(6000)
xtrunc = x[x > -1]
xtrunc.dist = fitdist(xtrunc, "norm")

round(cbind(sapply(list(x=x,xtrunc=xtrunc), function(x) c(mean=mean(x),sd=sd(x))),
      xtrunc.dist=xtrunc.dist$estimate),3)

          x xtrunc xtrunc.dist
mean -0.007  0.275       0.275
sd    1.009  0.806       0.806

您可以在下面的 ecdf 图中看到，截断数据和拟合截断数据的正态分布具有大致相同的第 10 个百分位数，而未截断数据的第 10 个百分位数（正如我们所期望的那样）发生了偏移向左。

从 R 中的数据的完全高斯拟合中获取分位数

obtaining quantiles from complete gaussian fit of data in R

r

normal-distribution

quantile