decompose by R
疑問
前回の最後の所でR言語を使って、血圧の年変動を調べた。最終的に４つのグラフが示され、それを觀たオイラーは、季節変動を取り出すにFFTを使ってると勝手に解釈してた。
でも後で考えるに、それはないだろうと思った。じゃ、自分がやるとしたらどうやる？ グラフだけじゃ最終結果しか解らないので、その前段階を見る事にする。
> bld <- read.csv("bld.csv", header = F)
colnames(bld) <- c("yyyymm","hi")
hi.ts <- ts(bld$hi, start=c(2012,1), frequency=12)
> res <- decompose(hi.ts)
> round(res$x, 1)
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2012 125.1 123.3 121.0 118.2 118.9 115.8 113.3 117.5 116.3 119.8 122.5 124.7
2013 124.0 119.8 116.8 117.0 118.6 120.6 118.9 119.4 118.8 115.2 119.4 119.8
2014 122.0 120.3 122.7 121.4 121.6 123.4 122.9 122.4 123.2 123.0 121.5 124.6
2015 127.7 125.4 125.8 124.0 126.4 128.3 126.1 130.1 128.1 125.8 128.5 123.8
2016 129.3 126.3 121.3 124.6 126.5 128.1 129.6 128.7 130.3 129.7 129.0 128.2
2017 127.5 128.9 128.9 127.1 129.5 133.5 129.7
> round(res$trend, 1)
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2012 NA NA NA NA NA NA 119.7 119.5 119.1 118.9 118.9 119.1
2013 119.5 119.8 120.0 119.9 119.6 119.2 119.0 118.9 119.2 119.6 119.9 120.1
2014 120.4 120.7 121.0 121.5 121.9 122.2 122.7 123.1 123.4 123.7 124.0 124.4
2015 124.7 125.2 125.7 126.0 126.4 126.7 126.7 126.8 126.7 126.5 126.6 126.5
2016 126.7 126.8 126.8 127.1 127.2 127.5 127.6 127.6 128.0 128.4 128.7 129.0
2017 129.2 NA NA NA NA NA NA
> round(res$seasonal, 2)
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2012 2.10 -0.05 -1.64 -1.77 -0.41 1.30 -0.84 0.56 0.17 -0.62 0.69 0.50
2013 2.10 -0.05 -1.64 -1.77 -0.41 1.30 -0.84 0.56 0.17 -0.62 0.69 0.50
2014 2.10 -0.05 -1.64 -1.77 -0.41 1.30 -0.84 0.56 0.17 -0.62 0.69 0.50
2015 2.10 -0.05 -1.64 -1.77 -0.41 1.30 -0.84 0.56 0.17 -0.62 0.69 0.50
2016 2.10 -0.05 -1.64 -1.77 -0.41 1.30 -0.84 0.56 0.17 -0.62 0.69 0.50
2017 2.10 -0.05 -1.64 -1.77 -0.41 1.30 -0.84
> round(res$random, 2)
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2012 NA NA NA NA NA NA -5.52 -2.48 -2.98 1.50 2.94 5.12
2013 2.44 0.03 -1.54 -1.09 -0.55 0.09 0.82 -0.03 -0.53 -3.73 -1.17 -0.85
2014 -0.54 -0.36 3.31 1.62 0.09 -0.11 1.09 -1.25 -0.39 -0.10 -3.17 -0.32
2015 0.91 0.29 1.71 -0.30 0.39 0.27 0.18 2.70 1.28 -0.13 1.22 -3.21
2016 0.54 -0.39 -3.91 -0.66 -0.35 -0.68 2.90 0.53 2.09 1.93 -0.35 -1.28
2017 -3.88 NA NA NA NA NA NA
> round(res$figure, 2)
[1] 2.10 -0.05 -1.64 -1.77 -0.41 1.30 -0.84 0.56 0.17 -0.62 0.69 0.50
$type
[1] "additive"
attr(,"class")
[1] "decomposed.ts"
元データは規模を少し縮小して、201201 - 201607 までとした。
それから、結果の表示はroundを使って少数点以下を適当に丸めてしまった(てか、80桁に収まるように調整したってのが正解)。widthは標準でそうなってるけど、options(width=200) とかやって調整出来る。
４つのグラフの元データと$figureってのがつまっていた。 特徴的なのは、res$trend と res$random が、元データの前後６つ分のデータがNAになっている事。これって、データ無しを意味する。
ここから考察するに、元データ $x をフィルタリングってか平均化して、滑からtrend を作っていると予想される。
この平均化されたデータを電気回路で言う直流とみなせば、実波形は、それに１年周期の交流成分(と雑音)が重畳してるとみなせる。
ならば、実データから平均化データを引き算すれば、交流成分が取り出せる。後は同位相(平たく言えば月)同士を加算して、それの平均を求めていく。そうすれば、雑音成分が相殺される。 かくして、周期成分seasonalが求まる。これは、全データでも不変なはずだから、NAとする必要はない。
後は雑音成分の抽出。これは、random = trend - seasonal で求まる。trendにNAが含まれているように当然NAも生じてくる。なお、$figureは、seasonalの１周期分だ。
こんな推測をしたけど、どうよ。
トレンドを求める時、時系列データと移動平均 とか、移動平均(by wiki) を使ってるぽい。心情的に、現在の値を求めるのに未来のデータも使うって納得出来無いぞ。
src
推測を確認する為、ソースをお取り寄せ。これがOSSの楽しい所だ。見なきゃもったいない。
[sakae@fb /tmp/R-4.2.1]$ emacs ./src/library/stats/R/HoltWinters.R
# decompose additive/multiplicative series into trend/seasonal figures/noise
decompose <-
function (x, type = c("additive", "multiplicative"), filter = NULL)
{
type <- match.arg(type)
l <- length(x)
f <- frequency(x)
if (f <= 1 || length(na.omit(x)) < 2 * f)
stop("time series has no or less than 2 periods")
:
## return values
structure(list(x = x,
seasonal = seasonal,
trend = trend,
random = if (type == "additive")
x - seasonal - trend
else
x / seasonal / trend,
figure = figure,
type = type),
class = "decomposed.ts")
}
なんだい、大事な部分を端折っているぞ。
help
何はなくともmanてかhelpに当たれ、です。
> help(decompose)
:
The additive model used is:
Y[t] = T[t] + S[t] + e[t]
The multiplicative model used is:
Y[t] = T[t] * S[t] * e[t]
The function first determines the trend component using a moving
average (if ‘filter’ is ‘NULL’, a symmetric window with equal
weights is used), and removes it from the time series. Then, the
seasonal figure is computed by averaging, for each time unit, over
all periods. The seasonal figure is then centered. Finally, the
error component is determined by removing trend and seasonal
figure (recycled as needed) from the original time series.
This only works well if ‘x’ covers an integer number of complete
periods.
Examples:
require(graphics)
m <- decompose(co2)
m$figure
plot(m)
## example taken from Kendall/Stuart
x <- c(-50, 175, 149, 214, 247, 237, 225, 329, 729, 809,
530, 489, 540, 457, 195, 176, 337, 239, 128, 102, 232, 429, 3,
98, 43, -141, -77, -13, 125, 361, -45, 184)
x <- ts(x, start = c(1951, 1), end = c(1958, 4), frequency = 4)
m <- decompose(x)
## seasonal figure: 6.25, 8.62, -8.84, -6.03
round(decompose(x)$figure / 10, 2)
なんとなく、分ったような解らないような。。。
それは同然の事。ABC予測を証明した望月先生によれば、数は、足し算でも表現出来るし、かけ算でも表せる。圧倒的に難しいのは、足し算らしい。ええ、凡人は、任意の数は、数の元である素数の合成(かけ算)で表せるってのを思い出すのがやっとですから。
by python
こういう時は別の舞台で調べてみる。ユーザーの多い言語ね。
Python: statsmodels で時系列データを基本成分に分解する
(Python編) 時系列データをサクッとSTLでトレンド・季節性に分解
なるほど。でも、使うだけで、ソースにあたる人はいない。もったいないお化けが でてきますよ。
使い倒す
複雑な波形を分解する時、波形の足し算になってるから、それを分解してねって方法がある。これは馴染だだ。そう、FFTだ。オイラーが脊髄反応した奴。
もうひとつは、かけ算になってるからねってのがある。振幅変調された波形は、被変調波を搬送波とのかけ算で出来上がっているってやつ。
かけ算方式での分解、フィルターの設定も追加。
pdf("mul.pdf")
plot(decompose(hi.ts, type="multiplicative"))
dev.off()
xx = c(0.25, 0.25, 0.25, 0.25)
pdf("mulxx.pdf")
plot(decompose(hi.ts, type="multiplicative", filter=xx))
dev.off()
肝は、ソースの中でフィルターを適用して、トレンドを抽出する部分かな。
## filter out seasonal components
if (is.null(filter))
filter <- if (!f %% 2)
c(0.5, rep_len(1, f - 1), 0.5) / f
else
rep_len(1, f) / f
trend <- filter(x, filter)
今回の例だと、f = 12 すなわち１年を１２に分割、それって月だな。
rep_len(x,n) は、xをn回繰り返すって関数。
> rep_len(2,5) [1] 2 2 2 2 2
Ｒのデータ型、データ構造 これを読むとoctaveを思い出すぞ。
> x = 1:10 > x [1] 1 2 3 4 5 6 7 8 9 10 > filter(x, rep(0.2, 5)) Time Series: Start = 1 End = 10 Frequency = 1 [1] NA NA 3 4 5 6 7 8 NA NA > filter(x, rep(0.2, 5), sides=1) Time Series: Start = 1 End = 10 Frequency = 1 [1] NA NA NA NA 3 4 5 6 7 8
平均を求める時、何も指定しないと、注目してるポイントを中心に、過去と未来の値が使われる。これって気持悪い。と、思うのはオイラーだけじゃないみたいで、そんな人の為にsides=1を指定すると、過去のデータだけを使うようになる。
高血圧とか糖尿病は生活習慣病と言うがごとく、過去の不摂生の積み重ね(運動不足とか酒の飲み過ぎとか)で、現在の状態になってる。だから、過去のデータを勘案するのが正しい態度と思える。
上では、フィルターの係数を同一にしたけど、より過去の影響は少く、現時点に近付くにつれて係数を大きくするって事も出来る。例えば、c(0.1, 0.2, 0.3, 0.4)
また、高機能な stl と言う、時系列分析関数も提供されてる。どれを使うか迷うな。
filter
フィルターは、電気屋さんならみんな大好きってか、何かしらお世話になってる。たとえば、電源の平滑回路。ACを整流してDCにした時、その交流分を平滑。巨大な電解コンデンサーと、チョークコイルの回路ね。今だと物理サイズを小くする為、ACを一旦高い周波数に変換してから、整流してる。これだど、コンデンサーもコイルも物理的に小型に出来る。
オーデォオ屋さんなら、イコライザーとか言うフィルターね。無線屋さんなら、バンドパスフィルターとかクリスタル・フィルターとかだ。
で、デジタルなフィルターはどうなってる？ ちょいとRのソースを当る。 src/library/stats/R/filter.R
if(method == "convolution") {
if(nfilt > n) stop("'filter' is longer than time series")
sides <- as.integer(sides)
if(is.na(sides) || (sides != 1L && sides != 2L))
stop("argument 'sides' must be 1 or 2")
circular <- as.logical(circular)
if (is.na(circular)) stop("'circular' must be logical and not NA")
if (is.matrix(x)) {
y <- matrix(NA, n, nser)
for (i in seq_len(nser))
y[, i] <- .Call(C_cfilter, x[, i], filter, sides, circular)
} else
y <- .Call(C_cfilter, x, filter, sides, circular)
} else {
スピードが必要になるんで、核心部分はC言語だ。それを追跡すると、 src/library/stats/src/filter.c
SEXP cfilter(SEXP sx, SEXP sfilter, SEXP ssides, SEXP scircular)
{
:
R_xlen_t nx = XLENGTH(sx), nf = XLENGTH(sfilter);
if(sides == 2) nshift = nf /2; else nshift = 0;
if(!circular) {
for(i = 0; i < nx; i++) {
z = 0;
if(i + nshift - (nf - 1) < 0 || i + nshift >= nx) {
out[i] = NA_REAL;
continue;
}
for(j = max(0, nshift + i - nx); j < min(nf, i + nshift + 1) ; j++) {
tmp = x[i + nshift - j];
if(my_isok(tmp)) z += filter[j] * tmp;
else { out[i] = NA_REAL; goto bad; }
}
out[i] = z;
bad:
continue;
}
scheme風な型宣言が出てきたぞ。フィルターは、かけ算と足し算で出来ている。