廿TT

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

ネットスラングの流行と衰退を追う微分方程式

今日の川柳

これは orz というネットスラングのグーグルトレンドです。

2004年ごろから急激に流行って、だんだん人気が衰えている。この動きの背景にどんなメカニズムがあるか考えます。

変数 S をネットスラングを使う人のポピュレーションとします。
変数 I をまだネットスラングを使ってないが、これから使う可能性のある人のポピュレーションとします。

S と I が接触すると、病気が感染するかのように I の人が S の人に変化すると考えます。
また一方で、ネットスラングは一定の割合でだんだん忘却されていくと考えます。

これをかんたんな微分方程式で記述すると以下のようになります。

 dS(t)/dt = \beta S(t) I(t)  - \alpha S(t)
 dI(t)/dt = -\beta S(t) I(t)

ここで t は時刻です。(ってこれSIRモデルと一緒やん。まじで後から気づいた。)

\betaネットスラングの感染に関するパラメータで、\alpha は忘却に関するパラメータです。

R を使って最小二乗法で曲線 S をグーグルトレンドのデータに当てはめてみます。

f:id:abrahamcow:20190103222027p:plain

こんな感じです。

推定されたパラメータとS、Iの初期値は以下のようになりました。

> print(exp(opt$par))
     alpha       beta          S          I 
0.01636698 1.60774059 0.09456285 0.70956493 

R のコードを貼ります。やはり R はいい……。

library(deSolve)
library(tidyverse)
library(cowplot)
dat <-read.csv("~/Downloads/multiTimeline.csv",skip=1)
dat <- dat %>% 
  set_names(c("time","interest")) %>% 
  mutate(time=as.Date(paste0(time,"-01")),interest=interest/100)

SImodel <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    dS <- beta*S*I - alpha*S
    dI <- -beta*S*I
    list(c(dS, dI))
  })
}

Nt <- nrow(dat)

times <- 1:Nt

costfunc <- function(par,y){
  parameters = exp(par[1:2])
  state <- exp(par[3:4])
  ode_out <- ode(y=state, times=times, func=SImodel, parms=parameters)
  sum((y-ode_out[,"S"])^2)
}

opt <- optim(c(alpha=-2,beta=-1,S=-2,I=0),costfunc,y=dat$interest)

print(exp(opt$par))

parameters = exp(opt$par[1:2])
state <- exp(opt$par[3:4])
state <- c(S=plogis(opt$par[3]),I=plogis(opt$par[3],lower.tail = TRUE))
ode_out <- ode(y=state, times=times, func=SImodel, parms=parameters)

ggplot(dat,aes(x=time,y=interest))+
  geom_line(aes(y=ode_out[,"S"],colour="S"),size=1)+
  geom_line(aes(y=ode_out[,"I"],colour="I"),size=1)+
  geom_point()+
  theme(legend.title = element_blank())

微分方程式で数学モデルを作ろう

微分方程式で数学モデルを作ろう

catindog.hatenablog.com