R: Доверительные интервалы нелинейной подгонки с неаналитической моделью.

Мне нужно сопоставить данные xy с моделью, которая не является аналитической. У меня есть функция f(x), которая численно вычисляет модель для каждого x, но нет аналитического уравнения. Для подгонки я использую optim в R. Я минимизирую RMS между моделью и данными. Он работает хорошо и возвращает разумные параметры.

Я хотел бы найти доверительные интервалы (или, по крайней мере, стандартные ошибки) для наиболее подходящих параметров. Я нашел в Интернете, что это можно сделать из матрицы Гессе, но только если максимизировать функцию логарифмического правдоподобия. Я не знаю, как это сделать, все, что у меня есть, это x, y и f(x), из которых я нахожу RMS. Увы, у меня нет хорошего способа оценить ошибки на y.

Как я могу найти доверительные интервалы для моих подходящих параметров?

Изменить: возможно, пример в R может помочь объяснить, о чем я прошу. В этом примере используется простая аналитическая функция для подбора данных, в моем реальном случае функция не аналитическая, поэтому я не могу использовать, например, nls.

set.seed(666)

# generate data
x <- seq(100) / 100
y <- 0.5 * x + rnorm(100, sd = 0.03) + 0.2

# function to fit
f <- function(x, a, b) {
  a * x + b
}

# error function to minimise: RMS
errfun <- function(par, x, y) {
  a <- par[1]
  b <- par[2]
  err <- sqrt(sum((f(x, a, b) - y)^2))
}

# use optim to fit the model to the data
par <- c(1, 0)
res <- optim(par, errfun, gr=NULL, x, y)

# best-fitting parameters
best_a <- res$par[1]
best_b <- res$par[2]

Наилучшие параметры подгонки: a = 0,50 и b = 0,20. Мне нужно найти 95% доверительные интервалы для них.


r non-linear-regression
person Marek Gierliński    schedule 08.01.2018    source источник
comment
Если нет (простого) аналитического выражения для вероятности, вам, вероятно, лучше использовать (непараметрический) бутстрап. Вычислите f(x*) для множества различных x*, выбранных с заменой из x.   -  person JDL    schedule 08.01.2018
comment
Как найти доверительные интервалы из бутстрапа? Я не очень четко выразил свои определения, функция, которую я подгоняю к данным, - это f (x; a, b, c), где a, b и c - параметры модели.   -  person Marek Gierliński    schedule 08.01.2018
comment
Я превратил его в полный ответ. Надеюсь, это более полезно (при условии, что вы используете стандартную запись y как результат, x как данные, f(x) как своего рода оценку и a,b,c как независимые параметры модели (например, количество итераций), которые не зависят от x) .   -  person JDL    schedule 08.01.2018

Ответы (1)


arrow_upward
0
arrow_downward

Это работа для бутстрапа:

(1) создать большое количество синтетических наборов данных x*. Они создаются путем выборки из x с заменой того же количества данных, что и в x. Например, если ваши данные (1,2,3,4,5,6), x* может быть (5,2,4,4,2,3) (обратите внимание, что значения могут появляться несколько раз или не появляться вовсе, потому что мы делаем выборку с заменой)

(2) Для каждого x* рассчитайте f(x*). Если есть другие параметры, не зависящие от данных, не меняйте их. (поэтому f(x,a,b,c) становится f(x*,a,b,c), пока a,b,c не зависит от x. Назовите эти величины f*.

(3) Из этих f* вы можете оценить все, что захотите. Если вам нужно стандартное отклонение f(x), возьмите стандартное отклонение f*. Если вам нужен доверительный интервал 95%, возьмите диапазон от 2,5 до 97,5 процентилей f*. Более формально, если вы хотите оценить g(f(x)), вы оцениваете его как g(f(x*)).

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

Чтобы применить это к примеру, который вы привели в своем коде:

x <- seq(100) / 100
y <- 0.5 * x + rnorm(100, sd = 0.03) + 0.2

# function to fit
f <- function(x, a, b) {
  a * x + b
}

# error function to minimise: RMS
errfun <- function(par, x, y) {
  a <- par[1]
  b <- par[2]
  err <- sqrt(sum((f(x, a, b) - y)^2))
}

# this is the part where we bootstrap
# use optim to fit the model to the data
best_a <- best_b <- numeric(10000)
for(i in 1:10000){
  j <- sample(100,replace=TRUE)
  x.boot <- x[j]; y.boot <- y[j]
par <- c(1, 0)
res <- optim(par, errfun, gr=NULL, x.boot, y.boot)

# best-fitting parameters
best_a[i] <- res$par[1]
best_b[i] <- res$par[2]
}
# now, we look at the *vector* best_a
# for example, if you want the standard deviation of a,
sd(best_a)
# or a 95% confidence interval for b,
quantile(best_b,c(0.025,0.975))
person JDL    schedule 08.01.2018
comment
Спасибо, JDL. Это очень ясное объяснение. Однако мне нужны доверительные интервалы для подходящих параметров a, b и c. То есть, помимо наиболее подходящего a, мне нужны a_lo и a_up, которые, например, составляют 95% CI для a. - person Marek Gierliński; 09.01.2018
comment
Тогда это часть результатов f(x). Поместив их в скобки, они выглядят так, как будто они указаны пользователем. f(x) разрешено иметь несколько элементов. - person JDL; 09.01.2018
comment
Извините, я все еще не понимаю, как начальная загрузка x может помочь мне с моей проблемой. Я расширил вопрос, включив в него простой пример. Возможно, вы могли бы объяснить, как применить загрузку на этом примере. - person Marek Gierliński; 10.01.2018
comment
Ах! Я вижу сейчас! Это выборка с заменой из данных (x, y), подгонка каждой выборки и нахождение параметров a* и b*. Я сделал довольно много начальной загрузки в своей жизни, но не видел такого подхода. Но теперь я вижу, что это имеет смысл. Большое спасибо за ваши усилия! - person Marek Gierliński; 11.01.2018
comment
Рад помочь :) - person JDL; 11.01.2018