2016-08-30 14 views
2

Я хотел бы создать обертку для функции нелинейного наименьшего квадрата Levenberg-Marquardt nls.lm (библиотека minpack.lm), аналогичную nls2 (библиотека nls2), чтобы дать метод грубой силы для оценки соответствия модели наблюдаемым данным.R - Использование вложенных данных для запуска функции с различными наборами параметров

Идея заключается в том, чтобы создать диапазон стартовых комбинаций значений и либо:

  • передать их в функцию, а затем сравнить выходную функцию для наблюдаемых данных для создания 2 значения в R^для каждого из начинать комбинации значений и запускать nls.lm с лучшим из них.

или

  • запустить nls.lm на всех комбинациях и выбрать лучший возвращенное подгонку.

Я хотел бы сделать это без зацикливания и после того, как вдохновения из here пытается использовать вложенные dataframes, с одной колонкой для ввода списка параметров, по одному для значений, возвращаемых моей функцией, один для 2 значений R^и один для лучших подходят модели, что-то вроде:

df 
# start_val fun_out  R^2 
# 1 {a=2,b=2} {22,24,26...} 0.8 
# 2 {a=3,b=5} {35,38,41...} 0.6 

Это код, который я до сих пор:

require(dplyr);require(tidyr) 

foo <- function(x,a,b) a*x^2+b # function I am fitting 
x <- 1:10 # independent variable 
y_obs <- foo(x,1.5,2.5) + rnorm(length(x),0,10) # observed data (dependent variable) 

start_range <- data.frame(a=c(1,2),b=c(2,3)) # range of allowed starting points for fitting 
reps <- 2 # number of starting points to generate 

# Create a data frame of starting points 
df<-as.data.frame(sapply(start_range, function(x) runif(reps,min=x[[1]],max=x[[2]]))) %>% 
    mutate(id=seq_len(reps)) %>% # fudge to make nest behave as I want 
    nest(1:ncol(start_range)) %>% 
    mutate(data=as.list(data)) %>% 
    as.data.frame() 

df 
# id    data 
# 1 1 1.316356, 2.662923 
# 2 2 1.059356, 2.723081 

зависание теперь пытается передать параметры в данных в функции foo(). Я попытался с помощью do.call(), и даже при использовании постоянных параметров следующее сообщение об ошибке появляется:

mutate(df,y=do.call(foo,list(x,1,2))) 
# Error: wrong result size (5), expected 2 or 1 

Есть ли способ, чтобы создать столбцы dataframe, которые содержат списки непосредственно без использования nest()?

Также при попытке создать список для перехода к do.call() с использованием столбцов dataframe, как вы создаете список, в котором первым элементом является вектор x, второй параметр a, а третий - параметр b? Последовательность разделяет список по колонке:

mutate(df,my_list=list(x,data)) 
# id    data        my_list 
# 1 1 1.316356, 2.662923   1, 2, 3, 4, 5, 6, 7, 8, 9, 10 
# 2 2 1.059356, 2.723081 1.316356, 2.662923, 1.059356, 2.723081 
+1

Вам нужно уловить ошибки из 'nls.lm' в вашей функции. Я предлагаю адаптировать исходный код 'nls2' (который, конечно же, не использует dplyr). – Roland

+0

Спасибо @ Роланд, этот подход работал. – lapsel

ответ

1

Бег nls2 использованием algorithm = "random-search" и all = TRUE и указанный maxiter будет оценивать foo на maxiter случайных точек и возврата starting_fits которые являются припадки в этих точках. Он состоит из набора объектов класса "nls", оцененных по каждому из случайно выбранных начальных значений. Он не выполняет оптимизацию с каждого из этих начальных значений, а просто возвращает объект "nls". То есть, nls является не перспективе. Теперь для каждого начального запуска введите nlsLM, давая fits, список nlsLM подходит и из них суммирует их в data (кадр данных с одной строкой за каждый прогон) и покажите как минимум.

Если мы хотим только выбрать наилучшее начальное значение и просто запустите nlsLM, то используйте альтернативный код рядом с концом.

library(nls2) 

fo <- y_obs ~ foo(x, a, b) 
starting_fits <- nls2(fo, algorithm = "random-search", 
start = start_range, control = nls.control(maxiter = reps), all = TRUE) 

fits <- lapply(starting_fits, function(fit) nlsLM(fo, start = coef(fit))) 

data <- data.frame(RSS = sapply(fits, deviance), t(sapply(fits, coef)), 
    start = t(sapply(starting_fits, coef))) 
# data$fits <- fits # optional to store each row's fitted object in that row 
subset(data, RSS == min(RSS)) # minimum(s) 

дает:

 RSS  a  b start.a start.b 
2 706.3956 1.396616 7.226525 1.681819 2.768374 

R квадрат используется для линейной регрессии. Он недействителен для нелинейной регрессии. Вместо этого показана остаточная сумма квадратов (RSS).

Альтернативно, если вы просто хотите выбрать наилучшее начальное значение и запустить nlsLM на этом, то просто опустите аргумент all=TRUE из звонка nls2 и сделайте это. Если вам нужны коэффициенты и RSS для последующего кода, попробуйте coef(fit) и deviance(fit).

starting_fit <- nls2(fo, algorithm = "random-search", 
start = start_range, control = nls.control(maxiter = reps)) 

fit <- nlsLM(fo, start = coef(starting_fit)) 

Примечание 1: Если вы получаете ошибки из nlsLM попробуйте заменить nlsLM(...) с try(nlsLM(...)). Это приведет к появлению сообщений об ошибках (используйте try(..., silent = TRUE), если они вам не нужны), но не прекратит обработку.

Примечание 2: Я предполагаю, что приведенный в вопросе foo является лишь примером, а реальная функция сложнее. Показанный foo является линейным по коэффициентам, поэтому для него можно использовать lm. Не требуется нелинейная оптимизация.

+0

Сделали некоторые обновления и исправления. –

2

Возможно, такой подход?

library(dplyr) 
library(purrr) 

foo2 <- function(x,data) data$a*x^2+data$b 
r2 <- function(e, o) 1 - sum((e - 0)^2)/sum((e - mean(e)^2)) 

df <- as.data.frame(sapply(start_range, function(x) runif(reps,min=x[[1]],max=x[[2]]))) %>% 
    mutate(id=seq_len(reps)) %>% # fudge to make nest behave as I want 
    nest(1:ncol(start_range)) 

df %>% 
    mutate(fun_out = map(data, foo2, x = x), 
     R2 = map(fun_out, o = y_obs, r2)) 

Результат:

# A tibble: 3 x 4 
    id    data fun_out  R2 
    <int>   <list>  <list> <list> 
1  1 <tibble [1 x 2]> <dbl [10]> <dbl [1]> 
2  2 <tibble [1 x 2]> <dbl [10]> <dbl [1]> 
3  3 <tibble [1 x 2]> <dbl [10]> <dbl [1]>