Оценка временного интервала нелинейной регрессии в R - диффузионных моделях - PullRequest
Новая
       60

Оценка временного интервала нелинейной регрессии в R - диффузионных моделях

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

Как вы оцениваете временные интервалы нелинейной регрессии в R?

Мы стремимся запустить модель диффузии Нортона Баса, в которой у нас есть три неизвестных параметра: m, p и q (потенциальный размер рынка,параметр инновации и параметр имитации).Мы бы хотели запустить расширенную регрессионную модель диффузии Баса (см. Рисунок 1 и 2).

Функция задается продажами = m1 * F1 (t) -m1 * F1 (t) * F2 (t-t2).

F (t) = ((1-е^ - (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)

Sales functions sought estimated

Где F (t) определяется как:

Definition of F

1 Ответ

0 голосов
/ 17 мая 2019

Я рассмотрел параметры от 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

enter image description here

...