R - 扩散模型中的非线性回归时间间隔估计答案

【问题标题】：Non-linear regression time interval estimation in R - diffusion modelsR - 扩散模型中的非线性回归时间间隔估计
【发布时间】：2019-07-24 20:42:12
【问题描述】：

您如何估计 R 中的非线性回归时间间隔？

我们正在寻求运行广义诺顿巴斯扩散模型，其中我们有三个未知参数：m、p 和 q（潜在市场规模、创新参数和模仿参数）。我们想运行扩展的 Bass 扩散模型回归（见图 1 和图 2）。

函数由sales = m1*F1(t)-m1*F1(t)*F2(t-t2)给出。

F(t) = ((1-e^-(p+g)*t)/((q/p)*e^-((p+g)*t)+1))

我们目前已运行以下代码，但不确定如何在回归中定义 F2(t-t2)？你会如何建议这样做？我们需要估计参数 m、q 和 p

GNB.model.s1 <- nls(s1 ~ 
                      M * (1 - (exp(-(P+Q) * t1)))/(1 + (Q/P) * (exp(-(P+Q) * t1)))
                    - M * (1 - (exp(-(P+Q) * t1)))/(1 + (Q/P) * (exp(-(P+Q) * t1)))
                    * ( (1 - (exp(-(P+Q) * t1)))/(1 + (Q/P) * (exp(-(P+Q) * t1)))
                        - (1 - (exp(-(P+Q) * t2)))/(1 + (Q/P) * (exp(-(P+Q) * t2)))),
                    start = list(M=20000, P=0.03, Q=0.38), trace = T)

其中 F(t) 由下式给出：

【问题讨论】：

标签： r integral nonlinear-optimization non-linear-regression

【解决方案1】：

我考虑了来自https://pubsonline.informs.org/doi/abs/10.1287/mnsc.1120.1529 的参数下面是使用 DE 的解决方案：

fgt1 = function(params,t){
  P=params[1];Q=params[2]
 pf=(1-exp(-(P+Q)*t))/((Q/P)*exp(-(P+Q)*t)+1)
 pf[t<0]=0
 pf
}
fgt2 = function(params,t){  
  P=params[1];Q=params[3]
  pf=(1-exp(-(P+Q)*t))/((Q/P)*exp(-(P+Q)*t)+1)
  pf[t<0]=0
  pf
}

st = 1:25  # time
pt2=11  # 2006-1984   1995-1984
Params0=c(0.009,0.33,0.5,5*10^7,21*10^7)

S1=function(params,t,t2){
  m1=params[4]
  m1*fgt1(params,t)-m1*fgt1(params,t)*fgt2(params,t-t2)
}

S2=function(params,t,t2){
  m1=params[4];m2=params[5]
  (m2 + m1*fgt1(params,t))*fgt2(params,t-t2)
}

#Simulate some data, use set.seed(324)
Sales = S1(Params0,st,pt2) + rnorm(length(st),sd=35)
Sales2 = S2(Params0,st,pt2) + ifelse(st<pt2,0,rnorm(length(st),sd=35))

sd=data.frame(Sales,Sales2)

library(NMOF)

algo1 <- list(printBar = FALSE,
              nP  = 200L,
              nG  = 1000L,
              F   = 0.50,
              CR  = 0.99,
              min = c(.01,.3,.4,10000,10000),
              max = c(.5,.5,.45,6*10^7,22*10^7))  # this appears critical

OF1 <- function(Param, data) {   
  t <- data$x   
  sg <- data$y
  t2 <- data$t2
  s1e <- data$model1(Param,t,t2);
  s2e <- data$model2(Param,t,t2);
  aux <- c(sg[,1],sg[,2]) - c(s1e,s2e); auxs <- sum(aux^2)
  if (is.na(auxs)) auxs <- 1e10 # for bad solutions!
  auxs
}

repair <- function(b,data) { # may be improved
  b[1:5]=abs(b[1:5])
  if(b[1]==0)b[1]=.01  # for Q/P
  b
}

algo1$repair=repair
data1 <- list(x = st, y = sd,  t2=11,model1 = S1,model2 = S2, ww = 1)
system.time(sol1 <- DEopt(OF = OF1, algo = algo1, data = data1))
sol1$xbest
OF1(sol1$xbest,data1)

plot(Sales,ylim=range(c(Sales,Sales2)),type="b",col=2)
points(data1$x[data1$x>=11],Sales2[data1$x>=11],col=3,pch=2)
lines(data1$x,data1$model1(sol1$xbest, data1$x,11),col=6,lwd=2)
lines(data1$x[data1$x>=11],data1$model2(sol1$xbest, data1$x[data1$x>=11],11),col=7,lwd=2)

#> sol1$xbest
#[1] 9.000012e-03 3.299996e-01 4.999998e-01 5.000003e+07 2.100000e+08

【讨论】：