1

背景：

我有一条曲线，其 Y 值由下面的小 R 函数产生（注释整齐）。如果你运行我的整个 R 代码，你会看到我的曲线（但请记住，它是一个函数，所以如果我更改参数值，我可以得到不同的曲线）：

在此处输入图像描述

问题：

显然，可以确定/假设许多区间将覆盖/占据该曲线下总面积的 95%。但是使用, optimize()，我怎样才能找到这些许多可能的 95% 间隔中的最短（以 x 值单位） ？那么这个最短 95% 区间的两端对应的 x 值是多少？

注意：像我这样的单峰曲线最短间隔的想法是有道理的。实际上，最短的将是趋向于高度（y 值）较大的中间的那个，因此对于预期的间隔覆盖/取 95，x 值不需要太大曲线下总面积的百分比。

这是我的 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 )
4

3 回答 3

2

根据您修改后的问题，我发现了最小化左右边界之间最短距离（以 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")

阴谋

于 2017-04-08T02:49:59.533 回答
2

如果我们认为这是试图计算面积最小的区间，我们可以开始计算我们正在绘制的每个区域的面积。然后我们可以找到最大的区域（大概在中心附近）并开始向外走，直到找到我们正在寻找的区域。

由于您已经计算了绘图的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都是等距的

于 2017-04-07T20:57:54.357 回答
0

我认为您不需要使用优化（除非这是未经批准的家庭作业的一部分）。相反，只需将累积总和归一化，并找出满足您的条件的点：

> 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")

诚然，这等效于在“密度”函数的域中构造一个具有归一化常数的积分函数。间隔都是相同维度的事实允许从高度乘以基数计算中省略“x”向量的差异。

在此处输入图像描述

我想可能还有另一种解释。这将要求我们发现需要多少个升序版本的值petasq才能求和到总和的 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%：

在此处输入图像描述

于 2017-04-07T21:23:27.857 回答