廿TT

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

変化点のあるポアソン過程へのパラメトリックな曲線あてはめ

今日の川柳

f:id:abrahamcow:20200219021313p:plain

Lambda <- function(x,a,b,c){a*x-b*log1p(exp(-x+c))+b*log1p(exp(c))}
lambda <- function(x,a,b,c){a+b/(1+exp(x-c))}
ll <- function(par,y,Tau){
  a <- exp(par[1])
  b <- exp(par[2])
  c <- exp(par[3])
  Lambda(Tau,a,b,c)-sum(log(lambda(y,a,b,c)))
}

dll <- function(par,y,Tau){
  a <- exp(par[1])
  b <- exp(par[2])
  c <- exp(par[3])
  lam <- lambda(y,a,b,c)
  da <- a*(Tau-sum(1/lam))
  db <- b*(-log1p(exp(-(Tau-c)))+log1p(exp(c))-
             sum(1/(1+exp(y-c))/lam))
  dc <- c*(-b/(1+exp(Tau-c))+b/(1+exp(-c))-
             sum((b*exp((y-c))/(1+exp(y-c))^2)/lam))
  c(da,db,dc)
}

mind <- min(coal$date)
Tau <- max(coal$date)-mind
opt <- optim(c(1,0,3),ll,y=coal$date-mind,Tau=Tau, hessian = TRUE)

plot(coal$date-mind,1:length(coal$date),xlab = "date",ylab = "cumulative counts",type="s",xaxt="n")
curve(Lambda(x,exp(opt$par[1]),exp(opt$par[2]),exp(opt$par[3])),add=TRUE,lty=2,col="royalblue",lwd=2)
axis(side=1,at=coal$date-mind,labels = round(coal$date))
abline(v=exp(opt$par[3]),lty=3)

点過程の時系列解析 (統計学One Point)

点過程の時系列解析 (統計学One Point)

追記

ちょっと一般化した。

Lambda <- function(x,a,b,c,d){a*x+b*c*log1p(exp(c*(x-d)))-c*b*log1p(exp(-c*d))}
lambda <- function(x,a,b,c,d){a+b*c^2/(1+exp(-c*(x-d)))}


ll <- function(par,y,Tau){
  a <- exp(par[1])
  b <- exp(par[2])
  c <- par[3]
  d <- exp(par[4])
  Lambda(Tau,a,b,c,d)-sum(log(lambda(y,a,b,c,d)))
}


dll <- function(par,y,Tau){
  a <- exp(par[1])
  b <- exp(par[2])
  c <- par[3]
  d <- exp(par[4])
  lam <- lambda(y,a,b,c,d)
  da <- a*(Tau-sum(1/lam))
  db <- b*(c*log1p(exp(c*(Tau-d)))-c*log1p(exp(-c*d))-
    sum(c^2/(1+exp(-c*(y-d)))/lam))
  dc <- (b*log1p(exp(c*(Tau-d))) + b*c*(Tau-d)/(1+exp(-c*(Tau-d)))-b*log1p(exp(-c*d)) + d*c*b/(1+exp(c*d))-
    sum((2*b*c/(1+exp(-c*(y-d)))+b*c^2*(y-d)*exp(-c*(y-d))/(1+exp(-c*(y-d)))^2)/lam))
  dd <- d*((-b*c^2/(1+exp(-c*(Tau-d)))+b*c^2/(1+exp(c*d)))-
    sum(-((b*c^3)*exp(-c*(y-d))/(1+exp(-c*(y-d)))^2)/lam))
  c(da,db,dc,dd)
}

# p <- runif(4,1,2)
# x <- sort(runif(10))
# numDeriv::grad(ll,p,y=x,Tau=1)
# dll(p,x,Tau=1)

n <- nrow(coal)
mind <- min(coal$date)
Tau <- max(coal$date)-mind
opt <- optim(c(1,1,-1,3),ll,dll,y=coal$date-mind,Tau=Tau, method = "BFGS")
opt$par
plot(coal$date-mind,1:length(coal$date),xlab = "date",ylab = "cumulative counts",type="s",xaxt="n")
curve(Lambda(x,exp(opt$par[1]),exp(opt$par[2]),opt$par[3],exp(opt$par[4])),
      add=TRUE,lty=2,col="royalblue",lwd=2)
axis(side=1,at=coal$date-mind,labels = round(coal$date))
abline(v=exp(opt$par[4]),lty=3)