Использование optimize() для поиска кратчайшего интервала, занимающего 95% площади под кривой в R.

Основываясь на вашем пересмотренном вопросе, я нашел оптимизацию, которая минимизирует КРАТЧАЙШЕЕ расстояние (в единицах 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")

Если мы представим это как попытку вычислить интервал с наименьшей площадью, мы можем начать вычислять площади каждой из областей, которые мы наносим на график. Затем мы можем найти самую большую область (которая, предположительно, будет недалеко от центра) и начать идти, пока не найдем область, которую ищем.

Поскольку вы уже вычислили значения 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")

Это, по общему признанию, эквивалентно построению функции интегрирования с константой нормализации в области определения функции «плотности». Тот факт, что все интервалы имеют одинаковую размерность, позволяет исключить разность вектора «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% от общей площади над нулевой линией: