廿TT

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

ロトカ・ヴォルテラの方程式のパラメータ推定

ロトカ・ヴォルテラの方程式については、

abrahamcow.hatenablog.com

を見てください。

hare–lynx データ

ここに一組のデータ・セットがあります。

Hare=c(30, 47.2, 70.2, 77.4, 36.3, 20.6, 18.1, 21.4, 22, 25.4, 27.1, 40.3, 57, 76.6,
       52.3, 19.5, 11.2, 7.6, 14.6, 16.2, 24.7)
Lynx=c(4, 6.1, 9.8, 35.2, 59.4, 41.7, 19, 13, 8.3, 9.1, 7.4, 8, 12.3, 19.5,
       45.7, 51.1, 29.7, 15.8, 9.7, 10.1, 8.6)
Year = 1900:1920
plot(Year,Hare,type="b",ylim=c(0,max(Hare)),col="red",lwd=2,ylab="Population")
lines(Year,Lynx,type="b",pch=4,col="blue",lwd=2)
legend("topright",c("hare","lynx"),col=c("red","blue"),lty=1,lwd=2,pch=c(1,4))

f:id:abrahamcow:20160609030905p:plain

これは取引された毛皮の数からウサギとヤマネコの個体数を推定したデータで、hare はウサギ、lynx はヤマネコです。

このデータはロトカ・ヴォルテラの方程式がよく適合する例としてしられています。

パラメータ推定

データがあり、モデルもあるのであとは当てはめるだけですがこれがなかなかうまくいかない。

微分方程式を数値的に解きつつ、最小二乗法でパラメータを推定するコードは以下です。

library(deSolve)
LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dx = (a-k*y)*x
    dy = (-b+l*x)*y
    return(list(c(dx,dy)))
  })
}
objfun <- function(pars){
  out1   <- ode(ini,1:n, LVmod, exp(pars))
  sum((Hare-out1[,"x"])^2) + sum((Lynx-out1[,"y"])^2)
}
n <- length(Lynx)
ini=c(x=Hare[1],y=Lynx[1])
res <-optim(log(c(a=0.4,b=0.018,k=0.8,l=0.023)),objfun,control = list(maxit=2000))
pars <- exp(res$par)
tim <- seq(1,n,by=0.1)
out1   <- ode(ini,tim, LVmod, pars)
Year2=seq(1900,1920,by=0.1)
plot(Year,Hare,ylim=c(0,max(Hare)))
lines(Year2,out1[,"x"])
plot(Year,Lynx,ylim = c(0,max(Lynx)))
lines(Year2,out1[,"y"])

当てはめプロットをみるとバラバラです。

f:id:abrahamcow:20160609032456p:plain

f:id:abrahamcow:20160609040604p:plain

中心差分近似

そこでロトカ・ヴォルテラを少し変形してみなおすと、微分形は単純な一次式になっていることに気づきます。

 \displaystyle \frac{1}{x}\frac{dx}{dt} = a - k y

 \displaystyle \frac{1}{y}\frac{dy}{dt} = -b + l x

離散的なデータは微分することができないので差分で近似することを考えます。

中心差分近似は、

 \displaystyle f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}

の式で与えられます。

こうして変形したデータに対して単回帰で係数を推定してやればパラメータが求まります。

dydty=(1/Lynx[2:(n-1)])*((Lynx[3:n]-Lynx[1:(n-2)])/2)
dxdtx=(1/Hare[2:(n-1)])*((Hare[3:n]-Hare[1:(n-2)])/2)

fit1 <-lm(dydty~Hare[2:(n-1)])
fit2 <-lm(dxdtx~Lynx[2:(n-1)])

b = unname(-fit1$coefficients[1])
l = unname(fit1$coefficients[2])
a = unname(fit2$coefficients[1])
k = unname(-fit2$coefficients[2])
pars <- c(a=a,k=k,b=b,l=l)
out1   <- ode(ini,tim, LVmod, pars)

plot(Year,Hare)
lines(Year2,out1[,"x"])
plot(Year,Lynx)
lines(Year2,out1[,"y"])

当てはめプロットはこちら。

f:id:abrahamcow:20160609035514p:plain

f:id:abrahamcow:20160609035522p:plain

こうして求めたパラメータを初期値にしてネルダー・ミード法を回すとさらにフィッティングがよくなります。

res <-optim(log(pars),objfun,control = list(maxit=2000))
pars <- exp(res$par)
out1   <- ode(ini,tim, LVmod, pars)

plot(Year,Hare)
lines(Year2,out1[,"x"])
plot(Year,Lynx)
lines(Year2,out1[,"y"])

f:id:abrahamcow:20160609110642p:plain

f:id:abrahamcow:20160609110651p:plain

参考文献

本エントリは全面的に、
http://www.math.tamu.edu/~glahodny/Math442/Curve%20Fitting.pdf
に依存しています。データもここからもらってきました。

あと、微分方程式の入門書として以下の本もおすすめです。

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

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

Hirsch・Smale・Devaney 力学系入門―微分方程式からカオスまで

Hirsch・Smale・Devaney 力学系入門―微分方程式からカオスまで