使用 optimize() 找到 R 中曲线下占 95% 面积的最短间隔
Using optimize() to find the shortest interval that takes 95% area under a curve in R
背景：
我有一条曲线，其 Y 值由我下面的小 R 函数生成（整齐地注释).如果你 运行 我的整个 R 代码，你会看到我的曲线（但请记住，它是一个函数，所以如果我更改参数值，我可以获得不同的曲线）：
问题：
显然，可以 determine/assume 多个间隔 将 cover/take 这条曲线下总面积的 95%。但是使用 optimize()，我怎样才能找到这么多可能的 95% 间隔中的 SHORTEST（以 x 值单位）？那么这个最短的95%区间两端对应的x值是多少呢？
注意：像我这样的单峰曲线的最短间隔的想法是有道理的。实际上，最短的将是趋向于高度（y 值）较大的中间的那个，因此 x 值对于预期间隔不需要太大 cover/take 曲线下总面积的95%。
这是我的 R 代码（请 运行 整个代码）：
ppp <- function(f, N, df1, df2, petasq, alpha, beta) {
pp <- function(petasq) dbeta(petasq, alpha, beta)
ll <- function(petasq) df(f, df1, df2, (petasq * N) / (1 - petasq) )
marg <- integrate(function(x) pp(x)*ll(x), 0, 1)[[1]]
po <- function(x) pp(x)*ll(x) / marg
return(po(petasq) )
}
## @@@ END OF MY R FUNCTION.
# Now I use my function above to get the y-values for my plot:
petasq <- seq(0, 1, by = .0001) ## These are X-values for my plot
f <- 30 # a function needed argument
df1 <- 3 # a function needed argument
df2 <- 108 # a function needed argument
N <- 120 # a function needed argument
alpha = 5 # a function needed argument
beta = 4 # a function needed argument
## Now use the ppp() function to get the Y-values for the X-value range above:
y.values <- ppp(f, N, df1, df2, petasq, alpha, beta)
## Finally plot petasq (as X-values) against the Y.values:
plot(petasq, y.values, ty="l", lwd = 3 )
如果我们将此视为尝试计算面积最小的区间，我们可以开始计算我们正在绘制的每个区域的面积。然后我们可以找到最大的区域（大概在中心附近）然后开始往外走，直到找到我们要找的区域。
由于您已经计算了绘图的 x 和 y 值，我将重复使用这些值以节省一些计算。这是该算法的实现
pseduoarea <- function(x, y, target=.95) {
dx <- diff(x)
areas <- dx * .5 * (head(y,-1) + tail(y, -1))
peak <- which.max(areas)
range <- c(peak, peak)
found <- areas[peak]
while(found < target) {
if(areas[range[1]-1] > areas[range[2]+1]) {
range[1] <- range[1]-1
found <- found + areas[range[1]-1]
} else {
range[2] <- range[2]+1
found <- found + areas[range[2]+1]
}
}
val<-x[range]
attr(val, "indexes")<-range
attr(val, "area")<-found
return(val)
}
我们称它为
pseduoarea(petasq, y.values)
# [1] 0.3194 0.5413
这确实假设 petasq 中的所有值都是等间距的
我认为您不需要使用优化（除非这是未经认可的家庭作业的一部分）。相反，只需将累积和归一化并找出满足您的标准的点：
> which(cusm.y >= 0.025)[1]
[1] 3163
> which(cusm.y >= 0.975)[1]
[1] 5375
您可以检查这些是用于从 petasq 向量中提取值的合理索引：
abline( v= c( petasq[ c( which(cusm.y >= 0.025)[1], which(cusm.y >= 0.975)[1])]),
col="red")
这无可否认等同于在 "density" 函数的域中构造一个具有归一化常数的积分函数。间隔都是等维的事实允许从高度乘以基础计算中省略 "x"-向量的差分。
我想还有另一种可能的解释。这将要求我们发现 petasq 的 ascending-sorted 版本需要多少个值才能总和达到总和的 95%。这给出了不同的策略，绘图显示了水平线与曲线相交的位置：
which( cumsum( sort( y.values, decreasing=TRUE) ) > 0.95* sum(y.values, na.rm=TRUE) )[1]
#[1] 2208
sort( y.values, decreasing=TRUE)[2208]
#[1] 1.059978
png()
plot(petasq, y.values, ty="l", lwd = 3 )
abline( h=sort( y.values, decreasing=TRUE)[2208], col="blue")
dev.off()
要获得 petasq 值，您需要确定第一个 y.values 超过该值，然后确定下一个 y.values 低于该值。这些可以通过以下方式获得：
order(y.values, decreasing=TRUE)[2208]
#[1] 3202
order(y.values, decreasing=TRUE)[2209]
#[1] 5410
然后情节看起来像：
png(); plot(petasq, y.values, ty="l", lwd = 3 )
abline( v= petasq[ c(3202, 5410)], col="blue", lty=3, lwd=2)
dev.off()
两条蓝色虚线之间的面积占零线以上总面积的95%：
根据您修改后的问题，我找到了使 LEFT 和 RIGHT 边界之间的最短距离（以 x 值为单位）最小化的优化：
ppp <- function(petasq, f, N, df1, df2, alpha, beta) {
pp <- function(petasq) dbeta(petasq, alpha, beta)
ll <- function(petasq) df(f, df1, df2, (petasq * N) / (1 - petasq) )
marg <- integrate(function(x) pp(x)*ll(x), 0, 1)[[1]]
po <- function(x) pp(x)*ll(x) / marg
return(po(petasq) )
}
petasq <- seq(0, 1, by = .0001) ## These are X-values for my plot
f <- 30 # a function needed argument
df1 <- 3 # a function needed argument
df2 <- 108 # a function needed argument
N <- 120 # a function needed argument
alpha = 5 # a function needed argument
beta = 4 # a function needed argument
optim_func <- function(x_left) {
int_function <- function(petasq) {
ppp(petasq, f=f, N=N, df1=df1, df2=df2, alpha=alpha, beta=beta)
}
# For every LEFT value, find the corresponding RIGHT value that gives 95% area.
find_95_right <- function(x_right) {
(0.95 - integrate(int_function, lower=x_left, upper=x_right, subdivisions = 10000)$value)^2
}
x_right_obj <- optimize(f=find_95_right, interval=c(0.5,1))
if(x_right_obj$objective > .Machine$double.eps^0.25) return(100)
#Return the DISTANCE BETWEEN LEFT AND RIGHT
return(x_right_obj$minimum - x_left)
}
#MINIMIZE THE DISTANCE BETWEEN LEFT AND RIGHT
x_left <- optimize(f=optim_func, interval=c(0.30,0.40))$minimum
find_95_right <- function(x_right) {
(0.95 - integrate(int_function, lower=x_left, upper=x_right, subdivisions = 10000)$value)^2
}
int_function <- function(petasq) {
ppp(petasq, f=f, N=N, df1=df1, df2=df2, alpha=alpha, beta=beta)
}
x_right <- optimize(f=find_95_right, interval=c(0.5,1))$minimum
查看代码中的注释。希望这最终能解决您的问题 :) 结果：
> x_right
[1] 0.5409488
> x_left
[1] 0.3201584
此外，您可以绘制 LEFT 和 RIGHT 之间的距离作为左边界的函数：
left_x_values <- seq(0.30, 0.335, 0.0001)
DISTANCE <- sapply(left_x_values, optim_func)
plot(left_x_values, DISTANCE, type="l")
背景：
我有一条曲线，其 Y 值由我下面的小 R 函数生成（整齐地注释).如果你 运行 我的整个 R 代码，你会看到我的曲线（但请记住，它是一个函数，所以如果我更改参数值，我可以获得不同的曲线）：
问题：
显然，可以 determine/assume 多个间隔 将 cover/take 这条曲线下总面积的 95%。但是使用 optimize()，我怎样才能找到这么多可能的 95% 间隔中的 SHORTEST（以 x 值单位）？那么这个最短的95%区间两端对应的x值是多少呢？
注意：像我这样的单峰曲线的最短间隔的想法是有道理的。实际上，最短的将是趋向于高度（y 值）较大的中间的那个，因此 x 值对于预期间隔不需要太大 cover/take 曲线下总面积的95%。
这是我的 R 代码（请 运行 整个代码）：
ppp <- function(f, N, df1, df2, petasq, alpha, beta) {
pp <- function(petasq) dbeta(petasq, alpha, beta)
ll <- function(petasq) df(f, df1, df2, (petasq * N) / (1 - petasq) )
marg <- integrate(function(x) pp(x)*ll(x), 0, 1)[[1]]
po <- function(x) pp(x)*ll(x) / marg
return(po(petasq) )
}
## @@@ END OF MY R FUNCTION.
# Now I use my function above to get the y-values for my plot:
petasq <- seq(0, 1, by = .0001) ## These are X-values for my plot
f <- 30 # a function needed argument
df1 <- 3 # a function needed argument
df2 <- 108 # a function needed argument
N <- 120 # a function needed argument
alpha = 5 # a function needed argument
beta = 4 # a function needed argument
## Now use the ppp() function to get the Y-values for the X-value range above:
y.values <- ppp(f, N, df1, df2, petasq, alpha, beta)
## Finally plot petasq (as X-values) against the Y.values:
plot(petasq, y.values, ty="l", lwd = 3 )
如果我们将此视为尝试计算面积最小的区间，我们可以开始计算我们正在绘制的每个区域的面积。然后我们可以找到最大的区域（大概在中心附近）然后开始往外走，直到找到我们要找的区域。
由于您已经计算了绘图的 x 和 y 值，我将重复使用这些值以节省一些计算。这是该算法的实现
pseduoarea <- function(x, y, target=.95) {
dx <- diff(x)
areas <- dx * .5 * (head(y,-1) + tail(y, -1))
peak <- which.max(areas)
range <- c(peak, peak)
found <- areas[peak]
while(found < target) {
if(areas[range[1]-1] > areas[range[2]+1]) {
range[1] <- range[1]-1
found <- found + areas[range[1]-1]
} else {
range[2] <- range[2]+1
found <- found + areas[range[2]+1]
}
}
val<-x[range]
attr(val, "indexes")<-range
attr(val, "area")<-found
return(val)
}
我们称它为
pseduoarea(petasq, y.values)
# [1] 0.3194 0.5413
这确实假设 petasq 中的所有值都是等间距的
我认为您不需要使用优化（除非这是未经认可的家庭作业的一部分）。相反，只需将累积和归一化并找出满足您的标准的点：
> which(cusm.y >= 0.025)[1]
[1] 3163
> which(cusm.y >= 0.975)[1]
[1] 5375
您可以检查这些是用于从 petasq 向量中提取值的合理索引：
abline( v= c( petasq[ c( which(cusm.y >= 0.025)[1], which(cusm.y >= 0.975)[1])]),
col="red")
这无可否认等同于在 "density" 函数的域中构造一个具有归一化常数的积分函数。间隔都是等维的事实允许从高度乘以基础计算中省略 "x"-向量的差分。
我想还有另一种可能的解释。这将要求我们发现 petasq 的 ascending-sorted 版本需要多少个值才能总和达到总和的 95%。这给出了不同的策略，绘图显示了水平线与曲线相交的位置：
which( cumsum( sort( y.values, decreasing=TRUE) ) > 0.95* sum(y.values, na.rm=TRUE) )[1]
#[1] 2208
sort( y.values, decreasing=TRUE)[2208]
#[1] 1.059978
png()
plot(petasq, y.values, ty="l", lwd = 3 )
abline( h=sort( y.values, decreasing=TRUE)[2208], col="blue")
dev.off()
要获得 petasq 值，您需要确定第一个 y.values 超过该值，然后确定下一个 y.values 低于该值。这些可以通过以下方式获得：
order(y.values, decreasing=TRUE)[2208]
#[1] 3202
order(y.values, decreasing=TRUE)[2209]
#[1] 5410
然后情节看起来像：
png(); plot(petasq, y.values, ty="l", lwd = 3 )
abline( v= petasq[ c(3202, 5410)], col="blue", lty=3, lwd=2)
dev.off()
两条蓝色虚线之间的面积占零线以上总面积的95%：
根据您修改后的问题，我找到了使 LEFT 和 RIGHT 边界之间的最短距离（以 x 值为单位）最小化的优化：
ppp <- function(petasq, f, N, df1, df2, alpha, beta) {
pp <- function(petasq) dbeta(petasq, alpha, beta)
ll <- function(petasq) df(f, df1, df2, (petasq * N) / (1 - petasq) )
marg <- integrate(function(x) pp(x)*ll(x), 0, 1)[[1]]
po <- function(x) pp(x)*ll(x) / marg
return(po(petasq) )
}
petasq <- seq(0, 1, by = .0001) ## These are X-values for my plot
f <- 30 # a function needed argument
df1 <- 3 # a function needed argument
df2 <- 108 # a function needed argument
N <- 120 # a function needed argument
alpha = 5 # a function needed argument
beta = 4 # a function needed argument
optim_func <- function(x_left) {
int_function <- function(petasq) {
ppp(petasq, f=f, N=N, df1=df1, df2=df2, alpha=alpha, beta=beta)
}
# For every LEFT value, find the corresponding RIGHT value that gives 95% area.
find_95_right <- function(x_right) {
(0.95 - integrate(int_function, lower=x_left, upper=x_right, subdivisions = 10000)$value)^2
}
x_right_obj <- optimize(f=find_95_right, interval=c(0.5,1))
if(x_right_obj$objective > .Machine$double.eps^0.25) return(100)
#Return the DISTANCE BETWEEN LEFT AND RIGHT
return(x_right_obj$minimum - x_left)
}
#MINIMIZE THE DISTANCE BETWEEN LEFT AND RIGHT
x_left <- optimize(f=optim_func, interval=c(0.30,0.40))$minimum
find_95_right <- function(x_right) {
(0.95 - integrate(int_function, lower=x_left, upper=x_right, subdivisions = 10000)$value)^2
}
int_function <- function(petasq) {
ppp(petasq, f=f, N=N, df1=df1, df2=df2, alpha=alpha, beta=beta)
}
x_right <- optimize(f=find_95_right, interval=c(0.5,1))$minimum
查看代码中的注释。希望这最终能解决您的问题 :) 结果：
> x_right
[1] 0.5409488
> x_left
[1] 0.3201584
此外，您可以绘制 LEFT 和 RIGHT 之间的距离作为左边界的函数：
left_x_values <- seq(0.30, 0.335, 0.0001)
DISTANCE <- sapply(left_x_values, optim_func)
plot(left_x_values, DISTANCE, type="l")