Почему nls.lm дает противоречивые результаты для эквивалентных уравнений? - PullRequest
Новая
       16

Почему nls.lm дает противоречивые результаты для эквивалентных уравнений?

0 голосов
/ 04 марта 2019

Ученик заметил, что при одинаковых данных и эквивалентных уравнениях nls.lm последовательно давал разные результаты.

В ближайшем вопросе, который я могу найти, была ошибка в уравнениях, поэтому ниже я приведу код, показывающий, что уравнения эквивалентны.(Q: [Differences in ^ and exp() notation in nls models][1])

Два уравнения, которые мы ожидаем, будут иметь идентичные результаты:

  1. D0 * pars $ A * exp (-pars $ mu1 * tt) +D0 * pars $ B * exp (-pars $ mu2 * tt)
  2. D0 * (pars $ A * exp (-pars $ mu1 * tt) + pars $ B * exp (-pars $ mu2 * tt)))

Я подозревал, что это может быть результатом нахождения различных локальных минимумов, поэтому я повторил оптимизацию для уравнений 1 и 2 (Fit1 и Fit2) 10 раз.Если бы это была разница в локальных минимумах, я бы ожидал несовместимых оптимальных параметров.Однако, когда я запустил этот код (ниже), я получил те же 2 набора оптимальных параметров для каждой итерации.

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

library(reshape)
library(ggplot2)
library(minpack.lm)

mydata <- data.frame("Time.hrs" = c(0,4,8,12,24,48),
                        "values" = c(5.27, 4.60, 4.15, 3.85, 3.40, 3.15))

optimizer_function <- function(Fit_function) {
  # define inputs to nls.lm
  costFun <- function(pars, observed, tt, D0) observed - Fit_function(pars, tt, D0)
  pars0 <- list(A = 1, B = 0.1, mu1 = 0.1, mu2 = 0.1 )
  lb <- c(0.01, 0.00001, 0.01, 0.0001)
  ub <- c(3, 1, 1, 0.15)  
  #optimize
  nls <- nls.lm(par = pars0, lower = lb, upper = ub, fn=costFun,
                    D0 =  mydata$values[1], observed = mydata$values, tt = mydata$Time.hrs,
                 control = nls.lm.control(ftol = 1e-20, ptol = 1e-15))

  return(nls)}

Fit1 <- function(pars, tt, D0) D0 * pars$A * exp(-pars$mu1 * tt) + D0 * pars$B * exp(-pars$mu2 * tt)
Fit2 <- function(pars, tt, D0) D0 * (pars$A * exp(-pars$mu1 * tt) + pars$B * exp(-pars$mu2 * tt))
nls1 <- optimizer_function(Fit1)
nls2 <- optimizer_function(Fit2)

nls1
nls2


results1 <- NULL
results2 <- NULL

## Optimization
for (n in c(1:100)) {
  ## Solve
  nls1 <- optimizer_function(Fit1)
  nls2 <- optimizer_function(Fit2)
  results1 <- as.data.frame(rbind(results1, nls1$par))
  results2 <- as.data.frame(rbind(results2, nls2$par))
}
results1
results2

Чтобы убедиться, что уравнения подгонки фактически эквивалентны, я смоделировал каждый набор оптимальных параметров для каждого уравнения (Fit1 и Fit2).debug_plot показывает, что оценка Fit1 и Fit2 для любого набора оптимальных параметров дает эквивалентные результаты.

resultslist <- list("ParamsfromFit1" = nls1$par, "ParamsfromFit2" = nls2$par)
fitslist <- list("Fit1" = Fit1, "Fit2" = Fit2)

simulation_function <- function(parameters, raw_data, this_fit) {
  t_data <- seq(from = 0, to = 50, by = 0.1)
  Dinitial <- raw_data$value[1]
  y_data <- as.numeric(this_fit(pars = parameters, tt = t_data, D0 = Dinitial))
  this_ydata <- cbind.data.frame(t_data, y_data)
  return(this_ydata)}

sim_long <- data.frame()
for (res in c(1:length(resultslist))) {
  this_parameters <- resultslist[res][[1]]
  for (fit in c(1:length(fitslist))) {
    this_fit <- fitslist[fit][[1]]
    this_data <- simulation_function(this_parameters, mydata, this_fit)
    sim_long <- rbind.data.frame(sim_long,
                                 cbind.data.frame(rep(names(resultslist)[res],nrow(this_data)),
                                                  rep(names(fitslist)[fit],nrow(this_data)),
                                                  this_data))
  }
}

names(sim_long) <- c("Parameters", "Fit", "t_data", "y_data")

debug_plot <- ggplot() +
  geom_line(data = sim_long, aes(x = t_data, y = y_data, color = Parameters, linetype = Fit)) + 
  facet_wrap(~Fit) +
  theme_bw() + #removes default gray plot background
  scale_x_continuous(expand = c(0,0)) + #removes default space between plot and axis
  scale_y_continuous(expand = c(0,0)) 

debug_plot
...