読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R を使ってバスモデルを当てはめてみた

R 微分方程式 時系列

バスモデルのなんたるかについては バスモデル - ORWiki を参照。

バスモデルは以下の微分方程式で記述される。

 \displaystyle dx_{t}/dt=\left(a+{\frac{b}{m}}x_{t}\right)(m-x_t)

閉じた形で解が求まる。

\displaystyle x_{t}=m[1-c_{0}\exp \{-(a+b)t\}]/[\frac{b}{a} c_{0}\exp \{-(a+b)t\}+1]

検算

deSolve パッケージを使って数値的に解いた値と解析解をくらべて、この解が正しいことを一応確かめた。

f:id:abrahamcow:20151017005346p:plain

丸が数値解、曲線が解析解。

#Rのコード
Bass <- function(t,m,a,b,c){
  m*(1-c*exp(-(a+b)*t))/((b/a)*c*exp(-(a+b)*t)+1)
}
library(deSolve)
deBass <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dx = (a+(b/m)*x)*(m-x)
    return(list(c(dx)))
  })
}
ini  <- c(x = 0)
times <- seq(0,100,by=1)
pars  <- c(a=0.05,b=0.01,m=1000)
out1   <- ode(ini,times, deBass, pars)
curve(Bass(x,1000,0.05,0.01,1),0,100)
points(out1)

例1:インターネット利用者数

続いてインターネット利用者数のデータ(総務省|平成25年版 情報通信白書|インターネットの利用状況)にバスモデルとロジスティック曲線をそれぞれ当てはめてみた。

f:id:abrahamcow:20160412204249p:plain

x <- c(6942,7730,7948,8529,8754,8811,9091,9408,9462,9610,9652)
tim <- 1:length(x)
logis <- function(x,pars){
  a <- pars[1]
  b <- pars[2]
  c <- pars[3]
  a/(b+exp(-c*x))
}
Bass <- function(t,pars){
  m <-pars[1]
  a <-pars[2]
  b <-pars[3]
  c <-pars[4]
  m*(1-c*exp(-(a+b)*t))/((b/a)*c*exp(-(a+b)*t)+1)
}
fit_logis <-optim(par=log(c(10000,0.2,0.3)),
                 fn=function(pars)sum((x-logis(tim,exp(pars)))^2),
                 control=list(maxit=2000))
fit_bass <-optim(par=log(c(10000,0.3,0.4,0.3)),
                  fn=function(pars)sum((x-Bass(tim,exp(pars)))^2),
                 control=list(maxit=2000))

plot(x,main="logistic vs Bass model")
curve(logis(x,exp(fit_logis$par)),add=TRUE,col="tomato",lwd=2)
curve(Bass(x,exp(fit_bass$par)),add=TRUE,col="royalblue",lwd=2)
legend("bottomright",c("logistic","Bass"),lwd=2,lty=1,col=c("tomato","royalblue"))

残差平方和(だけ)をみると、バスモデルのほうが当てはまりがいい。

> fit_logis$value
[1] 100152.9
> fit_bass$value
[1] 87822.71

係数 b を見るとほとんど 0 に近い値になっている。

> round(exp(fit_bass$par),3)
[1] 10033.285     0.208     0.000     0.370

例2:カラーテレビ普及率のデータ

バスモデルのおもしろいところは他人にまどわされない購入者(innovator)と既購入者数 x_t がふえてくると乗り遅れまいとする購入者(imitator) の和になっているところである。

dx_t/dt は、

 \displaystyle dx_{t}/dt=a(m-x_t)+\frac{b}{m}x_{t}(m-x_t)

とふたつの項に分解できる。

係数 b がほぼ 0 だとあまりおもしろくないので、今度はカラーテレビの普及率のデータ(第1章第2節3 1 情報通信機器の世帯普及率 : 平成18年版 情報通信白書)にバスモデルを当てはめてみる。

x <- c(0.3, 1.6, 5.4, 13.9, 26.3, 42.3, 61.1, 75.8, 85.9, 90.3,
       93.7, 95.4, 97.7, 97.8, 98.2, 98.5, 98.9, 98.8, 99.2, 99.1,
       98.9, 98.7, 99.0, 99.3, 99.4, 99.3, 99.0, 99.1, 99.0, 98.9,
       99.1, 99.2, 99.2, 98.9, 99.0, 99.2, 99.3, 99.4, 99.0, 99.3, 99.4)/100
#
tim <- 1:length(x)
Bass <- function(t,pars){
  a <-pars[1]
  b <-pars[2]
  (1-exp(-(a+b)*t))/((b/a)*exp(-(a+b)*t)+1)
}

fit_bass <-optim(par=log(c(0.1,0.2)),
                  fn=function(pars)sum((x-Bass(tim,exp(pars)))^2),
                 control=list(maxit=2000))
fit_bass <-optim(par=fit_bass$par,
                 fn=function(pars)sum((x-Bass(tim,exp(pars)))^2),
                 control=list(maxit=2000),method = "BFGS")
plot(x)
curve(Bass(x,exp(fit_bass$par)),add=TRUE,col="royalblue",lwd=2)

f:id:abrahamcow:20160525032449p:plain

差分系列と dx_t/dta(m-x_t)\frac{b}{m}x_{t}(m-x_t) を合わせてプロットしてみる。

