используя функцию R NLS с длинным списком переменных - PullRequest
0 голосов
/ 12 июня 2018

Я пытаюсь оценить нелинейную модель с функцией R nls.

Я проиллюстрирую мою проблему на примере из средства «help» для nls.

Treated <- Puromycin[Puromycin$state == "treated", ]
weighted.MM <- function(resp, conc, Vm, K)
{
    pred <- (Vm * conc)/(K + conc)
    (resp - pred) / sqrt(pred)
}

Pur.wt <- nls( ~ weighted.MM(rate, conc, Vm, K), data = Treated,
          start = list(Vm = 200, K = 0.1))

В этом примере weighted.MM - это функция с очень ограниченным числом аргументов,и у меня нет проблем с реализацией аналогичных подходов для типа модели, с которой я работаю.

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

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

MyArguments <- list(Vm, K)  

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

Error: object 'Vm' not found

В качестве альтернативы,

    MyArguments <- list("Vm", "K") 

weighted.MM1 <- function(resp, conc1, conc.1, thearguments)
{
  conc <- c(conc1, conc.1)
  pred <- (thearguments[[1]] * conc)/(thearguments[[2]] + conc)
  (resp - pred) / sqrt(pred)
}

Pur.wt1 <- nls( ~ weighted.MM1(rate, conc1, conc.1, MyArguments),
                data = lisTreat, start = list(Vm = 200, K = 0.1))

дает:

Error in thearguments[[1]] * conc : 
  non-numeric argument to binary operator
Called from: weighted.MM1(rate, conc1, conc.1, MyArguments)

Есть ли работавокруг для этого?

Ответы [ 2 ]

0 голосов
/ 13 июня 2018

Я пытался реализовать новое решение, основываясь на предложениях, любезно сделанных @ Onyambu.

Однако это создало новые проблемы.

Сначала я попытался реализовать решениес нлс.Вот фактический код, который я использовал:

DiffModel <- nls(COPERTFreq ~ CalculateProbaVehChoiceDiffusion(MyArguments, years = RegistYear),
                 data = DataSetForModel , start = list(MyArguments = c(ASC_Mat, ASC_No_Size)))

, где CalculateProbaVehChoiceDiffusion () - это нелинейная функция, которая была определена в другом месте, RegistYear - это константа, а MyArguments - это список коэффициентов для оценки, где c(ASC_Mat, ASC_No_Size) в качестве начальных значений.

Это приводит к следующему сообщению об ошибке:

Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

Теперь я где-то читал, что эту проблему можно решить с помощью взамен nlsLM.Это приводит к появлению нового сообщения об ошибке:

Error in `rownames<-`(`*tmp*`, value = "MyArguments") : 
  length of 'dimnames' [1] not equal to array extent

ОК, поэтому я снова запустил модель с nls.lm в режиме отладки.Это показывает, что сообщение об ошибке исходит из следующей строки кода:

names(out$par) <- rownames(out$hessian) <- colnames(out$hessian) <- names(out$diag) <- names(par)

Однако именно при проверке объекта «out» становится ясно, в чем заключается проблема.Во-первых, наш $ hessian - это просто скаляр, а я ожидаю, что его число строк и столбцов будет равно числу параметров.Во-вторых, из $ par $ MyArguments видно, что, за исключением первого элемента, значение MyArguments не меняется от одной итерации к другой.

Это известная ошибка или мне нужно изменить способ передачи MyArguments в вызов функции?

Обратите внимание, что, насколько я вижу, эта проблема также возникает, когда я применяю nlsLm к примеру, предоставленному @Onyambu:

> undebug(nls.lm)
> Treated <- Puromycin[Puromycin$state == "treated", ]
> 
> lisTreat <- with(Treated,
+                  list(conc1 = conc[1], conc.1 = conc[-1], rate = rate))
> 
> weighted.MM1 <- function(resp, conc1, conc.1, thearguments)
+ {
+   conc <- c(conc1, conc.1)
+   pred <- (thearguments[1] * conc)/(thearguments[2] + conc)
+   (resp - pred) / sqrt(pred)
+ }
> nls( ~ weighted.MM1(rate, conc1, conc.1, MyArguments),
+      data = lisTreat, start = list(MyArguments =c( 200, 0.1)))
Nonlinear regression model
  model: 0 ~ weighted.MM1(rate, conc1, conc.1, MyArguments)
   data: lisTreat
MyArguments1 MyArguments2 
   206.83468      0.05461 
 residual sum-of-squares: 14.6

Number of iterations to convergence: 5 
Achieved convergence tolerance: 3.858e-06
> 
> nlsLM( ~ weighted.MM1(rate, conc1, conc.1, MyArguments),
+        data = lisTreat, start = list(MyArguments =c( 200, 0.1)))
Error in `rownames<-`(`*tmp*`, value = "MyArguments") : 
  length of 'dimnames' [1] not equal to array extent

Таким образом, хотя nls работает с этим примером,nlsLM не делает, и выдает то же сообщение об ошибке, что и с моим кодом.Затем я снова запускаю nlsLM с nls.lm в режиме отладки.После следующей строки

out <- .Call("nls_lm", par, lower, upper, fn1, jac1, ctrl, 
    new.env(), PACKAGE = "minpack.lm")

я проверяю объект out и вижу:

$par
$par$MyArguments
[1] 244.5117   0.1000


$hessian
[1] 0.02859739

$fvec
 [1] -5.5215487 -0.9787468 -0.5543382 -1.5986605  0.4486608 -0.9651245  0.7020058  1.2419040  1.1430780  0.4488084  1.1445818  1.6121474

$info
[1] 1

$message
[1] "Relative error in the sum of squares is at most `ftol'."

$diag
$diag[[1]]
[1] 0.2077949


$niter
[1] 4

$rsstrace
[1] 112.59784  43.41211  42.89350  42.89349  42.89349

Таким образом, значение второго аргумента не изменилось после 4 итераций.Это может конечно быть правильным решением.Но я нахожу поразительным совпадением то, что то же самое происходит с моей моделью.

Окончательное редактирование : Я наконец-то решил вернуться к полному перечислению всех аргументов.Как я писал в постановке задачи, он не очень элегантен, но, по крайней мере, работает с nls.lm (хотя все еще не с nls).

0 голосов
/ 12 июня 2018

Поскольку вы ищете значения своих аргументов, вам не нужно определять их: посмотрите на код ниже:

weighted.MM1 <- function(resp, conc1, conc.1, thearguments)
{
  conc <- c(conc1, conc.1)
  pred <- (thearguments[1] * conc)/(thearguments[2] + conc)
  (resp - pred) / sqrt(pred)
}
nls( ~ weighted.MM1(rate, conc1, conc.1, MyArguments),
                data = lisTreat, start = list(MyArguments =c( 200, 0.1)))
Nonlinear regression model
  model: 0 ~ weighted.MM1(rate, conc1, conc.1, MyArguments)
   data: lisTreat
MyArguments1 MyArguments2 
   206.83468      0.05461 
 residual sum-of-squares: 14.6

Number of iterations to convergence: 5 
Achieved convergence tolerance: 3.858e-06
...