R言語所蔵のサンプルデータセットをABC順(大文字)優先で拝見させていただいとります。前回、前々回と太陽黒点関係のデータセットが続いています。正直飽きました。しかし今回は太陽黒点らしい解析など絶無。だただた時系列データを「いじる」ときに必要な小ネタ、TIPSの特集という感じです。勉強になったなあ。ホントか?
※「データのお砂場」投稿順Indexはこちら
sunspot.year
今回練習させていただきますのは、sunspots「3兄弟」の中で一番短いデータセットである sunspot.yearです。サンプルデータセットの解説ページが以下に
Yearly Sunspot Data, 1700–1988
本解説ページには処理例が掲載されております。既に太陽黒点周期みたいなものは前々回に少し観察しているので、今回は処理例(時系列データの小ネタ中心)をそのまま練習するの回です。なおいろいろ「テク」を御披露いただくのにsunspot.yearだけでは足りず、sunspotsも併用してます。
生データの確認
従前から data() でサンプルデータセットをロードして、str() でその構造を確認するというのを定石としてきました。ここの処理例も同じデス。
ただね、str呼び出すのに utils::strと指定してます。上記では処理例とは異なり utils:strの後で「何もなし」 str を使ってみてますが、特に違いは無いようです。環境(当方はRStudio使用)やバージョンの違いなんでしょうか?書いても問題ないし。
ともあれ、ターゲットのデータセットはいずれもTime-Seriesデータです。
時系列データの開始、終了
2本の時系列の開始と終了がズレているので、その重なっている部分を求めるための準備だと思うのです。処理例では以下のようにして、2系列のデータの開始終了を比べながら、t1に開始年と月、t2に終了年と月を求めてます。
よく見るとひっかかるのは、開始はstart(sm)などとしているのに、終了の方はend(sm)[1]などと第1要素(Rは1始まりデス)をことさら指定していることです。なんでじゃ?
実際にやってみるとend も start も返り値の形式は同じで、sunspotsデータのような年+月次データの場合、年 月 の2要素のベクトルを返してくるみたいです。sunspot.yearのような年だけの場合は年だけ1要素です。
それに対して max とると、ベクトルの長さの差などは踏みつぶして、上記のように年の最大値のみを取り出すことが可能です(月は1から12しかないもんね。)一方同じやり方で min とると、1月の1が一番小さいので、1などという答えが返ってきます。お間抜け。それでminの方はことさらに[1]などと第1要素を指定しているみたいです。正確を期するならば常に[1]書いていろよと思いますケド。
なお、(t1 <- …みたいに丸かっこで式全体を括っているのは、代入した結果を返してもらう(ここでは画面に表示)ためです。丸かっこしないと以下のように代入されても値は返ってきません。この辺はR言語のお作法の小ネタですな。
時系列データをwindow?
さて、2系列のデータの重なっている期間の開始、終了をt1, t2に確保できたので、期間をそろえた時系列データを抽出するステップです。使用するのは window()関数、これに元の時系列データと開始、終了を指定すれば良いと。
ただ結果をみると、monthの方は1984までになっているのにyearの方は1983までになってました。それにyearの方のend=t2[1]の[1]ってどゆこと?
t2は1983年12月だけれども、yearの方は1983年というのは1月相当にみなされ?ていて12月じゃないみたい。それにたいしてmonthの方はしっかり1983年12月まで入っているから?よくわかりませんな。
期間を揃えたデータの一致、不一致
その次は、上記で期間がそろった(といっても年次と月次の違いがあり)データを比較してデータが一致するかしないかを見たいということみたいです。stopifnot()関数で比較結果を判定してます。
問題(差異の指定範囲超え)なければ、上記のように何もかえってきませぬ。素気ないです。しかし、上記にはテク満載感あり。月次データの平均を取り、そしてあるtoleranceを指定してそのデータを年次データと比較(許容範囲内なら等しいとみなす)してます。かっこいいデス。しかし、toleranceで指定している0.0020などという値、どっから来たの?
そのお答えは以下のaggregateの計算にありました。多分、先にこちらをやってみて、toleranceを0.0020にすれば「イコール」になると分かってから上記の関数適用しているんじゃないかと。
しかし、重み付きの平均をとるのに、月毎の日数のベクタを作成しているところも何気にTIPS。特に2月の28.25は思いつかないデス。数百年分のデータの有効数字3桁くらいならこれでOKよな。