dBass1 <- function(t,pars){
  a <-pars[1]
  b <-pars[2]
  a*(1-Bass(t,pars))
}
dBass2 <- function(t,pars){
  a <-pars[1]
  b <-pars[2]
  b*(1-Bass(t,pars))*Bass(t,pars)
}
plot(diff(x))
curve(dBass1(x,exp(fit_bass$par))+dBass2(x,exp(fit_bass$par)),lwd=2,add=TRUE,col="royalblue",n=500)
curve(dBass1(x,exp(fit_bass$par)),col="forestgreen",lwd=2,add=TRUE,n=500)
curve(dBass2(x,exp(fit_bass$par)),col="orange",lwd=2,add=TRUE,n=500)
legend("topright",c("Adopters","Innovators","Imitators"),lwd=2,lty=1,
       col=c("royalblue","forestgreen","orange"))

f:id:abrahamcow:20160525032743p:plain

ついでにロジスティックカーブとの比較。

logis <- function(t,pars){
  b <- pars[1]
  c <- pars[2]
  1/(1+b*exp(-c*t))
}
derrivedlogis <- function(t,pars){
  b <- pars[1]
  c <- pars[2]
  (b*c*exp(c*t))/((b+exp(c*t))^2)
}
fit_logis <-optim(par=log(c(111,0.3)),
                  fn=function(pars)sum((x-logis(tim,exp(pars)))^2),
                  control=list(maxit=2000))
fit_logis <-optim(par=fit_logis$par,
                  fn=function(pars)sum((x-logis(tim,exp(pars)))^2),
                  control=list(maxit=2000))
fit_bass$value
fit_logis$value
plot(x)
curve(Bass(x,exp(fit_bass$par)),add=TRUE,col="royalblue",lwd=2)
curve(logis(x,exp(fit_logis$par)),add=TRUE,col="tomato",lwd=2)
legend("topright",c("Adopters","Innovators","Imitators","Logistic"),lwd=2,lty=1,
       col=c("royalblue","tomato"))
plot(diff(x))
curve(dBass1(x,exp(fit_bass$par))+dBass2(x,exp(fit_bass$par)),lwd=2,add=TRUE,col="royalblue",n=500)
curve(dBass1(x,exp(fit_bass$par)),col="forestgreen",lwd=2,add=TRUE,n=500)
curve(dBass2(x,exp(fit_bass$par)),col="orange",lwd=2,add=TRUE,n=500)
curve(derrivedlogis(x,exp(fit_logis$par)),add=TRUE,col="tomato",lwd=2,n=500)
legend("topright",c("Adopters","Innovators","Imitators","Logistic"),lwd=2,lty=1,
       col=c("royalblue","forestgreen","orange","tomato"))

f:id:abrahamcow:20160525032900p:plain

修正バスモデル

a(m-x_t)\frac{b}{m}x_{t}(m-x_t) をそれぞれ innovator、imitator と解釈したが、 m-x_tt 時点でまだ購入してない潜在的消費者)の中にはイノベーターもイミテーターも含まれているはずである。

イノベーターとイミテーターは分離されていない、という批判が、

によってなされている。

せっかくなので上述の論文の修正バスモデルもカラーテレビ普及率のデータに当てはめてみる。

modBass <- function(t,pars){
  beta1 <-pars[1]
  beta2 <-pars[2]
  xi <- pars[3]
  A <- -xi/beta2
  (beta2*A+xi*exp(-(beta2+xi*beta1)*t))/
    (beta1*A-exp(-(beta2+xi*beta1)*t))
}
fit_modBass <-optim(par=log(c(0.1,0.2,0.3)),
                 fn=function(pars)sum((x-modBass(tim,exp(pars)))^2),
                 control=list(maxit=2000))
fit_modBass <-optim(par=fit_modBass$par,
                 fn=function(pars)sum((x-modBass(tim,exp(pars)))^2),
                 control=list(maxit=2000),method = "BFGS")
plot(x)
curve(Bass(x,exp(fit_bass$par)),add=TRUE,col="royalblue",lwd=2)
curve(modBass(x,exp(fit_modBass$par)),add=TRUE,col="tomato",lwd=2)
legend("bottomright",c("Bass","mod-Bass"),lwd=2,lty=1,
       col=c("royalblue","tomato"))

f:id:abrahamcow:20160525034206p:plain

lambda0 <- function(t,pars){
  beta1 <-pars[1]
  beta2 <-pars[2]
  xi <- pars[3]
  xi*lambda2(t,pars)
}
lambda2 <- function(t,pars){
  beta1 <-pars[1]
  beta2 <-pars[2]
  xi <- pars[3]
  -beta1*modBass(t,pars)+beta2
}
plot(diff(x))
curve(lambda0(x,exp(fit_modBass$par))+lambda2(x,exp(fit_modBass$par))*modBass(x,exp(fit_modBass$par)),lwd=2,add=TRUE,col="royalblue")
curve(lambda0(x,exp(fit_modBass$par)),col="forestgreen",lwd=2,add=TRUE)
curve(lambda2(x,exp(fit_modBass$par))*modBass(x,exp(fit_modBass$par)),col="orange",lwd=2,add=TRUE)
legend("topright",c("Adopters","Innovators","Imitators"),lwd=2,lty=1,
       col=c("royalblue","forestgreen","orange"))

f:id:abrahamcow:20160525034232p:plain