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

廿TT

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

めったに起きないコンバージョンの成長の非定常ポアソン過程によるモデル

R 微分方程式 確率過程

弱小アフィブログの成長

これは当ブログのアマゾンアフィリエイトレポートの折れ線グラフです。

f:id:abrahamcow:20151114044249p:plain

横軸が日付, 縦軸がその日の注文数合計です。

注文数合計は 0 が多く, 折れ線グラフは下に潰れて, 注文があった日だけスパイクのようにとげとげつきだしています。

サイトへの訪問者に対して, CV(注文)の発生が少ないため, ポアソン分布でモデル化できそうです。

ただし2014年の後半以降をみるとスパイクの間隔が短くなっており, パラメータ λ が一定のポアソン過程を当てはめるのは無理がありそうです。

そこでポアソン過程の強度(レートパラメータ)が成長していくようなモデルを考えることにします。

成長率が

  1. これまでに成長した分
  2. 成長の限界に達するまでに残された伸びしろ

のおのおのに比例すると, 成長曲線はロジスティック曲線になります。

時刻 t における強度が  \lambda(t) の非定常ポアソン過程では, ある時点における CV の発生しやすさが  \lambda(t), そのプロセスを T まで観測したときの CV の発生数が \Lambda(T) = \int_{0}^{T}\lambda\left( u\right) \,du. になります。

尤度関数は以下のように書けます。

\displaystyle L= \prod_{i=1}^{n} \lambda ( t_{i}) \exp\left(-\Lambda(T)\right)

\Lambda(T) にロジスティック曲線をいれた非定常ポアソン過程が, 今回提案するモデルです。

R で当てはめる

データは report.txt · GitHub に置いておきます。

日付を 2014年01月01日からの経過日数になおしてやり, CV の発生した日を抜き出したのが cvday_all です。

loglik は対数尤度です。その中の lambda はロジスティック関数を微分して対数をとったものです。

dayly <-read.table("~/Downloads/report.txt",skip=2,stringsAsFactors = FALSE)
head(dayly)
plot(V4~as.Date(V1),data=dayly,type="l",xlab="date",ylab="cv")

cvday <-as.Date(dayly$V1[dayly$V4>1])
cvday_all <- as.integer(cvday - as.Date("2014/01/01"))

loglik <- function(par){
  A <-par[1]
  B <-par[2]
  C <- par[3]
  maxT <- as.integer(as.Date("2015/10/31") - as.Date("2014/01/01"))
  lambda <- log(A)+log(B)+log(C)+B*cvday_all-2*log(exp(B*cvday_all)+C)
  Lambda <- A/(1+C*exp(-B*maxT))
  sum(lambda) - Lambda
}

fit1 <-optim(par=c(1, 1, 1),loglik,control = list(fnscale=-1),hessian = TRUE)
y <-1:length(cvday_all)
plot(cvday_all,y, ylab="cumulative CVs")
mylogis<- function(x,A,B,C){
  A/(1+C*exp(-B*x))
}
curve(mylogis(x,fit1$par[1],fit1$par[2],fit1$par[3]),col="red3",lwd=3,add=TRUE)

最尤法でパラメータを推定し, \Lambda(T) を赤線で累積の CV 発生回数に重ねてプロットしたのが下の図です。

f:id:abrahamcow:20151114052310p:plain

よく当てはまっているんじゃないでしょうか。

今後の課題

ここで「それってふつうにロジスティックカーブを当てはめるのとなにが違うの?」という疑問が生じます。

実は、ふつうにロジスティック曲線を当てはめるのと, 当てはまりはたいして変わりません。

fit2 <-nls(y~mylogis(cvday_all,A,B,C),start=list(A=83,B=0.01,C=33))
summary(fit2)
pars<- coef(fit2)
curve(mylogis(x,pars[1],pars[2],pars[3]),col="royalblue",lwd=3,add=TRUE)

f:id:abrahamcow:20151114053333p:plain

青い線が最小二乗法で当てはめたロジスティック曲線です。

非定常ポアソンで考えるとなにか良いことがあるか, と改めて考えると, 特になにも思いつきません。

パラメータの信頼区間はだいぶ変わるみたいです。

v1 <-sqrt(diag(solve(-fit1$hessian)))
upper<-fit1$par+qnorm(0.975)*v1
lower<-fit1$par+qnorm(0.025)*v1
> cbind(lower,upper)
           lower        upper
[1,] 64.03350911 102.15447693
[2,]  0.00715367   0.01148647
[3,]  7.30375362  58.14519290
> confint(fit2)
Waiting for profiling to be done...
          2.5%        97.5%
A 82.106603560 87.811601461
B  0.008065798  0.008940331
C 22.169892489 27.357969741