2016-06-10 4 views
2

Я пытаюсь использовать optimize() найти минимальное значение п для следующей функции (Clopper-Пирсон нижней границы):Оптимизации для глобального минимума

f <- function (n, p=0.5) 
(1 + (n - p*n + 1)/
    (p*n*qf(p= .025, df1= 2*p, df2= 2*(n - p + 1))))^-1 

А вот как я попытался оптимизировать его:

n_clop <- optimize(f.1, c(300,400), maximum = FALSE, p=0.5) 
n_clop 

Я сделал это на интервале [300,400], потому что я подозреваю, что значение находится между ним, но в конечном итоге я хотел бы сделать оптимизацию между 0 и бесконечностью. Кажется, что эта команда производит локальный минимум, потому что независимо от того, какой интервал он дает нижнюю границу этого интервала как минимум, - это не то, что я подозреваю от clopper-pearson. Итак, мои два вопроса: как правильно найти глобальный минимум в R и как это сделать на любом интервале?

+0

Кроме того - я хотел бы, чтобы функция равнялась 0,5 (полуширина доверительного интервала для пропорции) и вычисляла n для этого. Не уверен, что моя настройка верна. –

+0

Кстати, политика SO заключается в том, что вам не нужно включать имя языка программирования (R) в заголовок вопроса - тэг [r] должен быть достаточным –

+0

Нет общего способа сделать это для каждая функция. Что вы готовы принять, чтобы сделать это разрешимым? Глобальная оптимизация - непростая задача для решения любого языка. – MrFlick

ответ

2

Я очень кратко просмотрел Wikipedia page you linked и не видят каких-либо очевидных опечаток в формуле (хотя я чувствую, что это должно быть 0,975 = 1-альфа/2, а не 0,025 = альфа/2?) , Однако оценка функции, которую вы закодировали в очень широком масштабе, предполагает, что нет никаких локальных минимумов, которые вас путают. Мое сильное предположение заключалось бы в том, что либо ваша логика неверна (т. Е. N-> 0 - действительно правильный ответ), либо что вы не кодировали то, что, по вашему мнению, кодируете, из-за опечатки (возможно, в статье в Википедии, хотя это кажется маловероятным) или мысли.

f <- function (n, p=0.5) 
(1 + (n - p*n + 1)/
    (p*n*qf(p= .025, df1= 2*p, df2= 2*(n - p + 1))))^-1 

Убедитесь, что вы получаете правильный ответ в течение интервала, который вы выбрали:

curve(f(x),c(300,400)) 

оценивающих более широкий диапазон (п = 0,00001 до 1000000.):

curve(f(10^x),c(-5,7)) 

enter image description here

Как предлагает @MrFlick, глобальная оптимизация сложна. Вы можете начать с optim(...method="SANN"), но лучший ответ определенно зависит от конкретного случая.

Смежные вопросы