如何使用nls（）拟合指数衰减模型中的多个常数？

Question

我正在处理这种关系：

y = h * R + x * v * h

其中：

x =（N-M）* exp（-Q * u）+ M

给出主要方程式：

y = h * R + v * h *（N-M）* exp（-Q * u）+ v * h * M

所有大写字母都是常量，所有小写字母都是变量。

我具有所有变量的真实数据，但是我要么不知道常数（R和Q）的值，要么想检查数据适合常数（N和M）值的能力。我想使用nls（）来使用变量数据拟合方程，以估计这些常数参数。

我如何使用nls（）函数编写代码以描绘主方程，以允许从模拟测量数据中估计参数R，N，Q和M？ （模拟测量数据=带_j后缀的小写字母，请参见下文。）

创建模拟数据：

library(dplyr)
library(ggplot2)

### Generate mock data

# Equations:
# y = h*R + x*v*h
# x = (N-M)*exp(-Q*u) + M
# y = h*R + ((N-M)*exp(-Q*u) + M)*v*h
# y = h*R + v*h*(N-M)*exp(-Q*u) + v*h*M

### Variables have varying periodicity,
# and so can be approximated via different functions,
# with unique noise added to each to simulate variability:

# Variability for each variable
n <- 1000 # number of data points
t <- seq(0,4*pi,length.out = 1000)
a <- 3
b <- 2

y.norm <- rnorm(n)
u.norm <- rnorm(n)
u.unif <- runif(n)
v.norm <- rnorm(n)
v.unif <- runif(n)
amp <- 1

# Create reasonable values of mock variable data for all variables except h;
# I will calculate from known fixed values for R, N, Q, and M.

y <- 1.5*a*sin(b*t)+y.norm*amp-10 # Gaussian/normal error
u <- ((1*a*sin(11*b*t)+u.norm*amp)+(0.5*a*sin(13*b*t)+u.unif*amp)+7)/2
v <- 1/((2*a*sin(11*b*t)+v.norm*amp)+(1*a*sin(13*b*t)+v.unif*amp)+20)*800-25

# Put vectors in dataframe
dat <- data.frame("t" = t, "y" = y, "u" = u, "v" = v)

### Create reasonable values for constants:

R=0.5
N=1.12
Q=0.8
M=1

### Define final variable based on these constants and the previous
# mock variable data:

dat$h = y/(R + v*(N-M)*exp(-Q*dat$u))

### Gather data to plot relationships:

dat_gathered <- dat %>%
  gather(-t, value = "value", key = "key")

### Plot data to check all mock variables:

ggplot(dat_gathered, aes(x = t, y = value, color = key)) + geom_line()

# Add small error (to simulate measurement error):

dat <- dat %>%
  mutate(h_j = h + rnorm(h, sd=0.05)/(1/h)) %>%
  mutate(u_j = u + rnorm(u, sd=0.05)/(1/u)) %>%
  mutate(v_j = v + rnorm(v, sd=0.05)/(1/v)) %>%
  mutate(y_j = y + rnorm(y, sd=0.05)/(1/y))

Answer 1

[nls似乎可以正常工作，但是解决方案（就参数而言）似乎是不唯一的...

## parameter values chosen haphazardly
n1 <- nls(y ~ h_j*(R + v_j*((N-M)*exp(-Q*u_j)+M)),
    start=list(R=1,N=2,M=1,Q=1),
    data=dat)

## starting from known true values
true_vals <- c(R=0.5,N=1.12,Q=0.8,M=1)
n2 <- update(n1, start=as.list(true_vals))

round(cbind(coef(n1),coef(n2),true_vals),3)
              true_vals
R 0.495 0.495      0.50
N 0.120 0.120      1.12
M 0.001 0.818      0.80
Q 0.818 0.001      1.00

在两个拟合上使用AIC()表示它们具有基本相同的拟合优度（并且预测几乎相同），这表明模型中存在一些对称性，可以互换M和Q。我还没有足够认真地考虑/研究方程式，以至于不知道为什么会这样。