数理人口学 稲葉寿先生 東京大学出版

#p.228 7.1基本的なＳＩＲモデル　－－学習中
#SIRモデルを最小２乗法で解く

S=c(21,16,13,11,12,10,9)
I=c(0,3,5,3,5,6,2)
R=c(0,2,3,7,4,5,10)
N=21

S=S[length(S):1];I=I[length(I):1];R=R[length(R):1]

t=c(1,2,3,4,5,6)

mu=1
beta=1
r=1

SX=S[1:(length(S)-1)]
IX=I[1:(length(I)-1)]
RX=R[1:(length(R)-1)]

SY=S[2:(length(S))]
IY=I[2:(length(I))]
RY=R[2:(length(R))]

ite=10000
eta=0.00001

for(l in 1:ite){

b=mean(SY+SX*(mu+beta-1))  
dmu=2*sum(SX*(SY-(b-SX*(mu+beta-1))))+2*sum(RX*(RY-(-(mu-1)*RX+r*IX)))+2*sum(IX*(IY-IX*(SX*beta-(mu+r-1))))  
dbeta=2*sum(SX*(SY-(b-SX*(mu+beta-1))))+2*sum((-SX*IX)*(IY-IX*(SX*beta-(mu+r-1))))
dr=2*sum(-IX*(RY-(-(mu-1)*RX+r*IX)))+2*sum(IX*(IY-(SX*IX*beta-(mu+r-1)*IX)))

mu=mu-eta*dmu
beta=beta-eta*dbeta
r=r-eta*dr

cost=sum((SY-(b-SX*(mu+beta-1)))^2)+sum((IY-IX*(SX*beta-(mu+r-1)))^2)+sum((RY-(-(mu-1)*RX+r*IX))^2)

print(cost)

}

R0=beta*N/(r+mu)

#定常解E1(感染者のいない状態：DFSS)
print(paste0("Sが",round(b/mu,2),"でIが",0))

#定常解E2(伝染病が人口に定着して共存するエンデミックな状態：ESS)
print(paste0("Sが",round(N/R0,2),"でIが",round(mu*(R0-1)/beta,2)))



init=c(19,2,0)