可靠地检索分位数函数的反函数
Reliably retrieve the reverse of the quantile function
我看过其他帖子(例如here)关于获取"reverse"分位数——即获取一系列值中某个值对应的百分位数。
但是,对于相同的数据系列 ,答案并没有给我与分位数 相同的值。
我也研究过分位数提供了9种不同的算法来计算百分位数。
所以我的问题是:是否有可靠的方法来获取分位数函数的反函数? ecdf 不采用 "type" 参数,因此似乎无法确保他们使用相同的方法。
可重现的例子:
# Simple data
x = 0:10
pcntile = 0.5
# Get value corresponding to a percentile using quantile
(pcntile_value <- quantile(x, pcntile))
# 50%
# 5 # returns 5 as expected for 50% percentile
# Get percentile corresponding to a value using ecdf function
(pcntile_rev <- ecdf(x)(5))
# [1] 0.5454545 #returns 54.54% as the percentile for the value 5
# Not the same answer as quantile produces
ecdf
给出文档中公式的结果。
x <- 0:10
Fn <- ecdf(x)
现在,对象 Fn
是一个插值阶跃函数。
str(Fn)
#function (v)
# - attr(*, "class")= chr [1:3] "ecdf" "stepfun" "function"
# - attr(*, "call")= language ecdf(x)
并保留原始 x
值和相应的 y
值。
environment(Fn)$x
# [1] 0 1 2 3 4 5 6 7 8 9 10
environment(Fn)$y
# [1] 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 0.54545455
# [7] 0.63636364 0.72727273 0.81818182 0.90909091 1.00000000
后者与文档中所说的用于计算它们的公式的结果完全相同。来自 help('ecdf')
:
For observations x= (x1,x2, ... xn), Fn is the fraction of
observations less or equal to t, i.e.,
Fn(t) = #{xi <= t}/n = 1/n sum(i=1,n) Indicator(xi <= t).
我将使用 seq_along
而不是 1:length(x)
。
seq_along(x)/length(x)
# [1] 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 0.54545455
# [7] 0.63636364 0.72727273 0.81818182 0.90909091 1.00000000
Fn(x)
# [1] 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 0.54545455
# [7] 0.63636364 0.72727273 0.81818182 0.90909091 1.00000000
link 中的答案非常好,但也许有帮助,看看 ecdf
只需 运行 以下代码:
# Simple data
x = 0:10
p0 = 0.5
# Get value corresponding to a percentile using quantile
sapply(c(1:7), function(i) quantile(x, p0, type = i))
# 50% 50% 50% 50% 50% 50% 50%
# 5.0 5.0 5.0 4.5 5.0 5.0 5.0
因此,这不是类型的问题。您可以使用 debug:
进入函数
# Get percentile corresponding to a value using ecdf function
debug(ecdf)
my_ecdf <- ecdf(x)
重点是
rval <- approxfun(vals, cumsum(tabulate(match(x, vals)))/n,
method = "constant", yleft = 0, yright = 1, f = 0, ties = "ordered")
之后你可以检查
data.frame(x = vals, y = round(cumsum(tabulate(match(x, vals)))/n, 3), stringsAsFactors = FALSE)
并且按照 n=11
进行划分,结果并不令人惊讶。如前所述,对于理论,请查看其他答案。
顺便说一句,你也可以绘制函数
plot(my_ecdf)
关于您的评论。我认为这不是可靠性的问题,而是如何定义 "inverse distribution function, if it does not exist":
的问题
广义逆的一个很好的参考:Paul Embrechts, Marius Hofert: "A note on generalized inverses", Math Meth Oper Res (2013) 77:423–432 DOI
我看过其他帖子(例如here)关于获取"reverse"分位数——即获取一系列值中某个值对应的百分位数。
但是,对于相同的数据系列 ,答案并没有给我与分位数 相同的值。
我也研究过分位数提供了9种不同的算法来计算百分位数。
所以我的问题是:是否有可靠的方法来获取分位数函数的反函数? ecdf 不采用 "type" 参数,因此似乎无法确保他们使用相同的方法。
可重现的例子:
# Simple data
x = 0:10
pcntile = 0.5
# Get value corresponding to a percentile using quantile
(pcntile_value <- quantile(x, pcntile))
# 50%
# 5 # returns 5 as expected for 50% percentile
# Get percentile corresponding to a value using ecdf function
(pcntile_rev <- ecdf(x)(5))
# [1] 0.5454545 #returns 54.54% as the percentile for the value 5
# Not the same answer as quantile produces
ecdf
给出文档中公式的结果。
x <- 0:10
Fn <- ecdf(x)
现在,对象 Fn
是一个插值阶跃函数。
str(Fn)
#function (v)
# - attr(*, "class")= chr [1:3] "ecdf" "stepfun" "function"
# - attr(*, "call")= language ecdf(x)
并保留原始 x
值和相应的 y
值。
environment(Fn)$x
# [1] 0 1 2 3 4 5 6 7 8 9 10
environment(Fn)$y
# [1] 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 0.54545455
# [7] 0.63636364 0.72727273 0.81818182 0.90909091 1.00000000
后者与文档中所说的用于计算它们的公式的结果完全相同。来自 help('ecdf')
:
For observations x= (x1,x2, ... xn), Fn is the fraction of
observations less or equal to t, i.e.,Fn(t) = #{xi <= t}/n = 1/n sum(i=1,n) Indicator(xi <= t).
我将使用 seq_along
而不是 1:length(x)
。
seq_along(x)/length(x)
# [1] 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 0.54545455
# [7] 0.63636364 0.72727273 0.81818182 0.90909091 1.00000000
Fn(x)
# [1] 0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 0.54545455
# [7] 0.63636364 0.72727273 0.81818182 0.90909091 1.00000000
link 中的答案非常好,但也许有帮助,看看 ecdf
只需 运行 以下代码:
# Simple data
x = 0:10
p0 = 0.5
# Get value corresponding to a percentile using quantile
sapply(c(1:7), function(i) quantile(x, p0, type = i))
# 50% 50% 50% 50% 50% 50% 50%
# 5.0 5.0 5.0 4.5 5.0 5.0 5.0
因此,这不是类型的问题。您可以使用 debug:
进入函数# Get percentile corresponding to a value using ecdf function
debug(ecdf)
my_ecdf <- ecdf(x)
重点是
rval <- approxfun(vals, cumsum(tabulate(match(x, vals)))/n,
method = "constant", yleft = 0, yright = 1, f = 0, ties = "ordered")
之后你可以检查
data.frame(x = vals, y = round(cumsum(tabulate(match(x, vals)))/n, 3), stringsAsFactors = FALSE)
并且按照 n=11
进行划分,结果并不令人惊讶。如前所述,对于理论,请查看其他答案。
顺便说一句,你也可以绘制函数
plot(my_ecdf)
关于您的评论。我认为这不是可靠性的问题,而是如何定义 "inverse distribution function, if it does not exist":
的问题广义逆的一个很好的参考:Paul Embrechts, Marius Hofert: "A note on generalized inverses", Math Meth Oper Res (2013) 77:423–432 DOI