Используя optimize(), чтобы найти кратчайший интервал, который занимает 95% площади под кривой в R
Фон:
У меня есть кривая, чьи значения Y производятся моей маленькой функцией R ниже (аккуратно аннотировано). Если вы запустите весь мой код R, вы увидите мою кривую (но помните, что это функция, поэтому, если я изменил значения аргумента, я мог бы получить другую кривую):
Вопрос:
Очевидно, что можно определить / предположить много интервалов, которые бы покрывали / занимали 95% общей площади под этой кривой. Но используя, optimize()Как я могу найти самый короткий (в единицах значения х) из этих многих возможных 95% интервалов? Каковы же будут соответствующие значения x для двух концов этого кратчайшего 95% -ого интервала?
Примечание: идея кратчайшего интервала для унимодальной кривой, подобной моей, имеет смысл. В действительности, самым коротким из них будет тот, который стремится к середине, где высота (значение y) больше, поэтому значение x не должно быть таким большим, чтобы намеченный интервал покрывал / принимал 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 )
3 ответа
Основываясь на вашем пересмотренном вопросе, я нашел оптимизацию, которая минимизирует самое короткое расстояние (в единицах значения 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% от общей площади над нулевой линией:
