2015-02-14

三角関数の整数係数多項式表現

ARIMA 多項式三角関数

ここ数日ARIMAをやっている
そこからの脱線で、ちょっとおもしろいことに気が付いた。三角関数の級数表現のこと。
ARIMAのar部分は、現在と過去とに影響されて次時点が決まるような確率過程である
周期性を持っているかに見える時系列データができることもわかる
さて。ARIMAから離れよう
周期性と言えば、三角関数。
三角関数は、
- $\frac{dx}{dt}=-p y$
- $\frac{dy}{dt}=q x$
というような連立常微分方程式の解である。
- $x=A \cos{pq t},y=A \sin{pq t}$
これでは、全然「1変数のARIMA」ではないのだが
- $x(t+1) = x(t) -py(t) = x(t)-p(y(t-1)+qx(t-1)$ と離散化して書き、どんどん、このy項をx項で置き換えていくと
- $x(t+1) = x(t) - (pq(\sum_{i=1}^{t-1} x(i))) -py(1)$ となる
- $-py(1)$ は定数項だから無視することにして
- $x(t+1) = x(t) - (P \sum_{i=1}^{t-1} x(i))$ となる( $P=pq$ と置いた)
- これがどんな数列になるかと言うと、周期数列になる

P <- 0.01
n <- 300
x <- rep(0,n)
x[1] <- 1 # 初期値は全体を等倍にするだけなので1として実行する
for(i in 2:n){
	x[i] <- x[i-1] - P * sum(x[1:(i-1)])
}
plot(x,type="l")

さて、三角関数になるはずの微分方程式の解の定幅時刻の値を算出する再帰的計算コードは、次のように、各要素がPの整数係数多項式になる
- $x(1)=1$
- $x(2)=1-P$
- $x(3)=x(2)-P(x(1)+x(2))=1-P-P(1+1-P)=1-3P+P^2$ ...
この整数係数は以下のようにループで計算することができる
- 要するに、x(t+1)の係数はx(t)の係数と、x(1)...x(t)の係数の和x(-1)をPの次数を一つ上げた項の係数との和にしている

n <- 10
coefs <- matrix(0,n,n+1)
coefs[1,1] <- 1

for(i in 2:n){
	coefs[i,] <- coefs[i-1,]
	tmp <- apply(matrix(coefs[1:(i-1),],ncol=n+1),2,sum)
	coefs[i,2:(n+1)] <- coefs[i,2:(n+1)] - tmp[1:n]
}
coefs

> coefs
      [,1] [,2] [,3] [,4] [,5]  [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]    1    0    0    0    0     0    0    0    0     0     0
 [2,]    1   -1    0    0    0     0    0    0    0     0     0
 [3,]    1   -3    1    0    0     0    0    0    0     0     0
 [4,]    1   -6    5   -1    0     0    0    0    0     0     0
 [5,]    1  -10   15   -7    1     0    0    0    0     0     0
 [6,]    1  -15   35  -28    9    -1    0    0    0     0     0
 [7,]    1  -21   70  -84   45   -11    1    0    0     0     0
 [8,]    1  -28  126 -210  165   -66   13   -1    0     0     0
 [9,]    1  -36  210 -462  495  -286   91  -15    1     0     0
[10,]    1  -45  330 -924 1287 -1001  455 -120   17    -1     0

- これは、行について、その係数の並びを見てもなんだかよくわからないけれど、縦に見ると、たとえば、３列目は1,5,15,30,70,...となるが、これはa(n) = n*(n-1)*(n-2)*(n-3)/24　であって、A000332。
- これをC(n,4)と書けば、実は、この係数行列の列は、C(n,1),C(n,2),C(n,4),C(n,8),....となる
確かにこの整数係数が先ほどの漸化式的数列計算と合っていることを確認する。もちろん、多項式表現の方は係数が大きくなりP^kが小さく/大きくなりすぎるので、計算は不正確になるのだが、それは計算上の話だけなので、今は少なくとも気にしないことにする

P <- 0.01
n <- 10
x <- rep(0,n)
x[1] <- 1 # 初期値は全体を等倍にするだけなので1として実行する
for(i in 2:n){
	x[i] <- x[i-1] - P * sum(x[1:(i-1)])
}
qs <- P^(0:n)
p <- coefs %*% qs
p
x

> p
           [,1]
 [1,] 1.0000000
 [2,] 0.9900000
 [3,] 0.9701000
 [4,] 0.9404990
 [5,] 0.9014930
 [6,] 0.8534721
 [7,] 0.7969164
 [8,] 0.7323916
 [9,] 0.6605429
[10,] 0.5820888
> x
 [1] 1.0000000 0.9900000 0.9701000 0.9404990 0.9014930 0.8534721 0.7969164
 [8] 0.7323916 0.6605429 0.5820888

このようにして、何かしら、コサイン関数の定幅値が整数係数多項式となった
では、この整数係数多項式は、どんなコサイン関数なのだろうか
コサイン関数は $A \cos(\theta t + \phi)$ というように、スケール係数 $A$ 、周期の係数 $\theta$ 、位相の係数 $\phi$ で決まる
理論的に攻めることもできるのだろうが、こんな風にしてが次のようになるらしいことがわかった
- $A=\frac{2}{\sqrt{4-P}}$
- $\phi = acos{\frac{\sqrt{4-P}}{2}$
- $\theta = 2\phi$
P=1,2,3のときは、xの値が単純な繰り返しになることを利用し、それぞれ
- $\frac{2}{\sqrt{3}} \cos{\frac{2\pi}{6} + \frac{\pi}{6}$
- $\frac{2}{\sqrt{2}} \cos{\frac{2\pi}{4} + \frac{\pi}{4}$
- $\frac{2}{\sqrt{1}} \cos{\frac{2\pi}{3} + \frac{\pi}{3}$
と一致することは確認できる
- また、P=0の場合は、x=1の水平線であることもわかる。このことから、 $A=\frac{2}{\sqrt{4-P}$ らしいことはわかる
そうすると、x(1)=1であることより、この時刻を0とすれば、位相 $\phi$ は $\cos{\phi} = \frac{1}{A}=\frac{1}{\frac{2}{\sqrt{4-P}}}$ となるはずであり、それをあてはめると上記の $A,\theta,\phi$ の定義となる
やってみよう

Ps <- sort(c(1:3,runif(5)*4))

n <- 20
t <- 0:(n-1)
t2 <- seq(from=0,to=n-1,length=10*n)
par(mfrow=c(2,4))
for(i in 1:length(Ps)){
	P <- Ps[i]
	A <- 2/sqrt(4-P)
	phi <- 2^(3-P) + 2
	phi <- acos(sqrt(4-P)/2)
	theta <- 2*phi
	x <- 1:n
	x[1]<-1
	for(i in 2:n){
		x[i] <- x[i-1] + (-P) * sum(x[1:(i-1)])
	}
	y <- A*cos(t2*theta+phi)
	ylim <- range(c(x,y))

	plot(t,x,type="l",ylim=ylim)
	points(t2,y,type="l",col=2)	
}
par(mfcol=c(1,1))