バスモデルのなんたるかについては バスモデル - ORWiki を参照。
バスモデルは以下の微分方程式で記述される。
閉じた形で解が求まる。
検算
deSolve パッケージを使って数値的に解いた値と解析解をくらべて、この解が正しいことを一応確かめた。
丸が数値解、曲線が解析解。
#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)
例１：インターネット利用者数
続いてインターネット利用者数のデータ（総務省｜平成25年版 情報通信白書｜インターネットの利用状況）にバスモデルとロジスティック曲線をそれぞれ当てはめてみた。
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
例２：カラーテレビ普及率のデータ
バスモデルのおもしろいところは他人にまどわされない購入者（innovator）と既購入者数 がふえてくると乗り遅れまいとする購入者（imitator) の和になっているところである。
は、
とふたつの項に分解できる。
係数 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)
差分系列と 、
、
を合わせてプロットしてみる。
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"))
ついでにロジスティックカーブとの比較。
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"))
修正バスモデル
、
をそれぞれ innovator、imitator と解釈したが、
（t 時点でまだ購入してない潜在的消費者）の中にはイノベーターもイミテーターも含まれているはずである。
イノベーターとイミテーターは分離されていない、という批判が、
によってなされている。
せっかくなので上述の論文の修正バスモデルもカラーテレビ普及率のデータに当てはめてみる。
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"))
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"))
