tag:blogger.com,1999:blog-87410219643261807502024-03-14T02:26:15.540+09:00CYBERer.NET番町のITおじさんが書く、ITエンジニアの教養ブログです。様々なプログラミング言語やソフトウェア工学についての知識、ITやIT業界の歴史、動向などを取り上げます。管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comBlogger99125tag:blogger.com,1999:blog-8741021964326180750.post-56694553891442322092022-04-19T01:55:00.000+09:002022-04-19T01:55:32.327+09:00JPX J-Quants API の Early access NFT に当選しました<div>日本取引所グループで、<a href="https://jpx.gitbook.io/j-quants-api-faq">J-Quants API</a> という取り組みがスタートしました。</div><div>個人投資家に向けて、投資情報を API で入手できるようにする取り組みで、4/18 から Early access が始まりました。</div><div><br /></div><h2 style="text-align: left;">アクセスキーをNFTで配布</h2><div>この Early access は、事前に応募した人から抽選で当たった人だけが参加できるのですが、当たった人には Early access へのアクセスキーが NFT として配布されます。</div><div><br /></div><div>私がもらった NFT は、<a href="https://opensea.io/assets/matic/0x929da6bb317a268e6953c2ea0ae4fad9f0ae11a0/121" target="_blank">https://opensea.io/assets/matic/0x929da6bb317a268e6953c2ea0ae4fad9f0ae11a0/121</a> です。シリアルナンバー #121 ですね。</div><div><br /></div><div><a href="https://opensea.io/collection/j-quantsnft-accesskey">https://opensea.io/collection/j-quantsnft-accesskey</a> を見ると、発行された NFT が全部見れます。シリアルナンバーが #1〜#30 と #100〜#139 に分断されていて、#30 は Early access 開始前の 4/14 に発行されていて、#100 以降は 4/18 発行なので、#1〜#30 は JPX 内部の方がお持ちなんじゃないかな、と思いました。レア物ですね。</div><div><br /></div><div>現時点では #139 まで発行されてますが、当選した人が登録したときに mint されるとのことと、当選者は1名 invite できるとのことで、まだ invite も終わっていないと思うので、まだ増えるでしょう。100名くらいになるんでしょうか。</div><div><br /></div><div>Early access 終了後は NFT はアクセスキーとしての役割は終え、記念品(?) として残るそうです。おもしろい。</div><div><br /></div><h2 style="text-align: left;">API 試してみた</h2><div>NFT を使って J-Quants API のページにログインすると、アクセスキーを取得できます。24時間しか有効期間が無いし、自動更新できないんですが、まぁ検証目的なのでしょうがないですね。</div><div>翌日決算発表予定を取得する API があったので、Python で叩いてみました。</div><div><br /></div>
コード:
<pre class="prettyprint"><div><div>import requests</div><div>import json</div><div><br /></div><div>import accesskey</div><div><br /></div><div>headers = {</div><div><span style="white-space: pre;"> </span>'Authorization': 'Bearer {}'.format(accesskey.ACCESS_KEY)</div><div>}</div><div>r = requests.get("https://api.jpx-jquants.com/v1/fins/announcement", headers=headers)</div><div>data = r.json()</div><div>print(json.dumps(data, indent=4, ensure_ascii=False))</div></div></pre>
結果:
<pre class="prettyprint"><div><div>{</div><div> "Code": "24110",</div><div> "Date": "20220415",</div><div> "CompanyName": "ゲンダイエージェンシー",</div><div> "FiscalYear": "3月31日",</div><div> "SectorName": "情報・通信業",</div><div> "FiscalQuater": "本決算",</div><div> "Section": "スタンダード"</div><div>}</div></div></pre>SectorName, FiscalQuater (Quarter のことかな? 報告しておこう), Section が日本語なのはちょっとつらい気もしますが、それでも API で取れるのはいいですね。いろんな自動化の夢が拡がります。<div><br /></div><div>6月からはβ版利用者が一般募集になるとのことなので、それまでアーリーアダプタを謳歌しようと思います。</div><div><br /></div><div><br /></div>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-32405242107104542732021-08-27T02:57:00.007+09:002021-08-27T03:04:19.262+09:00 M1 mac (macOS / AArch64) 用 JDK 17 の Release-Candidate Builds を動かしてみた<div>Java 17 は2021年9月リリース予定となっています。先立って <a href="https://jdk.java.net/17/">JDK 17</a> の Release-Candidate Builds が 2021年8月6日に公開されました。</div><div>今回は、M1 mac (macOS / AArch64) 版を動かしてみたので、さらっと手順を紹介します。</div><div><br /></div><h4 style="text-align: left;"><span><a name='more'></a></span>余談</h4><div>以前、本ブログで <a href="https://www.cyberer.net/2020/04/jdk-14-GA-release.html">Java 14 の新機能まとめ</a> を書きました。その後 Java 15 以降も書こうと思っていたのですが、書きそびれていて、とうとう Java 17 リリース間近になってしまいました。</div><div>JDK 17 は Long-time support (LTS) となる予定なので、さすがにそろそろ Java 15〜17 の新機能をまとめて取り上げておきたいなと思い、サンプルコードを書くための簡単な開発環境を用意してみました。</div><div><br /></div><h2 style="text-align: left;">ダウンロード</h2><div><a href="https://jdk.java.net/17/">JDK 17 のページ</a> から macOS / AArch64 の tar.gz をダウンロードします。正しくダウンロード出来ているか確認するため sha256 もチェックしておくと良いですね。ダウンロードしたファイルの sha256 はターミナルから次のような感じで算出できます。</div><div><br /></div><div><div><pre class="prettyprint">$ shasum -a 256 ~/Downloads/openjdk-17_macos-aarch64_bin.tar.gz
b5bf6377aabdc935bd72b36c494e178b12186b0e1f4be50f35134daa33bda052 /Users/xxxxxxxx/Downloads/openjdk-17_macos-aarch64_bin.tar.gz</pre><br />
</div></div><h2 style="text-align: left;">tar.gz アーカイブの展開</h2><div>RC だからだと思うのですが、dmg ではなく tar.gz 配布なので、展開して手動で配置してあげる必要があります。</div><div>tar.gz の展開は、Finder からダウンロードした tar.gz をダブルクリックするだけです。jdk-17.jdk というフォルダが展開されると思います。</div><div><br /></div><h2 style="text-align: left;">手動で配置</h2><div>展開されて出てきた jdk-17.jdk を Finder で ライブラリ - Java - JavaVirtualMachines (パスで言うと /Library/Java/JavaVirtualMachines) に移動します。</div><div><br /></div><h2 style="text-align: left;">ターミナルから動作確認</h2><div><pre class="prettyprint">$ /usr/libexec/java_home -v 17 --exec java -version
openjdk version "17" 2021-09-14
OpenJDK Runtime Environment (build 17+35-2724)
OpenJDK 64-Bit Server VM (build 17+35-2724, mixed mode, sharing)</pre>
</div><div style="text-align: left;"><br /></div><h2 style="text-align: left;">vscode の設定</h2><div>まず Language Support for Java(TM) by Red Hat エクステンションをインストールします。</div><div>その後 Settings を開き、jdk などで検索すると、Java: Home というエントリがあると思うので、その Edit in settings.json をクリックしてください。settings.json の java.home の値の箇所が開きますので、以下のように設定します。先程配置した jdk-17.jdk の中の Contents/Home が Java home になります。</div><div><pre class="prettyprint">"java.home": "/Library/Java/JavaVirtualMachines/jdk-17.jdk/Contents/Home",</pre><br />
</div><h2 style="text-align: left;">vscode で Java Project を作ってみる</h2><div>Language Support for Java(TM) by Red Hat エクステンション が入っていれば、EXPLORER の Create Java Project を押すことで Java Project が作れます。配置先パスとプロジェクト名を設定すれば、OK です。</div><div><br /></div><h2 style="text-align: left;">vscode から実行してみる</h2><div>上記の要領で Java Project を作ると、src/App.java にいわゆる Hello, World! が自動生成されていると思います。そのファイルを開き、右上の ▷ から Run Java を押せば、TERMINAL 内で実行されます。概ね以下のように TERMINAL に出力されると思います。java コマンドのパスが jdk-17.jdk になっていれば OK です。</div><div><pre class="prettyprint">$ /usr/bin/env /Library/Java/JavaVirtualMachines/jdk-17.jdk/Contents/Home/bin/java -XX:+ShowCodeDetailsInExceptionMessages -Dfile.encoding=UTF-8 -cp /Users/xxxxxxxx/jdk17sandbox/jdk17sandbox/bin App
Hello, World!</pre><br />
</div><div><br /></div><div><br /></div>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-85717115179928155602021-07-31T00:00:00.026+09:002021-08-27T03:09:19.674+09:00第2期決算出しました<p>第2期の決算処理を終え、確定申告・納税も済ませました。</p><p>昨年は「第2期は税理士さん頼まないとちょっと厳しいかもしれません」と言ってましたが、今回も何とか自力でやりました。今期はお金の動きが多かったのでちょっと大変でした。</p><p>今回は、電子申告にチャレンジし、無事出来たので、いい経験値でした。ですが、電子納税は手続きが間に合わず、窓口での払込みになってしまいました。次は電子納税にチャレンジです。</p><p>引き続きよろしくお願いいたします。</p><p><br /></p>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-80280236378278664302020-10-13T03:11:00.002+09:002020-10-13T03:11:40.684+09:00Windows 10 機の接続モニタより高い解像度で Chrome リモートデスクトップを使いたい (NVIDIA 機向け)<p>我が家の Windows 10 機は、グラフィックカードが NVIDIA GeForce GTX 1050 です。接続しているモニタは16年物くらいのえらい古い物なんですが、最近壊れてしまって、電源を入れてから1秒くらいしか画面が映りません……。でも、MacBookPro から Chrome リモートデスクトップで接続して使っているため、特に困っておらず、買い替えに至っていません。</p><p>ただちょっと不満だったのが、解像度が 1280×1024 になってしまうことでした。Windows 純正のリモートデスクトップでは、クライアント側の環境に合わせて解像度が変更できますが、Chrome リモートデスクトップではモニタに映っている通りの解像度にしかなりません。接続されているモニタの解像度が 1280×1024 のため、Chrome リモートデスクトップでも 1280×1024 で使うしかないのです。クライアント側の MacBookPro は、通常 1440×900 という解像度なのですが、私は <a href="http://displaymenu.milchimgemuesefach.de/">Display Menu</a> というアプリを使って 2880×1800 まで解像度を上げています。そのため 1280×1024 というのはとても狭く感じ、不満だったのです。</p><p>ここ数ヶ月、どうにかして 2880×1800 の解像度で Windows 10 機に接続できないかといろいろ試していたのですが、やっとできたので、その手順を記しておきたいと思います。</p><span><a name='more'></a></span><h2 style="text-align: left;">前提となる環境の説明</h2><p></p><ul style="text-align: left;"><li>ホスト側</li><ul><li>Windows 10 機</li><li>Windows 10 Home Edition を使っているため、純正リモートデスクトップは使えません</li><li>NVIDIA GeForce GTX 1050 を搭載していて、NVIDIA ドライバが動いています</li><li>接続されているモニタは最大解像度が 1280×1024 です</li></ul><li>クライアント側</li><ul><li>MacBookPro</li><li>Display Menu が動いていて論理解像度が 2880×1800 になっています</li><li>ここから Chrome リモートデスクトップでホスト側の Windows 10 機に接続します</li></ul></ul><div><br /></div><h2 style="text-align: left;">手順</h2><p></p><h3 style="text-align: left;">Chrome リモートデスクトップで接続する</h3><div>物理モニタの解像度を超えた設定にしてしまうので、物理モニタ越しに作業していると、不都合が起きる可能性があります。Chrome リモートデスクトップで接続して作業します。</div><div><br /></div><h3 style="text-align: left;">カスタム解像度を作成する</h3><p>Chrome リモートデスクトップは、物理モニタに表示されている通りの解像度でしか表示できませんので、Windows 10 のディスプレイ設定で目的の解像度を選択する必要があります。しかし選択できる解像度は、接続されているモニタの最大解像度が 1280×1024 なため、それ以下の解像度に制限されています。ここに、NVIDIA ドライバの機能を使ってカスタム解像度を追加してあげるのが、最初のステップです。</p><p>デスクトップを右クリックして出てくるコンテキストメニューから、「NVIDIA コントロールパネル」を起動します。左側のタスクの選択から、ディスプレイ - 解像度の変更 をクリックします。</p><p>「2. 解像度を選択します。」のところに「カスタマイズ」ボタンがあるので、それを押します。カスタマイズウィンドウが開きます。</p><p>「カスタム解像度の作成」ボタンを押すと、またウィンドウが開きます。ここで作成したい解像度を入力します。私は以下のようにしました。</p><p></p><ul style="text-align: left;"><li>水平ピクセル: 2880</li><li>垂直ライン: 1760 (Chrome リモートデスクトップのヘッダが画面上部に 40 ピクセルほどあるので、この分 1800 から減らしています)</li><li>リフレッシュレート: 75Hz</li></ul><div>「テスト」ボタンを押してみると、2880×1760 に解像度が切り替わるはずです。「この解像度を保存しますか?」と聞かれるので、「はい」と答えます。</div><div>ここで「はい」と答えても、解像度は元に戻ってしまいます。「2. 解像度を選択します。」の解像度の選択肢に 2880×1760 が追加されていますが、ここでそれを選択して「適用」ボタンを押しても、解像度は切り替わりません。それどころか選択肢から消えてしまいます。でも、ここは一旦これで問題ありません。NVIDIA コントロールパネルを終了してください。</div><div><br /></div><h3 style="text-align: left;">ディスプレイのモードを選択する</h3><div>続いて、デスクトップのコンテキストメニューから、「ディスプレイ設定」を起動します。下の方に「ディスプレイの詳細設定」というリンクがあるので、押してください。</div><div>「ディスプレイの情報」というセクションに接続されている物理モニタの情報が表示されているかと思います。その下の方に「ディスプレイ1のアダプターのプロパティを表示します」というリンクがあるので、押します。プロパティウィンドウが開きます。</div><div>「アダプター」タブの下の方に「モードの一覧」というボタンがあるので、押します。</div><div>「有効なモードの一覧」が表示されますが、この中に先程作成したカスタム解像度「2880 x 1760, True Color (32ビット), 75 ヘルツ」があるはずです。それを選択して「OK」を押します。</div><div>プロパティウィンドウの「OK」も押します。「このディスプレイ設定をそのままにしますか?」と聞かれるので、表示に問題が無ければ「変更を維持する」ボタンを押します。画面が乱れてしまったときは、しばらく待っていれば元に戻るはずです。</div><div>これで解像度の変更は完了です。</div><div><br /></div><h3 style="text-align: left;">表示スケールを調整する</h3><div>元の解像度との兼ね合いだと思うのですが、この時点では表示スケールが200%になっています。ディスプレイ設定で125%や150%など、お好みのスケールに変更しましょう。100%はちょっと字が小さくてつらいと思いますが、目のいい方はそれでも良いかと!</div><div><br /></div><div>以上で終わりです。同じような悩みの方の参考になれば幸いです。</div><div><br /></div><div><br /></div><div><br /></div><div><br /></div><div><br /></div><div><br /></div><div><br /></div><p></p>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-68479745102012463312020-08-31T02:03:00.006+09:002020-08-31T15:18:25.762+09:00Java 15 そろそろ来ますね<div style="text-align: left;"><div>OpenJDK 15 は、Final Release Candidate が出ており、2020/9/15 に GA となる予定です。もう間も無くですね。</div><div><br /></div><div>どんな機能が入るのか、ざっと紹介しておきます。</div><span><a name='more'></a></span><h2 style="text-align: left;">Java 15 新機能一覧</h2><ul style="text-align: left;"><li>JEP 339: <a href="https://openjdk.java.net/jeps/339">Edwards-Curve Digital Signature Algorithm (EdDSA)</a></li><ul><li>OpenSSH の ed25519 鍵などでも用いられている エドワーズ曲線デジタル署名アルゴリズム (EdDSA) が導入されるようです</li><li>TLS 1.3 でも EdDSA 証明書が使えるようになっているので、そういった局面でも使えそうです</li></ul><li>JEP 360: <a href="https://openjdk.java.net/jeps/360">Sealed Classes (Preview)</a></li><ul><li>いよいよ sealed クラスが Preview です</li><li>クラスの継承、インターフェースの実装を制限する機能ですが、record との組み合わせが強力です</li></ul><li>JEP 371: <a href="https://openjdk.java.net/jeps/371">Hidden Classes</a></li><ul><li>ライブラリの内部用のクラスを外部からきっちり隠蔽して、予期せず直接利用されることを防ぐ機構です</li></ul><li>JEP 372: <a href="https://openjdk.java.net/jeps/372">Remove the Nashorn JavaScript Engine</a></li><ul><li>JEP 335 (Java 11) から非推奨とされていた Nashorn JavaScript Engine がいよいよ削除されます</li><li>ECMAScript 仕様の変更が早すぎてついていけなくなったかららしく……</li><li>これからは<a href="https://github.com/graalvm/graaljs">GraalJS</a>を使え、とのことですが、アーキテクチャ大きく変わる気がします<br /></li></ul><li>JEP 373: <a href="https://openjdk.java.net/jeps/373">Reimplement the Legacy DatagramSocket API</a></li><ul><li>JDK 1.0/1.1 時代からの実装が残っていた領域だそうで、大幅書き直しとのこと</li></ul><li>JEP 374: <a href="https://openjdk.java.net/jeps/374">Disable and Deprecate Biased Locking</a></li><ul><li>バイアスロック機構を無効化します</li><li>バイアスロックは Hashtable, Vector などの Java 初期の頃の同期化コレクション向けの機構で、イマドキのアプリケーションでは邪魔でしかなく、今回の無効化で性能が向上する可能性があります</li></ul><li>JEP 375: <a href="https://openjdk.java.net/jeps/375">Pattern Matching for instanceof (Second Preview)</a></li><ul><li>JEP 305 (Java 14) で Preview として入っていましたが、変更無しで再プレビューし、もっとフィードバックを集めるとのことです</li></ul><li>JEP 377: <a href="https://openjdk.java.net/jeps/377">ZGC: A Scalable Low-Latency Garbage Collector</a></li><ul><li>ZGC は JEP 333 (Java 11) で Linux(x64) 用に、JEP 364, 365 (Java 14) で macOS/Windows 用に Experimental 版が入っていましたが、ついに Production 版になります</li><li>Linux(Aarch64) 用が出来て、テストも OK、といったところのようです</li></ul><li>JEP 378: <a href="https://openjdk.java.net/jeps/378">Text Blocks</a></li><ul><li>JEP 368 (Java 14) で Second Preview となっていましたが、今回リリースです</li></ul><li>JEP 379: <a href="https://openjdk.java.net/jeps/379">Shenandoah: A Low-Pause-Time Garbage Collector</a></li><ul><li>Shenandoah GC もリリースです</li><li>高スループットなら G1 GC、低レイテンシなら ZGC か Shenandoah GC という新時代に突入ですね</li></ul><li>JEP 381: <a href="https://openjdk.java.net/jeps/381">Remove the Solaris and SPARC Ports</a></li><ul><li>JEP 362 (Java 14) で予告されていましたが、ついに消えます</li><li>Solaris/SPARC, Solaris/x64, Linux/SPARC が消えます</li><li>時代ですね……</li></ul><li>JEP 383: <a href="https://openjdk.java.net/jeps/383">Foreign-Memory Access API (Second Incubator)</a></li><ul><li>JEP 370 (Java 14) で Incubator 版でしたが、今回 Second Incubator とのこと</li></ul><li>JEP 384: <a href="https://openjdk.java.net/jeps/384">Records (Second Preview)</a></li><ul><li>JEP 359 (Java 14) で Preview だった record が Second Preview へ進んでいます</li><li>ローカルレコードなど、機能追加があるようです</li></ul><li>JEP 385: <a href="https://openjdk.java.net/jeps/385">Deprecate RMI Activation for Removal</a></li><ul><li>RMI Activation を廃止に向けて非推奨にするとのこと</li><li>RMI Activation を除く RMI は残るそうです</li><li>RMI は昔使っていましたが、RMI Activation が何のことか分かりません</li><li>きっと何の影響もないのでしょう……</li></ul></ul></div>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-75152087948867809372020-08-29T18:41:00.003+09:002020-08-29T18:45:55.304+09:00R - 入門本編 5章 配列と行列<div>
R における「配列 (array)」は、ベクトルを多次元に拡張したものである。2次元のものは特に「行列 (matrix)」と呼ぶ。逆に言うと、R では3次元以上のものは行列とは呼ばない。
</div>
<a name='more'></a>
<h2>配列</h2>
<div>
配列の実体は、次元ベクトル dim を属性に持つベクトルである。次元ベクトルとは、正整数の数値ベクトルでできている。次元ベクトルの要素数が配列の次元数であり、次元ベクトルの各要素の値は、各次元の要素数 (添字の上限値) を示す。
</div><div><br /></div>
<h3>生成</h3>
<div>
例えば 2×3 の2次元配列を作るには、以下のような方法がある。
</div><div><br /></div>
<h4>配列を生成する関数 array() を用いた方法</h4>
<pre class="prettyprint lang-r">> a <- array(0, dim=c(2, 3))
> a
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
</pre>
<div>
array() 関数を用いると、データベクトルと次元ベクトルから配列が生成できる。データベクトルの要素数が次元ベクトルから計算される総要素数に足りないときは、例によってデータベクトルの要素が先頭から再利用される。データベクトルの方が大きいときは必要なところまでが使われる。分かりやすい例としては、array(0, c(3, 4, 2)) とすれば、全要素が 0 の配列になる (要素数1のベクトルが必要なだけ再利用されているということ)。
</div><div><br /></div>
<h4>2×3 = 6 要素のベクトルに次元ベクトルを付加する方法</h4>
<pre class="prettyprint lang-r" style="text-align: left;">> a <- rep(0, 6)
> dim(a) <- c(2, 3)
> a
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
</pre>
<div><br /></div><h4>dim 属性を dim() ではなく attr() で直接操作する方法</h4>
<div>
出来るというだけで、当然だがおすすめしない。
</div>
<pre class="prettyprint lang-r" style="text-align: left;">> a <- rep(0, 6)
> attr(a, "dim") <- c(2, 3)
> a
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
</pre>
<div style="text-align: left;"><br /></div><h3 style="text-align: left;">添字の順番</h3>
<div>
ベクトルに次元ベクトルを付加して配列とした場合に、添字が割り振られる順番を確認してみる。
</div>
<pre class="prettyprint lang-r">> a <- 1:6
> dim(a) <- c(2, 3)
> a
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
</pre>
<div>
a[1,1] a[2,1] a[1,2] a[2,2] a[1,3] a[2,3] の順に割り当てられていることが分かる。大雑把に言えば、最初の添字が先に変わる順番。これは FORTRAN と同じで、C, Java とは異なる。
</div><div><br /></div>
<h3>添字とサブセクション</h3>
<div>
配列の中身を参照するにあたっては、単純な正整数値の指定だけではなく、添字ベクトルが有効である。次元数分添字ベクトルを「,」区切りで並べる。空の添字ベクトルが与えられた次元は、無条件で全てを対象とする意味になる (TRUE を指定したのと同じ)。
</div>
<pre class="prettyprint lang-r">> a[1,]
[1] 1 3 5
> a[2,]
[1] 2 4 6
> a[,1]
[1] 1 2
> a[,2]
[1] 3 4
> a[,3]
[1] 5 6
> a[,]
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
</pre>
<div>
a[] や a[,] は配列全体の意味になり、つまり a と同じである。
</div>
<div>
添字に単一かつ配列ではないベクトルが指定された場合は、次元ベクトルが無視されデータベクトルに対する参照となる。配列を指定した場合は、以下に述べる。
</div><div><br /></div>
<h3>添字配列</h3>
<div>
添字として、添字ベクトルの代わりに「添字配列」を用いることができる。不規則に配列中の位置を指定して値を取得したり代入したりできる。
</div>
<div>
例えば、4×5 の2次元数値配列 X の中の X[1,3]、X[2,2]、X[3,1] の値を 0 に置き換えたいとする。X は以下のようなもの。
</div>
<pre class="prettyprint lang-r">> X <- array(1:20, dim=c(4, 5))
> X
[,1] [,2] [,3] [,4] [,5]
[1,] 1 5 9 13 17
[2,] 2 6 10 14 18
[3,] 3 7 11 15 19
[4,] 4 8 12 16 20
</pre>
<div>
操作するにあたって、[1, 3], [2, 2], [3, 1] を意味する添字配列 i を作る。
</div>
<pre class="prettyprint lang-r">> i <- array(c(c(1, 2, 3), c(3, 2, 1)), dim=c(3, 2))
> i
[,1] [,2]
[1,] 1 3
[2,] 2 2
[3,] 3 1
</pre>
<div>
X[i] を見てみる。
</div>
<pre class="prettyprint lang-r">> X[i]
[1] 9 6 3
</pre>
<div>
X[i] を 0 で置換する。
</div>
<pre class="prettyprint lang-r">> X[i] <- 0
> X
[,1] [,2] [,3] [,4] [,5]
[1,] 1 5 0 13 17
[2,] 2 0 10 14 18
[3,] 0 7 11 15 19
[4,] 4 8 12 16 20
</pre>
<div>
より複雑な例。blocks, varieties の2つの因子を持つブロック実験計画用の計画行列を作成する。実験はn回の反復がある。blocks は3水準、varieties は2水準、実験は12回とする。
</div>
<pre class="prettyprint lang-r">> blocks <- factor(1:3)
> varieties <- factor(1:2)
> n <- 12
</pre>
<div>
blocks, varieties のそれぞれについて、実験回数×水準数の行列を作り、Xb, Xv とする。
</div>
<pre class="prettyprint lang-r">> Xb <- matrix(0, n, nlevels(blocks))
> Xv <- matrix(0, n, nlevels(varieties))
</pre>
<div>
blocks, varieties のそれぞれについて、実験に水準を順番に割り当てた配列 ib, iv を作る。
</div>
<pre class="prettyprint lang-r">> ib <- cbind(1:n, blocks)
> iv <- cbind(1:n, varieties)
> ib
blocks
[1,] 1 1
[2,] 2 2
[3,] 3 3
[4,] 4 1
[5,] 5 2
[6,] 6 3
[7,] 7 1
[8,] 8 2
[9,] 9 3
[10,] 10 1
[11,] 11 2
[12,] 12 3
> iv
varieties
[1,] 1 1
[2,] 2 2
[3,] 3 1
[4,] 4 2
[5,] 5 1
[6,] 6 2
[7,] 7 1
[8,] 8 2
[9,] 9 1
[10,] 10 2
[11,] 11 1
[12,] 12 2
</pre>
<div>
ib, iv は、1列目に実験の回号、2列目に水準の入った2次元配列である。
</div>
<div>
ib, iv を添字配列として、Xb, Xv の該当箇所に 1 を入れる。
</div>
<pre class="prettyprint lang-r">> Xb[ib] <- 1
> Xv[iv] <- 1
> Xb
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
[4,] 1 0 0
[5,] 0 1 0
[6,] 0 0 1
[7,] 1 0 0
[8,] 0 1 0
[9,] 0 0 1
[10,] 1 0 0
[11,] 0 1 0
[12,] 0 0 1
> Xv
[,1] [,2]
[1,] 1 0
[2,] 0 1
[3,] 1 0
[4,] 0 1
[5,] 1 0
[6,] 0 1
[7,] 1 0
[8,] 0 1
[9,] 1 0
[10,] 0 1
[11,] 1 0
[12,] 0 1
</pre>
<div>
これで、Xb, Xv は各実験で用いる水準をフラグ 1 で示した配列となった。
</div>
<pre class="prettyprint lang-r">> X <- cbind(Xb, Xv)
> X
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 1 0
[2,] 0 1 0 0 1
[3,] 0 0 1 1 0
[4,] 1 0 0 0 1
[5,] 0 1 0 1 0
[6,] 0 0 1 0 1
[7,] 1 0 0 1 0
[8,] 0 1 0 0 1
[9,] 0 0 1 1 0
[10,] 1 0 0 0 1
[11,] 0 1 0 1 0
[12,] 0 0 1 0 1
</pre>
<div>
Xb, Xv を併せれば、求める計画行列になる。因子に名前を付けていないため分かりづらいが、[,1], [,2], [,3] が blocks、[,4], [,5] が varieties に対応している。
</div>
<div>
生起行列は Xb, Xv を用いれば、以下のように求まる。
</div>
<pre class="prettyprint lang-r">> N <- crossprod(Xb, Xv)
> N
[,1] [,2]
[1,] 2 2
[2,] 2 2
[3,] 2 2
</pre>
<div>
table() を用いれば、ib, iv から直接求めることもできる。
</div>
<pre class="prettyprint lang-r" style="text-align: left;">> N <- table(ib[,2], iv[,2])
> N
1 2
1 2 2
2 2 2
3 2 2
</pre>
<div><br /></div><h3>配列の演算</h3>
<div>
配列は数値演算式で用いることができ、データベクトルの要素毎に演算を行った結果の配列となる。基本的には被演算項も同じ次元ベクトルを持つ必要がある。
</div><div><br /></div>
<h4>2つの配列の外積</h4>
<div>
配列に対する重要な演算として「外積 (outer product)」がある。a と b の外積は、次元ベクトルとして a と b の次元ベクトルを連結したものを持ち、データベクトルには a と b の各要素全ての組み合わせの積が格納される。外積は演算子「%o%」または関数 outer() で演算子「*」を用いることで表現する。
</div>
<pre class="prettyprint lang-r">> a <- array(1:2, c(1, 2))
> b <- array(1:12, c(3, 4))
> ab <- a %o% b
> ab <- outer(a, b, "*")
> dim(a)
[1] 1 2
> dim(b)
[1] 3 4
> dim(ab)
[1] 1 2 3 4
</pre>
<div>
例として、各要素が 0~9 の整数値を取る2×2行列 $\left(
\begin{array}{cc}
a & b \\
c & d
\end{array}
\right)$ の行列式 $ad - bc$ の値について、各要素の値が独立に一様に選ばれた場合の行列式の値の分布を求める。
</div>
<pre class="prettyprint lang-r">> d <- outer(0:9, 0:9)
> d
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0 0 0 0 0
[2,] 0 1 2 3 4 5 6 7 8 9
[3,] 0 2 4 6 8 10 12 14 16 18
[4,] 0 3 6 9 12 15 18 21 24 27
[5,] 0 4 8 12 16 20 24 28 32 36
[6,] 0 5 10 15 20 25 30 35 40 45
[7,] 0 6 12 18 24 30 36 42 48 54
[8,] 0 7 14 21 28 35 42 49 56 63
[9,] 0 8 16 24 32 40 48 56 64 72
[10,] 0 9 18 27 36 45 54 63 72 81
</pre>
<div>
ベクトル 0:9 と 0:9 の outer() をとることで、0*0~9*9 の全ての値が含まれる配列ができる。これが $ad$ および $bc$ の値の候補になる。
</div>
<pre class="prettyprint lang-r">> fr <- table(outer(d, d, "-"))
> fr
-81 -80 -79 -78 -77 -76 -75 -74 -73 -72 -71 -70 -69 -68 -67 -66 -65 -64 -63 -62 -61 -60 -59 -58 -57 -56 -55 -54 -53 -52 -51 -50
19 1 2 2 3 2 4 2 4 41 4 4 8 6 6 10 7 27 49 8 8 17 8 12 18 53 13 60 12 18 22 16
-49 -48 -47 -46 -45 -44 -43 -42 -41 -40 -39 -38 -37 -36 -35 -34 -33 -32 -31 -30 -29 -28 -27 -26 -25 -24 -23 -22 -21 -20 -19 -18
35 70 22 24 66 28 18 72 22 75 37 34 26 111 63 36 45 84 34 94 36 93 97 50 53 156 42 60 103 107 50 168
-17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
51 140 112 116 59 191 65 126 156 185 115 206 117 179 153 156 111 570 111 156 153 179 117 206 115 185 156 126 65 191 59 116
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
112 140 51 168 50 107 103 60 42 156 53 50 97 93 36 94 34 84 45 36 63 111 26 34 37 75 22 72 18 28 66 24
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
22 70 35 16 22 18 12 60 13 53 18 12 8 17 8 8 49 27 7 10 6 6 8 4 4 41 4 2 4 2 3 2
79 80 81
2 1 19
</pre>
<div>
d と d の outer() を演算子「-」でとることで、$ad - bc$ としてありうる全ての値がとれる。table() で度数を求めている。上記の例では分かりにくいが、奇数行が出現した値で、偶数行はその上の行の値が出現した回数が表示されている。これを plot() でヒストグラムにするには、以下のようにする。
</div>
<pre class="prettyprint lang-r">> plot(as.numeric(names(fr)), fr, type="h", xlab="Determinant", ylab="Frequency")
</pre>
<div>
<div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-0veOji6YUZw/X0oja1pS5nI/AAAAAAAAAQQ/2seWb5Ea8U0hLpPAq_pOiuGtMbgEAh9HgCLcBGAsYHQ/s915/plot.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="317" data-original-width="915" height="222" src="https://1.bp.blogspot.com/-0veOji6YUZw/X0oja1pS5nI/AAAAAAAAAQQ/2seWb5Ea8U0hLpPAq_pOiuGtMbgEAh9HgCLcBGAsYHQ/w640-h222/plot.png" width="640" /></a></div><br />X 軸の値を as.numeric() で数値に強制変換しているのは、出現した数値が fr 内では names 属性の値として文字列になってしまっているため。仮にループで解こうとすると、以下のようになる。
</div>
<pre class="prettyprint lang-r" style="text-align: left;">> fr <- integer(); for (a in 0:9) for (b in 0:9) for (c in 0:9) for (d in 0:9) { det <- as.character(a*d - b*c); fr[det] <- ifelse(is.na(fr[det]), 0, fr[det]) + 1 }
</pre>
<div><br /></div><h4>配列の一般化転置</h4>
<div>
aperm(a, perm) は置換ベクトル perm を用いて a の次元を入れ替える。perm は 1~dim(a) の置換になっている必要がある。新しい配列の j 番目の次元に perm[j] を持ってくる。行列の転置の一般化である。つまり、行列 A の転置は aperm(A, c(2, 1)) で得られる。
</div>
<pre class="prettyprint lang-r" style="text-align: left;">> A <- matrix(1:6, c(2, 3))
> A
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> aperm(A, c(2, 1))
[,1] [,2]
[1,] 1 2
[2,] 3 4
[3,] 5 6
</pre>
<div><br /></div><h2>行列</h2>
<div>
「行列 (matrix)」とは R においては2次元の配列のことを指す。R ではその他の次元のものは行列とは呼ばない。行列には行列専用の便利な演算、関数が定義されている。
</div><div><br /></div>
<h3>転置、行数、列数</h3>
<div>
t() は行列の転置を行う。任意の次元の配列における一般化転置は aperm() を用いるが、行列においては t() で簡便に行うことができる。nrow(), ncol() はそれぞれ行数、列数を取得する関数である。これらは配列に対しても1次元目、2次元目の要素数を返す関数として機能する。
</div><div><br /></div>
<h3>作成</h3>
<div>
行列は matrix() で作る。
</div>
<pre class="prettyprint lang-r" style="text-align: left;">> A <- matrix(c(c(1, 2), c(3, 4)), 2)
> B <- matrix(c(c(0, -1), c(-1, 0)), 2)
> A
[,1] [,2]
[1,] 1 3
[2,] 2 4
> B
[,1] [,2]
[1,] 0 -1
[2,] -1 0
</pre>
<div><br /></div><h3>積</h3>
<div>
演算子「%*%」は行列の積を求める演算子である。ベクトルも与えることができて適宜行ベクトルまたは列ベクトルとして解釈される。n×1行列、1×n行列も行ベクトルまたは列ベクトルとして扱われる。ベクトルが行ベクトルか列ベクトルかの解釈は、例えば x %*% x といった式では曖昧になる。t(x) %*% x と x %*% t(x) のどちらに変換すべきかは自明ではない。R では、より小さい行列、すなわち t(x) %*% x のスカラー値が結果となる。他方の値を求めるには x %*% t(x) の他に cbind(a) %*% a や a %*% rbind(a) のようにいずれかのオペランドを行列として明示しても良い。もちろん cbind(a) %*% rbind(a) でも良い。
</div>
<pre class="prettyprint lang-r">> A * B
[,1] [,2]
[1,] 0 -3
[2,] -2 0
</pre>
<div>
A * B は単に成分毎の積からなる行列であり、行列の積ではない。
</div>
<pre class="prettyprint lang-r">> A %*% B
[,1] [,2]
[1,] -3 -1
[2,] -4 -2
</pre>
<div>
%*% で行列の積が求まる。
</div>
<pre class="prettyprint lang-r">> A %*% B[,1]
[,1]
[1,] -3
[2,] -4
> A %*% B[,2]
[,1]
[1,] -1
[2,] -2
> A[1,] %*% B
[,1] [,2]
[1,] -3 -1
> A[2,] %*% B
[,1] [,2]
[1,] -4 -2
</pre>
<div>
%*% で行列とベクトルの積も計算できる。
</div>
<div style="text-align: left;"><br /></div><h3>クロス積</h3>
<div>
crossprod() 関数は行列のクロス積を計算する。crossprod(A, B) は t(A) %*% B と等しいが、より効率的に動作する。
</div>
<pre class="prettyprint lang-r" style="text-align: left;">> crossprod(A, B)
[,1] [,2]
[1,] -2 -1
[2,] -4 -3
> t(A) %*% B
[,1] [,2]
[1,] -2 -1
[2,] -4 -3
</pre>
<div><br /></div><h3>その他行列演算</h3>
<div>
diag() 関数は、ベクトルを渡すとそのベクトルの値を対角成分に持つ行列を返す。行列を渡すと対角成分を抜き出したベクトルを返す。Matlab と同様の挙動。正整数 k を引数に渡すと、k×kの単位行列を返す。ややこしい。
</div>
<div>
solve(A, b) 関数は線形方程式 $Ax = b$ を解いて x を求める。
</div>
<pre class="prettyprint lang-r">> A <- matrix(c(3, 2, 4, 5), 2)
> b <- c(11, 12)
> solve(A, b)
[1] 1 2
</pre>
<div>
solve(A) は A の逆行列を求める。しかし、solve(A, b) は solve(A) %*% b よりも効率が良いので、二次形式を求めるときも x %*% solve(A) %*% x とするより x %*% solve(A, x) とした方が効率が良い。solve(A) を用いることは稀だ。
</div>
<div>
eigen(Sm) は対称行列 Sm の固有値と固有ベクトルを求める。固有値 values と固有ベクトル vectors のリストとして返ってくる。eigen(Sm)$values, eigen(Sm)$vectors などとすればそれぞれ固有値、固有ベクトルのみが取り出せる。$ は再帰的構造のオブジェクトから要素を取り出す演算子。
</div>
<div>
svd(M) は行列 M の特異値分解を行う。特異値分解とは、M を U %*% D %*% t(V) の形に分解することである。M が m 行 n 列の行列とすると、U は m 行 m 列のユニタリ行列、V は n 行 n 列のユニタリ行列、D は m 行 n 列で対角成分は非負、それ以外の成分は 0 の行列となる。svd(M) の結果は d, u, v の3つの要素を持つリストとして返される。ただし d は行列ではなく対角成分を並べたベクトルである (行列にしたければ diag(d))。M の行列式の絶対値は prod(svd(M)$d) で求まる。
</div>
<div>
正方行列のトレース (対角和) tr() という関数を定義してみる。
</div>
<pre class="prettyprint lang-r">> tr <- function(X) ifelse(is.matrix(X), sum(diag(X)), X[1])
> tr(A)
[1] 8
> tr(9.9)
[1] 9.9
> tr(-9)
[1] -9
> tr(c(1,2,3,4,4)
+ )
[1] 1
> tr(integer())
[1] NA
</pre>
<div>
行列が指定された場合は、diag() で対角成分を取り出し、sum() で総計する。それ以外が指定されたときはベクトルとみなして先頭の要素を返す。結果としてスカラーなら指定された値がそのまま返る。3次元以上の配列が指定されたときも先頭の要素が返るが、それが合理的かどうかはわからない。
</div>
<div>
関数 lsfit(X, y) は計画行列 X と観測値のベクトル y に対して最小二乗法によるフィッティングを行う (Least squares fit)。ヘルプに載っている例は以下のもの。
</div>
<pre class="prettyprint lang-r">> ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
> trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
> weight <- c(ctl, trt)
> lsD9 <- lsfit(x = unclass(gl(2,10)), y = weight)
> ls.print(lsD9)
Residual Standard Error=0.6964
R-Square=0.0731
F-statistic (df=1, 18)=1.4191
p-value=0.249
Estimate Std.Err t-value Pr(>|t|)
Intercept 5.403 0.4924 10.9723 0.000
X -0.371 0.3114 -1.1913 0.249
</pre>
<div>
lsfit(X, y) は内部で QR 分解を行っている。
</div>
<pre class="prettyprint lang-r">> X <- cbind(1, unclass(gl(2, 10)))
> Xplus <- qr(X)
</pre>
<div>
この例は、計画行列の QR 分解を行い、Xplus に付値している。
</div>
<div>
qr(X) で X の QR 分解を得る。Xplus はクラス qr で list モードのオブジェクト。qr, rank, qraux, pivot という要素を持つ。qr(X)$qr は X と同じ次元の行列で、上三角部分が R になっており、残りの要素と qr(X)$qraux が Q の情報をコンパクトな形式で持っている。qr.X(X) で元の行列が取り出せ、qr.Q(X), qr.R(X) でそれぞれ分解された直交行列と上三角行列が取り出せる。当然、qr.Q(X) %*% qr.R(X) == qr.X(X) == X。
</div>
<pre class="prettyprint lang-r">> b <- qr.coef(Xplus, y)
> fit <- qr.fitted(Xplus, y)
> res <- qr.resid(Xplus, y)
> b
X
5.403 -0.371
> fit
[1] 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 5.032 4.661 4.661 4.661 4.661 4.661 4.661 4.661 4.661 4.661 4.661
> res
[1] -0.862 0.548 0.148 1.078 -0.532 -0.422 0.138 -0.502 0.298 0.108 0.149 -0.491 -0.251 -1.071 1.209 -0.831 1.369
[18] 0.229 -0.341 0.029
</pre>
<div>
b, fit, res はそれぞれ、y の X の範囲への射影の係数ベクトル (coefficients)、直交射影、直交補空間への射影 (residuals)。X はフルランク (全ての行または列が線形独立) でなくても良く、冗長なものは取り除かれる。
</div><div><br /></div>
<h3>cbind() と rbind()</h3>
<div>
ベクトルや行列を合併して新しい行列を作る関数として、cbind(), rbind() がある。与えられたベクトルや行列を順番に cbind() では列方向に (右に)、rbind() では行方向に (下に) 並べる。
</div>
<div>
cbind() の引数に行列を指定する場合は、指定した全ての行列の行数が等しい必要がある。ベクトルの要素数は任意でよく、足りなければ先頭から再利用される。行列が指定されていない場合は要素数が最大のベクトルの要素数を行数とした行列が作られる。短いベクトルの要素の再利用もその要素数まで行われる。ここでの再利用は、ターゲットの行数が行数の少ない行列の行数や要素数の少ないベクトルの要素数の正整数倍でなければならない。そうでない場合はエラーとなる。
</div>
<div>
rbind() も行と列の位置付けが逆なことを除いて cbind() と同様。
</div><div><br /></div>
<h4>配列に c() を適用すると……</h4>
<div>
c() も配列に対して適用することができる。しかし、cbind(), rbind() が dim 属性を意識しているのに対し、c() は無視している。つまり、データベクトルそのままの素のベクトルとみなして処理する。配列 X に対して c(X) とするだけでもデータベクトルが取り出される。
</div><div><br /></div>
<h2>因子から度数分布表を作る</h2>
<div>
table() 関数は因子を用いて度数分布表を作成する関数である。4章の例を用いると、以下のようにできる。
</div>
<pre class="prettyprint lang-r">> table(statef)
statef
act nsw nt qld sa tas vic wa
2 6 2 5 4 2 5 4
</pre>
<div>
この例では、tapply(statef, statef, length) と同じことが行われている。
</div>
<div>
2つの因子を扱うこともできる。
</div>
<pre class="prettyprint lang-r">> table(statef, sex)
sex
statef f m
act 2 0
nsw 2 4
nt 2 0
qld 1 4
sa 1 3
tas 1 1
vic 4 1
wa 2 2
</pre>
<div>
新たに収入の値で区分する因子 incomef を定義する。
</div>
<pre class="prettyprint lang-r">> incomef <- factor(cut(incomes, breaks = 35 + 10 * (0:7)))
</pre>
<div>
cut(x, breaks) は x に与えられたベクトルを breaks に与えられたベクトルの値を境界にして分解する関数。これを使って3つの因子を指定した場合は以下のような表示になる。
</div>
<pre class="prettyprint lang-r">> table(statef, sex, incomef)
, , incomef = (35,45]
sex
statef f m
act 1 0
nsw 0 1
nt 0 0
qld 0 1
sa 0 0
tas 0 0
vic 1 0
wa 0 0
, , incomef = (45,55]
sex
statef f m
act 1 0
nsw 1 0
nt 1 0
qld 0 1
sa 1 1
tas 0 0
vic 1 0
wa 2 1
, , incomef = (55,65]
sex
statef f m
act 0 0
nsw 1 2
nt 1 0
qld 1 2
sa 0 2
tas 1 1
vic 1 1
wa 0 1
, , incomef = (65,75]
sex
statef f m
act 0 0
nsw 0 1
nt 0 0
qld 0 0
sa 0 0
tas 0 0
vic 1 0
wa 0 0
</pre>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-89534550844306393812020-08-29T14:32:00.002+09:002020-08-29T14:32:26.164+09:00R - 入門本編 4章 順序付き因子と順序無し因子<div>
R では質的変数、すなわち順序尺度、名義尺度を表現する方法として「因子 (factor)」が用意されています。モデル構築にあたっては、順序尺度は「順序付き因子 (ordered factor)」、名義尺度は「順序無し因子 (unordered factor)」として正しくデータを表現しておく必要がありますので、よく理解しておきましょう。
</div>
<a name='more'></a>
<h2>因子とは何か?</h2>
<div>
因子 (factor) は、質的変数を R で表現する方法である。質的変数には、順序尺度と名義尺度があるが、R ではそれぞれに対して順序付き因子と順序無し因子が対応する。因子ではない数値ベクトルは、量的変数 (間隔尺度や比例尺度) とみなされる。R では順序無し因子、順序付き因子、因子ではない数値ベクトルは、モデル構築において異なる取り扱い (各尺度として適切な取り扱い) がされるようになっている。
</div><div><br /></div>
<h3>因子の使用例</h3>
<div>
オーストラリアの会計士30人の本拠地所在州と収入のデータを分析するとする。state は州名略称の文字列ベクトル、incomes は収入の数値ベクトルとする。
</div>
<pre class="prettyprint lang-r">> state <- c("tas", "sa", "qld", "nsw", "nsw", "nt", "wa", "wa", "qld", "vic", "nsw", "vic", "qld", "qld", "sa",
+ "tas", "sa", "nt", "wa", "vic", "qld", "nsw", "nsw", "wa", "sa", "act", "nsw", "vic", "vic", "act")
> incomes <- c(60, 49, 40, 61, 64, 60, 59, 54, 62, 69, 70, 42, 56, 61, 61,
+ 61, 58, 51, 48, 65, 49, 49, 41, 48, 52, 46, 59, 46, 58, 43)
> state
[1] "tas" "sa" "qld" "nsw" "nsw" "nt" "wa" "wa" "qld" "vic" "nsw" "vic" "qld" "qld" "sa" "tas" "sa" "nt" "wa" "vic"
[21] "qld" "nsw" "nsw" "wa" "sa" "act" "nsw" "vic" "vic" "act"
> incomes
[1] 60 49 40 61 64 60 59 54 62 69 70 42 56 61 61 61 58 51 48 65 49 49 41 48 52 46 59 46 58 43
</pre>
<div>ちなみに略称と州名の関係は以下の通り。州と準州合わせて8つある。
<table>
<tbody><tr><td>
act</td><td>the Australian Capital Territory</td>
</tr><tr><td>
nsw</td><td>New South Wales</td>
</tr><tr><td>
nt</td><td>the Northern Territory (準州)</td>
</tr><tr><td>
qld</td><td>Queensland</td>
</tr><tr><td>
sa</td><td>South Australia</td>
</tr><tr><td>
tas</td><td>Tasmania</td>
</tr><tr><td>
vic</td><td>Victoria</td>
</tr><tr><td>
wa</td><td>Western Australia</td>
</tr>
</tbody></table>
</div>
<div style="text-align: left;"><br /></div><h4>因子の作成</h4>
<div>
州名は名義尺度なので、state から順序無し因子 statef を作る。
</div>
<pre class="prettyprint lang-r">> statef <- factor(state)
> statef
[1] tas sa qld nsw nsw nt wa wa qld vic nsw vic qld qld sa tas sa nt wa vic qld nsw nsw wa sa act nsw vic vic act
Levels: act nsw nt qld sa tas vic wa
</pre>
<div>statef を表示してみると、要素が "" で括られずに表示された。これは文字列ベクトルではなくなったためである。Levels 欄には「水準」が表示されている。因子の水準とは、その因子の取り得る値だ。水準はアルファベット順にソートされているように見える。
</div><div><br /></div>
<h4>因子の水準の取得</h4>
<pre class="prettyprint lang-r">> levels(statef)
[1] "act" "nsw" "nt" "qld" "sa" "tas" "vic" "wa"
> nlevels(statef)
[1] 8
</pre>
<div>
levels() で水準の一覧が、nlevels() で水準の数が取れる。
</div><div><br /></div>
<h4>因子の性質</h4>
<div>
因子の性質を把握するため、いろいろ表示してみる。
</div>
<pre class="prettyprint lang-r">> summary(statef)
act nsw nt qld sa tas vic wa
2 6 2 5 4 2 5 4
</pre>
<div>
summary() を使うと、各水準の度数が表示された。
</div>
<pre class="prettyprint lang-r">> mode(statef)
[1] "numeric"
</pre>
<div>
モードは numeric になっており、実体は数値ベクトルのようだ。
</div>
<pre class="prettyprint lang-r">> attributes(statef)
$levels
[1] "act" "nsw" "nt" "qld" "sa" "tas" "vic" "wa"
$class
[1] "factor"
</pre>
<div>
属性には levels に水準を持っており、class は factor だ。
</div>
<pre class="prettyprint lang-r">> unclass(statef)
[1] 6 5 4 2 2 3 8 8 4 7 2 7 4 4 5 6 5 3 8 7 4 2 2 8 5 1 2 7 7 1
attr(,"levels")
[1] "act" "nsw" "nt" "qld" "sa" "tas" "vic" "wa"
</pre>
<div>
unclass() してみると、実体は水準の添字を要素として持つ数値ベクトルになっているようだ。
</div>
<div>
つまり、因子は factor クラスのオブジェクトで、水準に辞書順に数値を振り、要素を水準の添字で置き換えた数値ベクトルということのようだ。
</div><div><br /></div>
<h4>因子を使ったデータ分析</h4>
<div>
因子を活用して、収入の分析を行ってみよう。まずは復習がてら、全標本の標本平均を出してみる。
</div>
<pre class="prettyprint lang-r">> mean(incomes)
[1] 54.73333
</pre>
<div>
続いて、州毎の標本平均を出してみる。tapply() という関数を用いる。
</div>
<pre class="prettyprint lang-r">> incmeans <- tapply(incomes, statef, mean)
> incmeans
act nsw nt qld sa tas vic wa
44.50000 57.33333 55.50000 53.60000 55.00000 60.50000 56.00000 52.25000
</pre>
<div>
tapply(X, INDEX, FUN) は、ベクトル X を因子 INDEX を用いて水準別に分解し、それぞれに関数 FUN を適用して、その結果を水準名で名前付けした配列で返すもの。incomes と statef は同じ要素数なので、incomes に因子 statef を適用できる。適用することで、incomes のどの要素が statef のどの水準に対応するかが決定される。それを用いて水準別に incomes の要素を分解し、水準ごとに1つの数値ベクトルを作る (8州あるので、計8個)。その8州ごとに分解されたベクトルをそれぞれ mean() に渡すことで、各州の標本平均が算出される。結果は配列に格納される。
</div>
<div>
州ごとの標本平均に対する標準誤差を計算するとする。標準誤差を計算する組込み関数がないので、独自に定義する (stderr() は標準エラー出力オブジェクトだ)。
</div>
<pre class="prettyprint lang-r">> stderror <- function(x) sqrt(var(x)/length(x))
> incstderror <- tapply(incomes, statef, stderror)
> incstderror
act nsw nt qld sa tas vic wa
1.500000 4.310195 4.500000 4.106093 2.738613 0.500000 5.244044 2.657536
</pre>
<div>
対象がベクトルのときは tapply() を用いるが、配列・行列の場合は apply() を用いる。また、tapply() の第2引数は factor でなくとも as.factor() で強制変換することで因子になりうるものであれば、指定できる。今回の例では state はその条件を満たしているので、tapply(incomes, state, mean) とすることも可能だった。
</div>
<div>
母集団の平均の95%信頼区間を求めてみる。母集団の平均が $\mu$、標本平均が $\mu_0$、標準偏差が $s$、標本数が $n$ のとき、$t = (\mu_0-\mu)/(\frac{s}{\sqrt{n}})$ が自由度 $n-1$ の $t$-分布になることを利用すると、$\mu$ の信頼区間は $\mu_0 \pm t s \sqrt{n}$ として求まる。
</div>
<div>
まずは全標本に対して。
</div>
<div>
95% にあたる t の値を求める。
</div>
<pre class="prettyprint lang-r">> t <- qt(1-(1-0.95)/2, length(incomes)-1)
> t
[1] 2.04523
</pre>
<div>
qt() は $t$-分布の Quantile function。
</div>
<pre class="prettyprint lang-r">> w <- t * sqrt(var(incomes)) / sqrt(length(incomes))
> w
[1] 3.117875
</pre>
<div>
信頼区間の幅 (片側分) を w として計算。
</div>
<pre class="prettyprint lang-r">> r <- c(mean(incomes) - w, mean(incomes) + w)
> r
[1] 51.61546 57.85121
</pre>
<div>
平均 $\pm$ 幅が求める信頼区間。
</div>
<div>
これを州ごとに行ってみる。数値ベクトルを引数に取ってその信頼区間を返す関数を定義しておいて、それを tapply() で適用すれば良い。
</div>
<pre class="prettyprint lang-r">> f <- function(x) {t <- qt(1-(1-0.95)/2, length(x)-1); w <- t * sqrt(var(x)) / sqrt(length(x)); c(mean(x)-w, mean(x)+w)}
> incrange <- tapply(incomes, statef, f)
> incrange
$act
[1] 25.44069 63.55931
$nsw
[1] 46.25363 68.41304
$nt
[1] -1.677921 112.677921
$qld
[1] 42.19966 65.00034
$sa
[1] 46.28451 63.71549
$tas
[1] 54.1469 66.8531
$vic
[1] 41.4402 70.5598
$wa
[1] 43.79253 60.70747
</pre>
<div>
複数のカテゴリを持つ場合にも tapply() は使用できる。性別の因子を付加してみる。
</div>
<pre class="prettyprint lang-r">> sex <- factor(c('m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f'))
> sex
[1] m f m f m f m f m f m f m f m f m f m f m f m f m f m f m f
Levels: f m
</pre>
<div>
tapply() に複数の因子を適用するには、因子を list にまとめれば良い。
</div>
<pre class="prettyprint lang-r">> tapply(incomes, list(statef, sex), mean)
f m
act 44.5 NA
nsw 55.0 58.50
nt 55.5 NA
qld 61.0 51.75
sa 49.0 57.00
tas 61.0 60.00
vic 55.5 58.00
wa 51.0 53.50
</pre>
<div>
この場合、結果は2次元配列で返ってくる。因子の指定の仕方によっては結果が不調和配列になっている可能性がある。
</div>
<div>
※ R 入門には不調和配列 (Ragged array) と書かれている。通常 Ragged array と言えば、列数の異なる行が集まっている配列のことだが、R の配列はそうはなり得ないように思われる。ここで Ragged array と言っているのは、因子の組み合わせに該当する要素が無いために結果に NA が含まれることを言っているのだと思う。
</div><div><br /></div><h2 style="text-align: left;">順序付き因子</h2>
<div>
順序付き因子を作るには factor() のかわりに ordered() を用いる。順序付き因子はクラス ordered でクラス factor も継承している。
</div>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-7876568294364710152020-08-23T03:46:00.003+09:002020-08-23T03:48:55.438+09:00R - 入門本編 3章 オブジェクト、そのモードと属性<div>
R の型システムの入門編。型をいちいち明示しない言語では、型システムをちゃんと理解しておかないとちょいちょいつまずくので、しっかり把握しておきたい。
</div><div><br /></div>
<a name='more'></a>
<h2>本質的属性 - モードと要素数</h2>
<div>
R では、処理の対象は全て「オブジェクト」である。オブジェクトには様々な種類があり、オブジェクトは「モード」と「要素数」という2つの本質的属性を持つ。
</div><div><br /></div>
<h3>ベクトルのモード</h3>
<div>
まずアトミックな構造として、「数値 (numeric)」、「複素数 (complex)」、「論理値 (logical)」、「文字列 (character)」の各ベクトルがある。この numeric, complex, logical, character を「モード」と呼んでいる。ベクトルは全て同じモードの値からなる。空のベクトルもモードを持つ。
</div>
<div><br /></div><div>
mode() は引数に与えられたオブジェクトのモードを返す。
</div>
<pre class="prettyprint lang-r">> mode(1)
[1] "numeric"
> mode(1i)
[1] "complex"
> mode(TRUE)
[1] "logical"
> mode("1")
[1] "character"
</pre>
<div><br /></div><div>
空のベクトルは以下のように作れる。引数の 0 は要素数 0 を意味する。
</div>
<pre class="prettyprint lang-r">> mode(double(0))
[1] "numeric"
> mode(complex(0))
[1] "complex"
> mode(logical(0))
[1] "logical"
> mode(character(0))
[1] "character"
</pre>
<div><br /></div><div>
NA や NaN, Inf もモードを持つ。
</div>
<pre class="prettyprint lang-r">> mode(NA)
[1] "logical"
> mode(NaN)
[1] "numeric"
> mode(Inf)
[1] "numeric"
</pre>
<div><br /></div><div>
しかし、NA は numeric モードのベクトルの要素に記述できる。
</div>
<pre class="prettyprint lang-r" style="text-align: left;">> mode(c(1, 2, NA))
[1] "numeric"
</pre>
<div><br /></div><h3>リスト、関数、式のモード</h3>
<div>
R には、「リスト」と呼ばれる種類のオブジェクトがある。任意のモードの値を複数持つ順序付きの列で、モードは list である。ベクトルがアトミックな構造であるのに対し、リストは再帰的な構造を持つ。簡単に言えば、リストの要素としてリストを持つことができる。
</div>
<div>
他の再帰的な構造としては、「関数 (function)」、「式 (expression)」がある。式がオブジェクトであるのは、R の大きな特徴である。式の一種としてモデル記述の際に用いられる「モデル式 (model formula)」というものがある。
</div><div><br /></div>
<h3>モードの変換</h3>
<div>
モードは変換することができる。例えば numeric モードのベクトルを character モードに変換すると、次のようになる。
</div>
<pre class="prettyprint lang-r">> z <- 0:9
> z
[1] 0 1 2 3 4 5 6 7 8 9
> z <- as.character(z)
> z
[1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9"
</pre>
<div>
as.integer() を用いれば整数型に変換され、モードは numeric に戻る。
</div>
<pre class="prettyprint lang-r">> z <- as.integer(z)
> z
[1] 0 1 2 3 4 5 6 7 8 9
</pre>
<div>
他にも as.xxx() 形式のモード変換関数がいろいろある。
</div><div><br /></div>
<h3>要素数</h3>
<div>
オブジェクトの本質的属性のもう1つが「要素数 (length)」である。length() で取得できる。オブジェクトの要素数は、今の要素数を越える添字値の位置に値を入れることで変化する。
</div>
<div><br /></div><div>
まず、numeric モードの空のベクトルを作る。
</div>
<pre class="prettyprint lang-r">> e <- numeric()
> length(e)
[1] 0
</pre>
<div><br /></div><div>
添字値 3 に値を入れる。
</div>
<pre class="prettyprint lang-r">> e[3] <- 17
> e
[1] NA NA 17
> length(e)
[1] 3
</pre>
<div>
要素数は 3 になり、値を入れていない箇所には NA が入った。ファイルからベクトルやリストに値を読み込む scan() 関数などでも内部的にこのような要素数調整が行われている。
</div>
<div>
要素数を切り詰めるには、添字ベクトルで切り取って再付値すれば良い。
</div><div><br /></div>
<h2>属性</h2>
<div>
オブジェクトには、本質的属性であるモードと要素数以外にも、属性を付加することができる。非本質的属性は、attributes() で取得したり更新したりすることができる。例えば自分で定義した関数には srcref という属性がついていて、関数の定義式を保持している。名前を付けたベクトルには names という属性に名前を保持している。</div><div><br /></div>
<pre class="prettyprint lang-r">> g <- function(x, y) x + y
> attributes(g)
$srcref
function(x, y) x + y
> fruit <- c(5, 10, 1, 20)
> names(fruit) <- c("orange", "banana", "apple", "peach")
> attributes(fruit)
$names
[1] "orange" "banana" "apple" "peach"
</pre>
<div><br /></div><div>
属性の取得・更新は、attributes(x) で x の全属性を名前付きリストとして取得し、これを操作すれば良い。また、attr() という関数で指定した名前の属性を操作することもできる。
</div>
<div>
names() で付けた名前が names という属性として保持されているなど、R の体系の一部をなしているので、付加したり除去したりする際には注意が必要である。例えば、行列オブジェクトは dim という属性で行数、列数を保持しているので、dim を変更することで行数、列数が変化する。
</div>
<pre class="prettyprint lang-r">> z <- matrix(nrow=4, ncol=4)
> z
[,1] [,2] [,3] [,4]
[1,] NA NA NA NA
[2,] NA NA NA NA
[3,] NA NA NA NA
[4,] NA NA NA NA
> attr(z, "dim") <- c(2, 8)
> z
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] NA NA NA NA NA NA NA NA
[2,] NA NA NA NA NA NA NA NA
> attr(z, "dim") <- c(8, 8)
以下にエラー attr(z, "dim") <- c(8, 8) :
dims [product 64] はオブジェクト [16] の長さに整合しません
</pre>
<div>
要素数は本質的属性でありこの方法では変更できないため、要素数が変化するような指定はエラーになっている。
</div><div><br /></div>
<h2>オブジェクトのクラス</h2>
<div>
R ではオブジェクト指向スタイルのプログラミングが可能となっている。class という特別な属性を用いる。例えばあるオブジェクトが data.frame というクラスを持っている場合、print(), plot(), summary() 等の総称的関数が data.frame というクラスを意識した振舞いをする。unclass() 関数によって一時的にクラスの影響を取り外すことができる。
</div>
<div><br /></div><div>
例として、データフレームを作ってみる。
</div>
<pre class="prettyprint lang-r">> df <- data.frame(cbind(x=1, y=1:10))
> df
x y
1 1 1
2 1 2
3 1 3
4 1 4
5 1 5
6 1 6
7 1 7
8 1 8
9 1 9
10 1 10
</pre>
<div>
cbind() で列名 x と y の2列を持つ10行の行列を作り、そこからデータフレームを作っている。表示すると、行列らしく綺麗に表示されている。
</div><div><br /></div>
<pre class="prettyprint lang-r">> class(df)
[1] "data.frame"
</pre>
<div>
class() でクラス名を見ると、data.frame となっている。
</div><div><br /></div>
<pre class="prettyprint lang-r">> attr(df, "class")
[1] "data.frame"
</pre>
<div>
class という属性に格納されていることがわかる。
</div><div><br /></div>
<pre class="prettyprint lang-r">> unclass(df)
$x
[1] 1 1 1 1 1 1 1 1 1 1
$y
[1] 1 2 3 4 5 6 7 8 9 10
attr(,"row.names")
[1] 1 2 3 4 5 6 7 8 9 10
</pre>
<div>
unclass() でクラスの影響を取り外して表示してみると、上記のようにリストとして取り扱われて表示される。データフレームの正体は列ごとのベクトルの名前付きリストに row.names という属性がついたものだと分かる。
</div>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-85583368402922729482020-08-20T01:32:00.003+09:002020-08-20T01:40:32.807+09:00R - 入門本編 2章 数字とベクトル<div style="text-align: left;">R を使うにあたって、ベクトルの理解は必須です。ベクトルを基本的なデータ構造とし、様々な演算子や関数がベクトルに対して定義されることにより、データ分析のためのプログラミングが非常に簡潔に記述できるようになっています。</div><div style="text-align: left;"><a href="https://cran.r-project.org/doc/manuals/r-release/R-intro.pdf">An Introduction to R</a> の2章に沿いつつ (古い版ですが……)、若干の+αがあります。<br /></div><div style="text-align: left;"><br /></div><span><a name='more'></a></span><h2 style="text-align: left;">ベクトルと付値</h2>
<div>
R の最も基礎的なデータ構造は「ベクトル」。順序付きの、値の列。ベクトルを作るには関数 c() を用いる (combine の c)。<br />
<pre class="prettyprint lang-r">> c(3.1, 4.1, 5.9, 2.6, 5.3)
</pre>
とやれば、5つの数値からなるベクトルが作られる。このように式だけを書いて実行すると、結果の値が表示され、捨てられる (実際には次の式が評価されるまでは .Last.value という名前で参照できる)。<br />
単独で現れた数値は、長さ1のベクトルと見なされる。
</div>
<div>
値に名前を付けるには、「付値」を行う。付値には演算子「<-」または「->」を用いる。<br />
<pre class="prettyprint lang-r">> x <- c(3.1, 4.1, 5.9, 2.6, 5.3)
</pre>
または<br />
<pre class="prettyprint lang-r">> c(3.1, 4.1, 5.9, 2.6, 5.3) -> x
</pre>
で、5つの数値からなるベクトルに x という名前が付く。ちなみに付値の式にも結果の値は存在する (付値された値になる) が、付値を行ったときは、結果の値は表示されない。<br />
「<-」も「->」も意味は同じで、左辺と右辺の役割が逆になるだけ。
</div>
<div>
ところで、他の言語では「値 v を変数 x に代入する」と言うところを、R では「値 v を名前 x に付値する」と言っている。関数型言語ではよく「束縛」という言葉も使われる。<br />
「変数への代入」よりも、「名前への付値・束縛」の方が抽象度が高い。「代入」では、同じメモリ空間を占め続けていることが暗に仮定されているが、「付値・束縛」ではそうではない。<br />
R では、「変数」と「名前」の違いや、「代入」と「付値」の違いはあまり意識する必要はなさそうだが、for 構文のカウンタは変数と呼ばれていたりするので、同じメモリ空間を占め続けて値が変化していくものは変数と呼び、そうでないものは単に名前と呼んでいるようだ。<br />
なお、演算子「=」は別の役割があるため、付値には使えない。
</div>
<div>
付値は関数 assign() を用いて行うことも可能。上記の5要素のベクトルの付値と同じことを assign() を用いて行うと以下のようになる。<br />
<pre class="prettyprint lang-r">> assign("x", c(3.1, 4.1, 5.9, 2.6, 5.3))
</pre>
assign() を用いると、構文解釈上、識別子になれない名前に付値することも可能。<br />
<pre class="prettyprint lang-r">> assign("a b", 7)
</pre>
値を参照するには get() を用いる。<br />
<pre class="prettyprint lang-r">> get("a b")
[1] 7
</pre>
`` で括っても識別子とみなされるようだ。<br />
<pre class="prettyprint lang-r">> `a b`
[1] 7
</pre>
</div>
<div>
上記のように x に付値してある状態で、以下のような式を評価してみる。<br />
<pre class="prettyprint lang-r">> 1/x
[1] 0.3225806 0.2439024 0.1694915 0.3846154 0.1886792
</pre>
ベクトルの各要素の逆数が計算され、付値ではない式の評価なので、結果の値が表示された。結果もまた5要素のベクトルである。ベクトル演算については次節で詳しく見る。
</div>
<div>
c() はベクトルの連結も行える。<br />
<pre class="prettyprint lang-r">> y <- c(x, 0, x)
> y
[1] 3.1 4.1 5.9 2.6 5.3 0.0 3.1 4.1 5.9 2.6 5.3
</pre>
</div>
<div>
c() の機能はもっと奥が深そうだが、深堀りは後の章にて。<br />
</div><div><br /></div>
<h2>ベクトル演算</h2>
<div>
算術式内にベクトルを記述すると、その演算はベクトルの各要素に対して行われる。ベクトルの要素数は異なっていてもよく、最も要素数の多いベクトルと同じ要素数のベクトルが返り値となる。要素数の少ないベクトルは足りない分必要なだけ要素が再利用される。<br />
例えば、<br />
<pre class="prettyprint lang-r">> x <- c(3.1, 4.1, 5.9, 2.6, 5.3)
> y <- c(x, 0, x)
> v <- 2*x + y + 1
</pre>
とすると、一番要素数の多い y の 11 にあわせて、x は2回と1要素繰り返され、1 (これは要素数 1 のベクトル) は 11 回繰り返される。<br />
と、あるのだが、実際やってみると、<br />
<pre class="prettyprint lang-r"> 警告メッセージ:
In 2 * x + y :
長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません
</pre>
と言われてしまう。確かに中途半端に再利用されてほしいことって実用上無いだろうから、倍数になっていないときは何らかミスってるよね。<br />
しかし新しい An Introduction to R でも同じ説明なんだよなぁ。改訂しないのかな。<br />
y の要素数を 10 にしてやりなおしてみる。<br />
<pre class="prettyprint lang-r">> y <- c(x, x)
> v <- 2*x + y + 1
> v
[1] 10.3 13.3 18.7 8.8 16.9 10.3 13.3 18.7 8.8 16.9
</pre>
無事計算できた。
</div>
<div>
+, -, *, /, ^(べき乗), log(), exp(), sin(), cos(), tan(), sqrt() なんかはみな同様にベクトルを受けとり、同じ要素数のベクトルを返す。
</div>
<div>
その他ベクトル演算の代表的なものとしては、max(), min() (それぞれ最大値、最小値)、range(x) (c(min(x), max(x))、length() (要素数)、sum()、prod() (全要素を掛けた値) などがある。<br />
統計関数として、mean(x) (sum(x)/length(x) つまり標本平均)、var(x) (sum((x-mean(x))^2)/(length(x)-1) つまり標本分散) も用意されている。var() は n × p 行列を引数にとると、それを n 個の独立な p-変量標本ベクトルとみなし、p × p 標本共分散行列を返す。<br />
sort() は昇順ソート。sort.list() はその拡張。order() は値の小さい順に1から順番に番号付けを行ったベクトルを返す。
</div>
<div>
max(), min() は引数に複数のベクトルを与えた場合、全体を通して最大、最小の値を返すが、pmax(), pmin() は引数に与えられた各ベクトルの要素位置ごとに最大値、最小値を求め、それを並べたベクトルを返す。<br />
<pre class="prettyprint lang-r">> pmax(c(0,1,0,1),c(1,0,1,0))
[1] 1 1 1 1
</pre>
足りない要素は例によって補完され、結果ベクトルの要素数は引数のベクトルのうち最大の要素数のものと同じ。<br />
<pre class="prettyprint lang-r">> pmax(c(2,0),c(1,0,1,0))
[1] 2 0 2 0
</pre>
</div>
<div>
数値には、整数、実数、複素数がある。数値リテラルは実数型 (double) になる。数値の後ろに i を付けると複素数の虚部として扱われ、虚部を含む数値は複素数型 (complex) になる。<br />
整数型 (integer) は integer() (0 で初期化された整数型のベクトルを作る), as.integer() (引数の値を整数型に変換したものを返す) などで明示することによって利用できる。<br />
内部的な計算は倍精度浮動小数点数または倍精度複素数で行われている。<br />
複素数を取り扱いたい場合は、虚部を明示的に与えなければならない。例えば、複素数に拡張された sqrt() を使いたい場合は以下のようになる。
</div>
<div>
実数を引数にした場合:<br />
<pre class="prettyprint lang-r">> sqrt(-17)
[1] NaN
警告メッセージ:
In sqrt(-17) : 計算結果が NaN になりました
</pre>
</div>
<div>
複素数を引数にした場合:<br />
<pre class="prettyprint lang-r" style="text-align: left;">> sqrt(-17+0i)
[1] 0+4.123106i
</pre>
</div>
<div><br /></div><h2>規則的な数列の生成</h2>
<div>
演算子「:」で数列を生成できる。1:30 は c(1, 2, 3, ..., 28, 29, 30) と同じ。<br />
「:」は結合優先順位が最も高いので、2*1:15 は 2*(1:15) と評価され、c(2, 4, 6, ..., 26, 28, 30) を生成する。おそらく数列リテラルのように扱えることを想定して最強優先順位にしたのだろう。<br />
<pre class="prettyprint lang-r">> n <- 10
> 1:n-1
[1] 0 1 2 3 4 5 6 7 8 9
> 1:(n-1)
[1] 1 2 3 4 5 6 7 8 9
</pre>
30:1 とすれば降順の数列も生成できる (c(30, 29, 28, ..., 3, 2, 1))。
</div>
<div>
seq() はより一般的な数列生成関数。seq() は double 型で、seq.int() という integer 型のバージョンもある。全部で5つの引数を持ち、指定の仕方によって様々な振舞いがある。
</div>
<div>
seq(to) は 1:to と同じ。<br />
<pre class="prettyprint lang-r">> seq(10)
[1] 1 2 3 4 5 6 7 8 9 10
</pre>
seq(from, to) は from:to と同じ。<br />
<pre class="prettyprint lang-r">> seq(3, 12)
[1] 3 4 5 6 7 8 9 10 11 12
</pre>
seq(from, to, by) は、from から to まで増分 by の数列を生成する。<br />
<pre class="prettyprint lang-r">> seq(3, 12, 2)
[1] 3 5 7 9 11
</pre>
seq(length.out=length) は 1 から増分 1 で length 個の要素の数列を生成する。<br />
<pre class="prettyprint lang-r">> seq(length.out=10)
[1] 1 2 3 4 5 6 7 8 9 10
</pre>
seq(from, to, length.out=length) は、from から to までを length 個の要素で等間隔に変化させた数列を生成する。<br />
<pre class="prettyprint lang-r">> seq(2, 20, length.out=10)
[1] 2 4 6 8 10 12 14 16 18 20
</pre>
seq(along.with=x) は、1 から増分 1 で length(x) 個の要素の数列を生成する。<br />
<pre class="prettyprint lang-r">> seq(along.with=c(3.1, 4.1, 5.9, 2.6, 5.3))
[1] 1 2 3 4 5
</pre>
</div>
<div>
ここで出てきた length.out, along.with は、タグと呼ばれるもので、呼び出し時にパラメータの名前を付けることができるもの。他の言語では名前付き引数などと呼ばれているもの。関数のオーバーロード機構がある言語では、引数の型、個数が異なるものは別の実装を持つことができるが、名前付き引数が無い場合、同じ型、同じ個数の引数で別の振舞いをさせたい場合は、違う関数名を付けるしかない (そうしろという話もあるが)。名前付き引数を使えば、同じ関数に相乗りできる (そのかわり当然その関数の仕様はややこしくなる)。<br />
タグ名はタグが特定できる限り短縮して指定することもできる。seq(length.out=length) の場合、seq(length=length)、seq(len=length)、seq(l=length) まで縮めてもタグが特定できるため、正しく呼び出せる。<br />
R での関数呼び出し時の引数の照合は以下のように行われている。<br /><ol style="text-align: left;"><li>
タグに関する正確な照合</li><ul><li>
各名前付き引数 (呼び出し側) に対して、名前が正確に一致する形式的引数 (宣言側) が探索される</li><li>
複数見つかった場合はエラー (同じタグが複数回使われている)</li></ul><li>
タグに関する部分的な照合</li><ul><li>
残りの各名前付き引数が残りの形式的引数に対し部分的照合により比較される</li><li>
もし名前付き引数が、正確にある単独の形式的引数の最初の部分に一致すれば、両者は一致すると見なされる</li><li>
複数の名前付き引数に部分一致した場合はエラーになる</li><ul><li>例えば、fumble, fooey という2つの名前付き形式的引数を持つ関数 f() があったとする</li><li>f(f = 1, fo = 2) は2つ目の引数が fooey のみに一致するものの、1つ目の実引数は fumble, fooey 双方に部分一致するため、エラーとなる (fo が fooey のみに合致するので、f は fumble にあたるといった再帰的な照合はなされない)</li><li>f(f = 1, fooey =2) は2つ目の引数が正確に一致し、部分的照合の対象から外されるため、正しい構文になる</li></ul><li>形式的引数が ‘...’ を含む場合は、部分的照合はそれに先立つ引数に対してのみ適用される</li></ul><li>位置に関する照合</li><ul><li>未照合の全ての形式的引数は、宣言の順に名前無しの引数と結びつけられる</li><li>‘...’ を含む場合は、残りの引数全てが名前の有無にかかわらず引き渡される</li></ul></ol>
seq() の場合は、from, to, by, length.out, along.with というタグ名になっているので、以下のような呼び出しもできる。<br />
<pre class="prettyprint lang-r">> seq(to=10, from=1)
[1] 1 2 3 4 5 6 7 8 9 10
</pre>
</div>
<div>
関数 rep() は様々な方法でベクトルを複製する関数である。<br />
rep(x) は x そのまま。<br />
<pre class="prettyprint lang-r">> rep(x)
[1] 3.1 4.1 5.9 2.6 5.3
</pre>
rep(x, times) は、x を times 回反復したものを返す。<br />
<pre class="prettyprint lang-r">> rep(x, times=2)
[1] 3.1 4.1 5.9 2.6 5.3 3.1 4.1 5.9 2.6 5.3
</pre>
rep(x, times, length.out) は、rep(x, times) の length.out 個目まで。<br />
<pre class="prettyprint lang-r">> rep(x, times=2, length.out=7)
[1] 3.1 4.1 5.9 2.6 5.3 3.1 4.1
</pre>
rep(x, each) は、x の各要素を each 個ずつ反復。<br />
<pre class="prettyprint lang-r">> rep(x, each=3)
[1] 3.1 3.1 3.1 4.1 4.1 4.1 5.9 5.9 5.9 2.6 2.6 2.6 5.3 5.3 5.3
</pre>
rep(x, times, each) は rep(x, each) を times 回反復。<br />
<pre class="prettyprint lang-r" style="text-align: left;">> rep(x, times=2, each=3)
[1] 3.1 3.1 3.1 4.1 4.1 4.1 5.9 5.9 5.9 2.6 2.6 2.6 5.3 5.3 5.3 3.1 3.1 3.1 4.1 4.1 4.1 5.9 5.9 5.9 2.6 2.6 2.6 5.3 5.3 5.3
</pre>
</div>
<div><br /></div><h2>論理ベクトル</h2>
<div>
R では論理ベクトルという形で論理量を取り扱うことができる。論理ベクトルの要素は TRUE, FALSE, NA (欠損値) の3種類の値をとる (クリーネの三値論理相当)。<br />
TRUE, FALSE のかわりにそれぞれ T, F と略記することが可能だが、T, F は既定の変数であって予約語ではない。上書き可能なので、なるべく使用しない方が良い。<br />
<, <=, >, >=, ==, != といった条件演算子による演算結果が論理ベクトルとなる。また、論理演算子として &, |, ! もある。C 言語などと同じ。<br />
論理ベクトルを算術演算の文脈で使用した場合は FALSE は 0、TRUE は 1 に強制変換される。<br />
</div><div><br /></div>
<h2>欠損値</h2>
<div>
論理値が「完全には知られていない」とか、「利用不能」とか、「欠損値である」とか言うとき、NA という識別子を用いることができる。一般的な処理は NA に対して NA を返す (分からないものは分からない、的な)。<br />
そうではない例としては、is.na(x) は x と同じ要素数のベクトルを返し、対応する x の要素が NA のところには TRUE が、そうでないところには FALSE が入る。<br />
<pre class="prettyprint lang-r">> is.na(c(1,2,3,NA))
[1] FALSE FALSE FALSE TRUE
</pre>
NA は厳密には「値」ではないので、NA == NA も TRUE にはならず、NA である。だから is.na(x) と同じことをやろうとして x == NA としてもうまくいかない。<br />
<pre class="prettyprint lang-r">> c(1,2,3,NA) == NA
[1] NA NA NA NA
</pre>
</div>
<div>
数値計算における欠損値として「非数 (Not a Number)」というものもある。いわゆる NaN である。<br />
0/0 や Inf - Inf は NaN を返す (エラーにはならない)。<br />
NaN に対する演算は基本的に全て NaN となる。is.nan() で NaN かどうかのチェックが可能。
</div>
<div>
is.na() は NA に対しても NaN に対しても TRUE となるが、is.nan() は NaN に対してのみ TRUE となる。<br />
</div><div><br /></div>
<h2>文字列ベクトル</h2>
<div>
文字列リテラルは "" (ダブルクオート) で括って表現する。'' (シングルクオート) も可。C 言語スタイルのエスケープシーケンスが使用可能。<br />
c() を使えば複数要素の文字ベクトルにできるのは数値や論理値と同様。<br />
paste() は、任意の数のベクトルを引数にとり、1つの文字列に連結する関数。数値も自明な形で (1 なら "1" に、等) 変換されて連結される。<br />
<pre class="prettyprint lang-r">> paste(c("A", "B"), 1:5)
[1] "A 1" "B 2" "A 3" "B 4" "A 5"
</pre>
戻り値は一番長いベクトルと同じ要素数になっている。短い引数は要素が再利用されている。ここでは要素数が倍数になっている必要はないようだ。<br />
デフォルトの区切り文字は半角スペースだが、名前付き引数 sep を用いると区切り文字を変えることができる。<br />
<pre class="prettyprint lang-r">> paste(c("A", "B"), 1:5, sep="-")
[1] "A-1" "B-2" "A-3" "B-4" "A-5"
</pre>
collapse パラメータを指定すると、結果のベクトルをさらに指定した文字を挟んで連結する。<br />
<pre class="prettyprint lang-r">> paste(c("A", "B"), 1:5, sep="-", collapse="|")
[1] "A-1|B-2|A-3|B-4|A-5"
</pre>
</div>
<div>
ちょっと興味があったので、A, B と 1, 2, 3, 4, 5 の全組み合わせができるか試してみた。<br />
<pre class="prettyprint lang-r">> paste(c(outer(c("A", "B"), 1:5, function(x, y) paste(x, y, sep="-"))), collapse="|")
[1] "A-1|B-1|A-2|B-2|A-3|B-3|A-4|B-4|A-5|B-5"
</pre>
outer() とか匿名関数とか先取りしちゃってるのはすまない。<br />
</div><div><br /></div>
<h2>添字ベクトル</h2>
<div>
ベクトルの後ろに「[]」で括られた「添字ベクトル」をつけることで、ベクトルの一部分を切り出すことができる。<br />
添字ベクトルには4パターンある。<br /><ol style="text-align: left;"><li>
論理ベクトル</li><ul><li>
対象のベクトルと同じ要素数である必要がある</li><li>
添字ベクトル中で TRUE になっている箇所に対応する要素が選び出される</li></ul><li>
正の整数値ベクトル</li><ul><li>
対象のベクトルを x としたときに、1:length(x) の部分集合が指定可能</li><li>
対応する要素をその順番で選び出す</li><li>
添字ベクトルの要素数は任意で、結果の要素数は添字ベクトルの要素数に等しい</li></ul><li>
負の整数値ベクトル</li><ul><li>
正の整数値ベクトルでは選び出す要素を指定したが、負値にすると除外する要素を指定した意味になる</li></ul><li>
文字列ベクトル</li><ul><li>
名前付きベクトル (後述) から要素を選び出すのに用いる</li><li>
指定された文字列に合致する名前を持つ要素が選び出される</li><li>
あとは正の整数値ベクトルと同様</li></ul></ol>
以下、使用例。<br />
<pre class="prettyprint lang-r">> x <- 101:110
> x
[1] 101 102 103 104 105 106 107 108 109 110
</pre>
x は 101~110 の数値ベクトル。
</div>
<div>
<pre class="prettyprint lang-r">> x %% 2 == 0
[1] FALSE TRUE FALSE TRUE FALSE TRUE FALSE TRUE FALSE TRUE
</pre>
偶数のとき TRUE になる論理式。%% は剰余演算。<br />
これを添字ベクトルとして使うと、次のように偶数の要素のみ選び出すことができる。<br />
<pre class="prettyprint lang-r">> x[x %% 2 == 0]
[1] 102 104 106 108 110
</pre>
</div>
<div>
正の整数値ベクトルの場合。<br />
<pre class="prettyprint lang-r">> x[1]
[1] 101
</pre>
これは C 言語などの配列の添字みたいな使い勝手。ただし 1 から始まる。
</div>
<div>
<pre class="prettyprint lang-r">> x[2:4]
[1] 102 103 104
> x[c(1,3,5,7,9)]
[1] 101 103 105 107 109
</pre>
複数要素のベクトルで指定した場合はこんな感じ。
</div>
<div>
次に負の整数値ベクトルの場合。<br />
<pre class="prettyprint lang-r">> x[-c(1,3,5,7,9)]
[1] 102 104 106 108 110
</pre>
1, 3, 5, 7, 9 番目の要素が除外されている。
</div>
<div>
文字列ベクトルの場合。<br />
まずはベクトルに名前を付加する。<br />
<pre class="prettyprint lang-r">> names(x) <- c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J")
> x
A B C D E F G H I J
101 102 103 104 105 106 107 108 109 110
</pre>
names(x) に付値することで x に名前を付加できる。名前のベクトルは x より短かくても良い (足りない要素には名前が付かない) が、多いとエラーとなる。<br />
名前を付けると、連想配列のような使い勝手になる。<br />
<pre class="prettyprint lang-r">> x["A"]
A
101
</pre>
複数要素の選び出しも可。<br />
<pre class="prettyprint lang-r">> x[c("D", "H")]
D H
104 108
</pre>
</div>
<div>
添字ベクトルは、付値される側のベクトルにも付けることができ、該当する要素を入れ替えるような挙動をする。<br />
<pre class="prettyprint lang-r">> x <- -5:5 / -5:5
> x
[1] 1 1 1 1 1 NaN 1 1 1 1 1
> x[is.na(x)] <- 0
> x
[1] 1 1 1 1 1 0 1 1 1 1 1
</pre>
</div>
<div>
<pre class="prettyprint lang-r" style="text-align: left;">> y <- -5:5
> y
[1] -5 -4 -3 -2 -1 0 1 2 3 4 5
> y[y < 0] <- -y[y < 0]
> y
[1] 5 4 3 2 1 0 1 2 3 4 5
</pre>
</div>
<div><br /></div><h2>他の型のオブジェクト</h2>
<div>
ベクトル以外にも以下のような型がある。
</div>
<div><ul style="text-align: left;"><li>
行列、配列 (ベクトルの多次元化)</li><li>
因子 (カテゴリ化)</li><li>
リスト (ベクトルの一般化で、異なる型の要素が混在できたり、ベクトルやリストを要素に持てたりする)</li><li>
データフレーム (行列に似た構造だが、異なる型の列が混在でき、実験データの保持に向いたデータ構造になっている)</li><li>
関数 (関数自体もオブジェクトである)
</li></ul></div>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-10874091289388781772020-08-16T14:27:00.005+09:002020-08-16T18:29:27.545+09:00jnethack-3.6.6-0.1 を macOS Catalina 10.15.6 で動かしてみた<p><a href="https://jnethack.osdn.jp/">JNetHack</a> をやりたくて、ビルドしてみたのですが、ちょっと手こずったので、記録しておきます。参考になれば幸いです。</p><h2 style="text-align: left;">JNetHack って何?</h2><p>知らない方も多いと思うので、簡単に説明しておきますと、JNetHack は <a href="https://www.nethack.org/">NetHack</a> を日本語化したものです。</p><p>NetHack はいわゆる<a href="https://ja.wikipedia.org/wiki/%E3%83%AD%E3%83%BC%E3%82%B0%E3%83%A9%E3%82%A4%E3%82%AF%E3%82%B2%E3%83%BC%E3%83%A0">ローグライクゲーム</a>です。rogue → Hack → NetHack と発展してきた、rogue の直系の子孫的存在です。</p><span><a name='more'></a></span><h2 style="text-align: left;">手順の背景と概要</h2><h3 style="text-align: left;">背景</h3><p>JNetHack は、NetHack のソースコードに対するパッチと、Windows 版のバイナリの2形態で配布されています。macOS 版のバイナリは配布されていません。NetHack は macOS 上でビルドして実行することを考慮しているので、JNetHack も基本的にはビルド・実行が可能です。ただし、いくつかの注意点があります。</p><p></p><ul style="text-align: left;"><li>JNetHack は内部的に日本語が CP932 (Microsoft 系の Shift JIS) エンコーディングになっていることを前提としている</li><li>最近の gcc はソースコードが CP932 で書かれていることに対応していないため、ソースコードは UTF-8 にしておく必要がある</li><li>JNetHack のパッチは CP932 で記述されている</li><li>JNetHack のパッチは macOS 10.10 上でのビルドしか考慮していないが、最新の NetHack は macOS 10.14 上のビルドを考慮できているので、そちらでビルドした方が良さそう</li></ul><p></p><p>JNetHack のビルドは、バージョンが上がると同じ手順では出来なくなることがあるので、こういった背景をおさえておくと、対処しやすいかと思います。</p><p><br /></p><h3 style="text-align: left;">概要</h3><p>おおまかな手順は以下の通りです。</p><p></p><ol style="text-align: left;"><li>brew で必要なツールをインストールする</li><li>NetHack のソースコード、JNetHack のパッチを入手する</li><li>NetHack のソースコードを展開する</li><li>JNetHack のパッチを適用する</li><li>JNetHack 用の変更点を macOS 10.14 用のビルドスクリプトに適用する</li><li>JNetHack 用のビルド前準備スクリプトを実行する</li><li>ビルド実行</li><li>インストール</li></ol><div><br /></div><h2 style="text-align: left;">詳細手順</h2><h3 style="text-align: left;">brew で必要なツールをインストールする</h3><p></p><div>そもそもですが、brew をインストールしていない方はそこからです。</div><div><a href="https://brew.sh/index_ja">https://brew.sh/index_ja</a></div><div><br /></div><div>GCC 10 と libiconv が必要なので、インストールしましょう。</div><div><br /></div>
<pre class="prettyprint lang-sh">brew install gcc@10
brew install libiconv
brew install nkf
</pre>
<div><br /></div><h3 style="text-align: left;">NetHack のソースコード、JNetHack のパッチを入手する</h3><div>NetHack のソースコードは<a href="https://www.nethack.org/v366/download-src.html">ここ</a>から入手してください。nethack-366-src.tgz です。</div><div>JNetHack のパッチは<a href="https://osdn.net/projects/jnethack/releases/72780">ここ</a>の jnethack-3.6.6-0.1.diff.gz です。</div><div>以下、~/Downloads にダウンロードしてあるものとして説明を続けます。</div><div><br /></div><div><br /></div><h3 style="text-align: left;">NetHack のソースコードを展開する</h3><div>まず、作業用のディレクトリを作って、そこに入りましょう。</div><div><br /></div>
<pre class="prettyprint lang-sh">mkdir -p ~/tmp/jnethack
cd ~/tmp/jnethack
</pre>
<div><br /></div><div>そして展開します。</div><div><br /></div>
<pre class="prettyprint lang-sh">tar xvzf ~/Downloads/nethack-366-src.tgz
</pre>
<div><br /></div><div>NetHack-NetHack-3.6.6_Released というディレクトリが作られて、中にいろいろ展開されるはずです。</div><div><br /></div><div><br /></div><h3 style="text-align: left;">JNetHack のパッチを適用する</h3><div>パッチをあてます。</div><div><br /></div>
<pre class="prettyprint lang-sh">cd NetHack-NetHack-3.6.6_Released
~/Downloads/jnethack-3.6.6-0.1.diff.gz |patch -p1
</pre>
<div><br /></div><div><br /></div><h3 style="text-align: left;">JNetHack 用の変更点を macOS 10.14 用のビルドスクリプトに適用する</h3><div>macOS 10.14 用のビルドスクリプトを JNetHack 用に修正します。</div><div><br /></div><pre class="prettyprint lang-sh">vi sys/unix/hints/macosx10.14</pre><div><br /></div><div>修正箇所は以下の diff を見てください。</div><div><br /></div>
<pre class="prettyprint lang-xml">61c61,62
< CC=gcc
---
> CC=gcc-10
> CFLAGS+=-fexec-charset=cp932 -Wno-comment -Wno-unused-parameter -Wno-overlength-strings -Wno-format-overflow
96c97
< WINTTYLIB=-lncurses
---
> WINTTYLIB=-lncurses -liconv</pre><br />
<div><br /></div><div>日本語でも説明しておきます。</div><div><ul style="text-align: left;"><li>61行目の CC を gcc-10 に変更</li><li>61行目の次の行に CLAGS+=... を追記</li><li>96行目 (1行増やしたことを考慮すると97行目) の WINTTYLIB に -liconv を追記</li></ul><div><br /></div></div><div>もう1つ修正するファイルがあります。</div><div><br /></div><pre class="prettyprint lang-sh">vi japanese/set_mac.sh</pre><div><br /></div><div>diff は以下の通り。</div>
<pre class="prettyprint">1c1
< brew link --overwrite gcc@9
---
> brew link --overwrite gcc@10
3c3
< sh sys/unix/setup.sh sys/unix/hints/macosx10.10
---
> sh sys/unix/setup.sh sys/unix/hints/macosx10.14</pre><br />
<br />
<div>要するに GCC 10 を使い、Makefile の元は macosx10.14 を使うように修正します。</div><div><br /></div><div><br /></div><h3 style="text-align: left;">JNetHack 用のビルド前準備スクリプトを実行する</h3>
<pre class="prettyprint lang-sh">sh configure</pre>
<div><br /></div><div><br /></div><h3 style="text-align: left;">ビルド実行</h3><pre class="prettyprint lang-sh">make all</pre><div><br /></div><div><br /></div><h3 style="text-align: left;">インストール</h3><pre class="prettyprint lang-sh">make install</pre><div><br /></div><div>~/bin/jnethack に実行ファイルがインストールされます。データ類は ~/nethackdir に置かれます。</div><div><br /></div><div><br /></div><h2 style="text-align: left;">起動してみる</h2><pre class="prettyprint lang-sh">~/bin/jnethack</pre><div><br /></div><div>以下のように進めることができれば、まずは成功かと思います。</div><div><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-CVWaofOjk2I/XzjBsYnENbI/AAAAAAAAAPY/8UGKtqJiTk8Xl2B1ik_HAf3YX9CuRdajQCLcBGAsYHQ/s937/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-08-16%2B14.17.49.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="664" data-original-width="937" src="https://1.bp.blogspot.com/-CVWaofOjk2I/XzjBsYnENbI/AAAAAAAAAPY/8UGKtqJiTk8Xl2B1ik_HAf3YX9CuRdajQCLcBGAsYHQ/s640/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-08-16%2B14.17.49.png" width="640" /></a></div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-4YLHAA2V1ng/XzjBsRNh3pI/AAAAAAAAAPc/ool0OLDl4vkIMLzEANFD9lCmwfEJBpXyQCLcBGAsYHQ/s937/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-08-16%2B14.17.55.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="664" data-original-width="937" src="https://1.bp.blogspot.com/-4YLHAA2V1ng/XzjBsRNh3pI/AAAAAAAAAPc/ool0OLDl4vkIMLzEANFD9lCmwfEJBpXyQCLcBGAsYHQ/s640/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-08-16%2B14.17.55.png" width="640" /></a></div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-QNtxXqYsa7E/XzjBtIQTYRI/AAAAAAAAAPg/Ixo7fjGEFOcJTw8WhppJ-KLTBCsicjEzQCLcBGAsYHQ/s937/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-08-16%2B14.18.05.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="664" data-original-width="937" src="https://1.bp.blogspot.com/-QNtxXqYsa7E/XzjBtIQTYRI/AAAAAAAAAPg/Ixo7fjGEFOcJTw8WhppJ-KLTBCsicjEzQCLcBGAsYHQ/s640/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-08-16%2B14.18.05.png" width="640" /></a></div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-SWmVBmk_BOA/XzjBtHjQ7aI/AAAAAAAAAPk/WP0sXFYnvGQxPlOksUrObysCpVwbEjBqQCLcBGAsYHQ/s937/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-08-16%2B14.18.08.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="664" data-original-width="937" src="https://1.bp.blogspot.com/-SWmVBmk_BOA/XzjBtHjQ7aI/AAAAAAAAAPk/WP0sXFYnvGQxPlOksUrObysCpVwbEjBqQCLcBGAsYHQ/s640/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-08-16%2B14.18.08.png" width="640" /></a></div><div><br /></div><div><br /></div>
管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-4565417534196637012020-07-28T01:16:00.001+09:002020-07-28T01:17:48.013+09:00第1期決算出しました決算処理を終えました。確定申告も済ませました。<div>今年は難しいお金の動きはなかったので、何とか自力で終えられました。よかったです。</div><div>第2期は税理士さん頼まないとちょっと厳しいかもしれません。いっぱい稼がないとですね。</div><div><br /></div><div>決算で忙しかったので、記事が全然進んでいませんが、下書きは多少書き溜まっているので、随時出していこうと思います。</div><div>引き続きよろしくお願いいたします。</div><div><br /></div>管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-51571041469927559682020-06-06T16:55:00.000+09:002020-06-06T16:56:01.046+09:00R - 実践編 5 多変量時系列のモデル推定R の vars パッケージや dse パッケージを使って、多変量時系列を VAR モデルや状態空間モデルにフィッティングさせてみました。<br />
<br />
<a name='more'></a><h2>
対象データの基本統計量</h2>
<br />
今回は、多変量時系列データに対して、VAR モデルのパラメータを推定してみましょう。TS という多変量時系列データを用意しました。TS は、X と Y の2つの変量を持っています。基本統計量は以下の通り。
<br />
<pre class="prettyprint lang-r">> class(TS)
[1] "xts" "zoo"
> nrow(TS)
[1] 85496
> frequency(TS)
[1] 1
> summary(TS)
Index X Y
Min. :1970-01-01 00:00:00.00 Min. :-61884.00 Min. :-3.423e-04
1st Qu.:1970-01-01 05:56:13.75 1st Qu.: 0.00 1st Qu.: 0.000e+00
Median :1970-01-01 11:52:27.50 Median : 0.00 Median : 0.000e+00
Mean :1970-01-01 11:52:27.50 Mean : -0.01 Mean :-9.850e-08
3rd Qu.:1970-01-01 17:48:41.25 3rd Qu.: 0.00 3rd Qu.: 0.000e+00
Max. :1970-01-01 23:44:55.00 Max. : 65967.00 Max. : 4.028e-04 </pre>
frequency が 1 の xts データで、85496 行あります。秒刻みで1日分に少し足りないくらいのデータだと思ってください。X も Y も四分位点 Q1, Q2 (Median), Q3 がほぼ 0 なのに最大、最小の絶対値はずいぶん大きい、explosive な感じのするデータです。
<br />
<br />
<h2>
単位根検定</h2>
<div>
X, Y それぞれに対し、Philllips-Perron 検定による単位根検定をしてみます。PP.test() を使います。</div>
<pre class="prettyprint lang-r">> PP.test(TS$X)
Phillips-Perron Unit Root Test
data: TS$X
Dickey-Fuller = -637.5825, Truncation lag parameter = 21, p-value =
0.01
> PP.test(TS$Y)
Phillips-Perron Unit Root Test
data: TS$Y
Dickey-Fuller = -263.3456, Truncation lag parameter = 21, p-value =
0.01</pre>
Phillips-Perron 検定では、X, Y ともに p-value が小さく、単位根があるとは言い難いとの結論になっています。<br />
<br />
※ 実は X も Y も単位根過程らしき時系列の1階の階差です。<br />
<br />
<h2>
標本自己相関の構造</h2>
自己相関を確認しておきましょう。<br />
<pre class="prettyprint lang-r">> acf(TS, lag.max=10)</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-iybIkyWylNM/XttD4acQfNI/AAAAAAAAAOk/Qml4ZE8zXkQSph3vFGwq7jsx0JquLeEagCLcBGAsYHQ/s1600/vacf.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="480" data-original-width="640" height="480" src="https://1.bp.blogspot.com/-iybIkyWylNM/XttD4acQfNI/AAAAAAAAAOk/Qml4ZE8zXkQSph3vFGwq7jsx0JquLeEagCLcBGAsYHQ/s640/vacf.png" width="640" /></a></div>
<br />
<pre class="prettyprint lang-r">> pacf(TS, lag.max=10)</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-va0XTkwfCPY/XttD8YQiHoI/AAAAAAAAAOo/ImvqHB46f_8pWa5AZbcfO9diCKo-BltPACLcBGAsYHQ/s1600/vpacf.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="480" data-original-width="640" height="480" src="https://1.bp.blogspot.com/-va0XTkwfCPY/XttD8YQiHoI/AAAAAAAAAOo/ImvqHB46f_8pWa5AZbcfO9diCKo-BltPACLcBGAsYHQ/s640/vpacf.png" width="640" /></a></div>
<br />
acf() を見ると、X, Y ともに、3~4次あたりまで自己相関を持っているように見えます。減衰傾向もあります。pacf() を見ると、X と Y の間の偏自己相関が多少ありそうにも見えます。<br />
一応、VAR モデルを当ててみる価値はありそうです。<br />
<br />
<h2>
VAR モデルの選択</h2>
R で VAR モデルを扱うには、vars パッケージを使います。
VAR モデルの選択は VARselect() によって行うことができます。100次まで探索させてみましょう。
<br />
<pre class="prettyprint lang-r">> VARselect(TS, 100)
$selection
AIC(n) HQ(n) SC(n) FPE(n)
96 74 48 96
$criteria
1 2 3 4 5
AIC(n) -8.1011924695 -8.1383711520 -8.1530109267 -8.1631402110 -8.1740660487
HQ(n) -8.1009915707 -8.1380363206 -8.1525421628 -8.1625375145 -8.1733294197
SC(n) -8.1005351749 -8.1372756610 -8.1514772394 -8.1611683273 -8.1716559686
FPE(n) 0.0003031774 0.0002921126 0.0002878673 0.0002849661 0.0002818696
(略)
46 47 48 49 50
AIC(n) -8.2401103055 -8.2402657046 -8.2410426603 -8.241164891 -8.2413203039
HQ(n) -8.2338824420 -8.2339039086 -8.2345469317 -8.234535230 -8.2345567103
SC(n) -8.2197341738 -8.2194513765 -8.2197901358 -8.219474170 -8.2191913867
FPE(n) 0.0002638551 0.0002638141 0.0002636093 0.000263577 0.0002635361
(略)
71 72 73 74 75
AIC(n) -8.2458743893 -8.2463429275 -8.2464925460 -8.2466685854 -8.2466740306
HQ(n) -8.2362982121 -8.2366328178 -8.2366485037 -8.2366906106 -8.2365621233
SC(n) -8.2145433480 -8.2145736898 -8.2142851120 -8.2140229550 -8.2135902038
FPE(n) 0.0002623386 0.0002622157 0.0002621765 0.0002621304 0.0002621289
(略)
96 97 98 99 100
AIC(n) -8.2483858131 -8.2483204228 -8.2483074155 -8.2483632402 -8.2483824885
HQ(n) -8.2354613222 -8.2352619994 -8.2351150595 -8.2350369518 -8.2349222675
SC(n) -8.2060998623 -8.2055962756 -8.2051450719 -8.2047627003 -8.2043437522
FPE(n) 0.0002616806 0.0002616977 0.0002617011 0.0002616865 0.0002616815</pre>
$selection が各種情報量規準に照らしたときの最適次数を示しています。AIC では 96、SC (SIC と同じです) では 48 が最適次数とのことになりました。HQ はハンナン=クイン情報量規準という、また少しペナルティ関数の違う情報量規準です。FPE は赤池最終予測誤差のことで、AIC の原型的なものです。<br />
<br />
次数は 48 か 96 が良さそうとなったので、VAR(48), VAR(96) それぞれのフィッティングを行い、予測値を出させてみましょう。
<br />
<pre class="prettyprint lang-r">> TSvar48 <- VAR(TS, 48)
> TSvar96 <- VAR(TS, 96)
> predict(TSvar48)
$X
fcst lower upper CI
[1,] 28.90594185 -1603.502 1661.314 1632.408
[2,] -80.28080566 -1852.209 1691.648 1771.929
[3,] 83.28752937 -1712.583 1879.158 1795.871
[4,] 20.07676960 -1779.406 1819.560 1799.483
[5,] -28.53638886 -1831.497 1774.424 1802.961
[6,] -5.67486663 -1811.140 1799.790 1805.465
[7,] -5.53840932 -1811.336 1800.259 1805.798
[8,] 20.36609880 -1785.637 1826.369 1806.003
[9,] -14.88094661 -1821.103 1791.341 1806.222
[10,] -0.09496698 -1807.907 1807.717 1807.812
$Y
fcst lower upper CI
[1,] -4.179815e-07 -3.857079e-05 3.773483e-05 3.815281e-05
[2,] 1.886498e-06 -3.654026e-05 4.031326e-05 3.842676e-05
[3,] -1.274244e-06 -3.972669e-05 3.717821e-05 3.845245e-05
[4,] -2.815922e-07 -3.875275e-05 3.818957e-05 3.847116e-05
[5,] -7.214299e-07 -3.919800e-05 3.775514e-05 3.847657e-05
[6,] -7.825950e-07 -3.925958e-05 3.769439e-05 3.847698e-05
[7,] 1.102569e-07 -3.837391e-05 3.859442e-05 3.848416e-05
[8,] -2.166975e-07 -3.870116e-05 3.826777e-05 3.848447e-05
[9,] -5.326749e-08 -3.854202e-05 3.843549e-05 3.848875e-05
[10,] -6.753158e-07 -3.918076e-05 3.783013e-05 3.850545e-05
> predict(TSvar96)
$X
fcst lower upper CI
[1,] 33.770859 -1592.898 1660.440 1626.669
[2,] -63.865710 -1833.106 1705.375 1769.241
[3,] 77.636350 -1716.435 1871.708 1794.072
[4,] 14.859046 -1783.052 1812.770 1797.911
[5,] -94.458953 -1895.953 1707.035 1801.494
[6,] 6.612096 -1797.559 1810.783 1804.171
[7,] 55.755082 -1748.780 1860.290 1804.535
[8,] 40.833596 -1763.947 1845.614 1804.780
[9,] -74.713039 -1879.741 1730.315 1805.028
[10,] -10.848328 -1817.510 1795.813 1806.661
$Y
fcst lower upper CI
[1,] 1.514097e-06 -3.660781e-05 3.963600e-05 3.812191e-05
[2,] 3.693606e-06 -3.470163e-05 4.208884e-05 3.839524e-05
[3,] -1.648351e-06 -4.006920e-05 3.677250e-05 3.842085e-05
[4,] 1.439977e-06 -3.699968e-05 3.987963e-05 3.843965e-05
[5,] -2.497274e-06 -4.094177e-05 3.594722e-05 3.844450e-05
[6,] -2.431337e-07 -3.868807e-05 3.820180e-05 3.844493e-05
[7,] -4.750051e-07 -3.892673e-05 3.797672e-05 3.845172e-05
[8,] -5.696669e-07 -3.902169e-05 3.788235e-05 3.845202e-05
[9,] 1.452356e-07 -3.831102e-05 3.860149e-05 3.845625e-05
[10,] -6.050617e-07 -3.907691e-05 3.786678e-05 3.847184e-05</pre>
X の予測値は似てるけど、Y は全然違いますね。係数行列を見てみても、高次のパラメータがあまり減衰していません。
残差の自己相関は無くなっていますが、モデル化としてはあまり良くない結果に思えます。<br />
<br />
<h2>
状態空間モデルの選択</h2>
dse パッケージの estBlackBox() を用いると、状態空間モデルの選択をしてくれます。状態空間モデルでは、時系列を入力と出力に分けて考える必要があります。dse では、TSdata というオブジェクトで入出力のある時系列データを保持することができるので、まず TS から TSdata を作りましょう。今回の分析対象である TS は、実は Y が入力で X は応答であることを期待しています。以下のように input, output を指定して変換します。<br />
<pre class="prettyprint lang-r">> DSETS <- TSdata(input=TS$Y, output=TS$X)</pre>
<br />
そして、以下のようにモデルの選択を実行します。<br />
<pre class="prettyprint lang-r">> DSEEST <- estBlackBox(DSETS)
> summary(DSEEST)
neg. log likelihood = 697537.7 sample length = 85496
X
RMSE 845.3734
innovations form state space: nested model a la Mittnik
inputs : Y
outputs: X
input dimension = 1 state dimension = 12 output dimension = 1
theoretical parameter space dimension = 36
180 actual parameters 0 non-zero constants
Initial values not specified.</pre>
<br />
上記のようなモデルが選択されました。しかし、残差の ACF, PACF を見てみると、かなり強い自己相関が残っていて、モデル化できたとは言えませんでした。残念。管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-44979167221772780022020-05-28T10:26:00.000+09:002020-05-28T10:26:13.536+09:00R - 実践編 4 時系列のモデル推定時系列データに当てはまるモデルの推定を R で行ってみましょう。<a href="https://www.cyberer.net/2020/05/r-unit-root-test.html">実践編 3</a> で用いたものと同じ時系列 TS を対象として、モデルを推定した様子を紹介したいと思います。<br />
<br />
<a name='more'></a><h2>
基本統計量の確認</h2>
モデルの推定に先立って、基本統計量を確認しておく必要があります。平均、分散の値や単位根が無いこと、explosive であることは、<a href="https://www.cyberer.net/2020/05/r-unit-root-test.html">実践編 3</a> ですでに見ている通りなので、ここでは、標本自己相関関数、標本偏自己相関関数を確認します。<br />
R には、acf() という関数が用意されています。時系列データから標本自己相関関数を計算、返却しつつ、コレログラムをプロットしてくれる関数です。また、pacf() 関数もあり、こちらは標本偏自己相関関数の計算・プロットです。なお、acf() に type = "partial" を指定しても pacf() を呼んでくれます。<br />
<br />
どの程度の次数まで見れば良いか分からないので、とりあえずデータの件数分全部見てみることにします。<br />
<pre class="prettyprint lang-r">> acf(TS, lag.max=nrow(TS))
</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-nxEJnuKxoC4/Xs8LYKOhGiI/AAAAAAAAAOM/2OQHFhI3PbwWxWhqp8JzK1yC1_z1NSsxACLcBGAsYHQ/s1600/acf.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="480" data-original-width="640" height="480" src="https://1.bp.blogspot.com/-nxEJnuKxoC4/Xs8LYKOhGiI/AAAAAAAAAOM/2OQHFhI3PbwWxWhqp8JzK1yC1_z1NSsxACLcBGAsYHQ/s640/acf.png" width="640" /></a></div>
<pre class="prettyprint lang-r">> pacf(TS, lag.max=nrow(TS))</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-oO2Wj8vAc28/Xs8Lfuep96I/AAAAAAAAAOQ/taqbPupxeXYqseWDa_00xwWPy8ojkXI4wCLcBGAsYHQ/s1600/pacf.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="480" data-original-width="640" height="480" src="https://1.bp.blogspot.com/-oO2Wj8vAc28/Xs8Lfuep96I/AAAAAAAAAOQ/taqbPupxeXYqseWDa_00xwWPy8ojkXI4wCLcBGAsYHQ/s640/pacf.png" width="640" /></a></div>
ちょっと見えづらいですが、青破線が「自己相関や偏自己相関が 0 である」という帰無仮説の有意水準5%での棄却点を示しています。標本自己相関関数は長期間に渡ってゆるやかに振動しながら減衰しており、標本偏自己相関関数は低い水準で激しく振動しています。少ない次数では偏自己相関があるようにも見えます。<br />
<br />
forecast パッケージには、tsdisplay() 関数があり、これは原系列と acf, pacf を1つのプロットとして出力するものです。<br />
<pre class="prettyprint lang-r">> tsdisplay(TS)</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-3xvMCe3RYAQ/Xs8LnWBj_NI/AAAAAAAAAOU/lpzqGvg7YtcXdWrKfObsiEOm8zLmTnxhgCLcBGAsYHQ/s1600/tsdisplay.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="480" data-original-width="640" height="480" src="https://1.bp.blogspot.com/-3xvMCe3RYAQ/Xs8LnWBj_NI/AAAAAAAAAOU/lpzqGvg7YtcXdWrKfObsiEOm8zLmTnxhgCLcBGAsYHQ/s640/tsdisplay.png" width="640" /></a></div>
原系列は諸事情によりボカしています。すみません。lag.max を指定していないので、適当なラグが選択されています。これで見ると標本偏自己相関関数については50次あたりまで減衰しながらの正の相関を持っているように見えています。<br />
<br />
<h2>
モデルの選択</h2>
基本統計量の概観を元に、モデルを選択してみます。自己相関に切断が見られないため、MA モデルは適しないと思われます。切断無く減衰傾向があることから、AR または ARMA モデルの適用ができる可能性があります。また、偏自己相関にも切断は見られないため、AR モデルよりは ARMA モデルの適用を試みる価値が高いのではないかと思われます。<br />
ARMA モデルの場合、次数の目安はここからは分からないので、いくつかフィッティングを行った上で情報量規準による選択を行うことになります。<br />
<br />
<h2>
ARMA (ARIMA) モデルのフィッティング</h2>
stats::ar() は AR モデルの次数を AIC に基づいて選択する機能を持っているのですが、ARMA モデルについての同様の関数は標準にはありません。<br />
ARMA モデルの最適次数の選択は、forecast パッケージにある auto.arima() を使います。この関数は自己回帰<b>和分</b>移動平均 (Autoregressive <b>integrated</b> moving average / ARIMA) モデルへのフィッティングを行うものです。ARIMA モデルは非定常過程のためのモデルで、階差に対して ARMA モデルを適用するものになっています。$ARIMA(p, d, q)$は$d$階の階差に$ARMA(p, q)$を適用するモデルを意味します。$ARMA(p, q)$は$ARIMA(p, 0, q)$と同等なので、auto.arima() が選ぶモデルの中には ARMA が選択対象に含まれています。ちなみに$p, q$についても$0$が選択され得るので、AR, MA も含まれています。つまり、ARIMA モデルは、AR, MA, ARMA の選択や階差の取り方を全てパラメータとして押し込んだモデルになっているので、ARIMA のパラメータ推定を情報量規準に基づいて行うことで、モデルの選択も含めてワンストップでできることになります。<br />
auto.arima() を実行してみます。<br />
<pre class="prettyprint lang-r">> TSarima <- auto.arima(coredata(TS))
> TSarima
Series: coredata(TS)
ARIMA(3,1,2)
Coefficients:
ar1 ar2 ar3 ma1 ma2
0.7562 -0.0543 0.0405 -1.1877 0.1992
s.e. 0.0705 0.0390 0.0078 0.0706 0.0694
sigma^2 estimated as 692537: log likelihood=-696194.6
AIC=1392401 AICc=1392401 BIC=1392457</pre>
<br />
auto.arima() は配列を渡さないとうまく動作しないので、coredata(TS) を渡しています。次数$p, q$は初期値では$2$から$5$の値を試すようになっています。この範囲では、$ARIMA(3, 1, 2)$が最適ということになったようです。<br />
もう少し次数の幅を広げて試してみましょう。<br />
<pre class="prettyprint lang-r">> TSarima <- auto.arima(coredata(TS), start.p = 0, max.p = 50, start.q = 0, max.q = 50)
> TSarima
Series: coredata(TS)
ARIMA(1,1,3)
Coefficients:
ar1 ma1 ma2 ma3
0.7719 -1.2032 0.1512 0.0619
s.e. 0.0056 0.0067 0.0057 0.0048
sigma^2 estimated as 692511: log likelihood=-696193
AIC=1392396 AICc=1392396 BIC=1392443
</pre>
<br />
$ARIMA(1, 1, 3)$が、$ARIMA(3, 1, 2)$よりもわずかに AIC の値が小さく、最適とのことになりました。ic="bic" を付けると、AIC にかえて BIC で評価するようになるが、これでも同じく$ARIMA(1, 1, 3)$が選択されました。<br />
従って、$ARIMA(1, 1, 3)$がモデル化の有力な候補ということになります。<br />
<br />
<h2>
モデルの診断</h2>
$ARIMA(1, 1, 3)$が一番フィッティングが良かったからといって、これが正解のモデルとは限りません。あくまでモデルの候補が見付かったに過ぎません。モデルの診断を行う必要があります。自己相関のモデル化をしているわけですから、残差に自己相関が無いことを確認することになります。<a href="https://www.cyberer.net/2020/05/timeseries-analysis.html#toc_headline_17">こちらの記事</a>でも取りあげましたが、かばん検定といいます。R では Box.test(data, lag=m, fitdf=k) という関数が用意されています。$Q(m)$値が$\chi^2(m-k)$で検定されます。<br />
今回は$ARIMA(1, 1, 3)$なのでパラメータ数$k$は$5$です。$Q(6)$から$Q(50)$まで確認してみることにします。<br />
<pre class="prettyprint lang-r">> sapply(6:50, function(l) Box.test(TSarima$residuals, lag=l, type="Ljung-Box", fitdf=5)$p.value)
[1] 6.119297e-07 2.136468e-07 8.512888e-08 0.000000e+00 0.000000e+00
[6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[26] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[31] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[41] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
</pre>
<br />
p値が極めて小さく、Ljung-Box 検定の帰無仮説である「自己相関が無いこと」は棄却されてしまいました。つまり自己相関が残っているとみられます。<br />
<br />
結果として、ARIMA モデルへのフィッティングは失敗だったということになります。acf, pacf から自己相関があることはほぼ間違いないので、ARIMA の範囲でモデル化できるような構造の自己相関ではなかったということです。より複雑な他のモデルを検討する必要があるでしょう。<br />
<br />
<br />管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-713833790475973652020-05-28T09:39:00.001+09:002020-05-28T09:41:29.934+09:00R - 実践編 3 時系列の単位根検定R で実際に時系列データの単位根検定をしてみましょう。Phillips-Perron 検定と Augmented Dickey-Fuller 検定が使用できます。<br />
<br />
<a name='more'></a><h2>
ランダムウォークと単位根</h2>
ランダムウォークと単位根については、詳しくは<a href="https://www.cyberer.net/2020/05/timeseries-analysis.html">こちらの記事</a>を見ていただければと思います。簡単に触れておくと、時系列データを$y_t = ay_{t-1} + \epsilon_t, \epsilon_t \sim iid(0, \sigma^2)$とモデル化できるとした場合、$|a| = 1$であれば、この$y_t$はランダムウォークと呼ばれます。また同時に、この時系列は「単位根を持つ」と言われます。<br />
ランダムウォークは非定常過程ですが、単位根を持つ場合、1階の階差を取ると$iid$過程になるため、弱定常過程に変化します。<br />
<br />
<h2>
単位根検定</h2>
時系列がランダムウォークであるかどうか (すなわち単位根を持つかどうか) を調べることを単位根検定と言います。単位根検定とは、「$|a| = 1$である」という帰無仮説を置いた帰無仮説検定です。<br />
R を使って時系列の単位根検定を行ってみましょう。以下のような時系列 TS を題材にします。<br />
<pre class="prettyprint lang-r">> class(TS)
[1] "xts" "zoo"
> names(TS)
[1] "X"
> summary(TS)
Index X
Min. :****-**-** **:**:** Min. :-77287.0
1st Qu.:****-**-** **:**:** 1st Qu.: -927.0
Median :****-**-** **:**:** Median : -348.0
Mean :****-**-** **:**:** Mean : -355.2
3rd Qu.:****-**-** **:**:** 3rd Qu.: 366.0
Max. :****-**-** **:**:** Max. : 70648.0
> length(TS)
[1] 85497
> sqrt(var(TS))
X
X 1481.54
</pre>
<br />
TS は xts クラスの時系列データで、X という列1つだけが含まれています。X の分布は summary(TS) の結果として表示されている通りです。件数は85,497件、標準偏差は1481.54となっています。たかだか85,497件しかないのに、$-52.17\sigma\sim47.69\sigma$の範囲に幅広く分布しています。なかなか恐ろしいデータです。<br />
<br />
<h3>
Phillips-Perron 検定</h3>
R では、単位根検定の一種である、Phillips-Perron 検定を行う関数として、PP.test() が用意されています。stats パッケージに入っています。実行してみましょう。<br />
<pre class="prettyprint lang-r">> PP.test(TS)
Phillips-Perron Unit Root Test
data: TS
Dickey-Fuller = -117.1077999999999974534, Truncation lag parameter =
21, p-value = 0.01</pre>
<br />
p-value = 0.01 なので、「単位根がある」という帰無仮説は棄却されたことになります。つまり、TS には単位根は無いであろう、という結果です。<br />
<br />
<h3>
Augmented Dickey-Fuller 検定</h3>
tseries パッケージには、Augmented Dickey-Fuller 検定を行う adf.test() 関数があります。帰無仮説は「単位根を持つ」で、対立仮説は「$|a| < 1$ (stationary)」か「$|a| > 1$ (explosive)」を選ぶことができます。実行してみましょう。<br />
<pre class="prettyprint lang-r">> adf.test(TS)
Augmented Dickey-Fuller Test
data: TS
Dickey-Fuller = 210.4179000000000030468, Lag order = 44, p-value = 0.99
alternative hypothesis: stationary
警告メッセージ:
In adf.test(TS) : p-value greater than printed p-value
</pre>
<br />
デフォルトでは対立仮説は stationary になっています。こちらでは、p-value = 0.99 ですので、帰無仮説が棄却されませんでした。対立仮説を explosive に変えてみましょう。<br />
<pre class="prettyprint lang-r">> adf.test(TS, alternative="explosive")
Augmented Dickey-Fuller Test
data: TS
Dickey-Fuller = 210.4179000000000030468, Lag order = 44, p-value = 0.01
alternative hypothesis: explosive
警告メッセージ:
In adf.test(TS, alternative = "explosive") :
p-value greater than printed p-value
</pre>
<br />
今度は p-value = 0.01 で帰無仮説は棄却されました。併せて考えると、「単位根は無い」ということになり、PP.test() と同じ結論になりました。<br />
<br />管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-48171933094526828732020-05-28T08:43:00.001+09:002020-05-28T08:50:24.425+09:00R - 実践編 2 データフレームにおける計算、集計分析対象データをデータフレームに読み込んだあと、本格的な分析に入る前に、計算や集計などを行う必要があることが多いかと思います。R での代表的なやり方を説明します。<br />
<br />
<a name='more'></a><h2>
既存列に基づく計算結果を新しい列に持たせる</h2>
データフレーム df の列 a, b を元に、何らかの計算 f() を行ってその結果を列 c に持たせるのであれば、基本的には以下のようにします。<br />
<pre class="prettyprint lang-r">df$c <- f(a, b)
</pre>
<br />
この方法では、f が同じ長さのベクトル2つを引数にとり、同じ長さのベクトルを結果として返すように定義されていることを想定しています。R では多くの算術演算子や算術関数がそのように振舞いますので、こういった簡単な書き方が可能なのです。<br />
<br />
1 ステップで書きづらいものや、無理に 1 ステップで書こうと思うと効率が悪くなるものは、次のような方法もあります。例えば、データフレーム df の列 a の値をキーとして、別のデータフレーム df2 から行を選択し、その行から列 b の値に応じて df2 の列 X か列 Y かのいずれかの値を取って列 c に格納する、といった例を考えてみましょう。<br />
<pre class="prettyprint lang-r">df$c <- local({
tmp <- df2[df2$a == df$a,]
ifelse(df$b == "X", tmp$X, tmp$Y)
})
</pre>
<br />
df2 から選択した行は一度 tmp に付値しておくことで、df2 への添字ベクトルの適用を一度で済ませています。そうしないと、列 X や列 Y を取得するときにも添字ベクトルによる行選択が必要となり、非効率です。local() は、識別子 tmp が名前空間を汚さないようにするためのものです。<br />
<br />
<h2>
二分探索法</h2>
別の例を見てみましょう。「df2\$a == df\$a」という同値条件に代えて、「df\$a の値を越えない最大の df2\$a の値を持つ行」を抽出してみましょう。<br />
単純に考えると、以下のようになります。<br />
<pre class="prettyprint lang-r">df$c <- local({
fun <- function(x) df2[which(df2$a == max(df2$a[df2$a <= x])),]
tmp <- lapply(df$a, fun)
ifelse(df$b == "X", tmp$X, tmp$Y)
})
</pre>
<br />
ですがこの方法は、df2 からの行抽出が何度も行われ、とても効率が悪いです。<br />
こういった場合は、df2 を列 a をキーに整列しておき、二分探索法を使います。まず以下のように二分探索を行う関数 bsearch() を作りましょう。標準で合ってもよさそうなものなのですが、見当たらないので……。<br />
<pre class="prettyprint lang-r">bsearch <- function(val, tab, L=1L, H=length(tab)) {
b <- cbind(L=rep(L, length(val)), H=rep(H, length(val)))
i0 <- seq_along(val)
repeat {
updt <- M <- b[i0, "L"] + (b[i0, "H"] - b[i0, "L"]) %/% 2L
tabM <- tab[M]
val0 <- val[i0]
i <- tabM < val0
updt[i] <- M[i] + 1L
i <- tabM > val0
updt[i] <- M[i] - 1L
b[i0 + i * length(val)] <- updt
i0 <- which(b[i0, "H"] >= b[i0, "L"])
if (!length(i0)) break;
}
b[,"L"] - 1L
}</pre>
<br />
この bsearch() を用いて、次のようにすると、とても速いです。<br />
<pre class="prettyprint lang-r">df$c <- local({
tmp <- df2[bsearch(df$a, df2$a),]
ifelse(df$b == "X", tmp$X, tmp$Y)
})</pre>
<br />
<h2>
データフレームの集計</h2>
データフレームを集計するには by() を使います。例えば df\$d に何らかの因子が入っているとして、その水準ごとに df\$c の平均、標準偏差、個数を集計するとしましょう。<br />
<pre class="prettyprint lang-r">summ <- by(df$c, df$d, function(x) { x <- x[[1]]; c(mean(x), sd(x), length(x)) })
</pre>
<br />
by() の戻り値は by クラスのオブジェクトで、df\$d の水準名が要素名となり、function(x) ... の戻り値が値となったリストの形をしています。<br />
<br />管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-27122428902318778162020-05-28T07:48:00.003+09:002020-05-28T07:48:32.781+09:00R - 実践編 1 ファイルの読み込み、データフレームの加工実際に R でデータ解析を行うために、CSV ファイル、TSV ファイルの読み込みを行ったときの知見をまとめました。<br />
<br />
<a name='more'></a><h2>
CSV ファイルの読み込みは read.csv()</h2>
CSV ファイルの読み込みは、read.csv() を用いて行います。read.csv() は、read.table() の引数の既定値を CSV ファイル向けに変更したものなので、動作は read.table() そのものです。read.table() については<a href="https://www.cyberer.net/2020/05/r-reading-from-data-files.html#toc_headline_3">こちらの記事</a>を見てください。引数の既定値は、以下のように変更されています。<br />
<br />
<ul>
<li>header = TRUE (1行目は列数によらずヘッダ行と見做す)</li>
<li>sep = "," (列のセパレータはカンマ)</li>
<li>quote = "\"" (引用符はダブルクオートのみ)</li>
<li>fill = TRUE (列数が異なっていても、補完する)</li>
<li>comment.char = "" (コメントは使えない)</li>
</ul>
<div>
以下の説明は、read.csv() を使ったときの知見の紹介ですが、多くは read.table() にも当て嵌まります。</div>
<div>
<br /></div>
<h3>
列のモードの指定</h3>
列のモードは read.table() の仕様に従って適宜自動的に判断されますが、それが望ましくない場合は colClasses で各列のクラスを指定することができます。<br />
一旦は何も考えずに読み込んで、後から強制変換すれば良いだろうとも思っていたのですが、実際にやってみると、こんな問題がありました。<br />
<br />
<h4>
桁数の多い数値列</h4>
桁数が多い数字列は、double として読み込まれてしまうのですが、このとき倍精度実数の最大精度を越えていると、丸められてしまうということが分かりました。<br />
多くの場合、桁数の多い数字列というのは、数値として扱いたいわけではなく、コード値のような扱いをするものだと思います。そういう場合は、colClasses で character や factor に指定しておけば、丸められずに読み込めます。本当に大きな数値を扱いたい場合は、いわゆる多倍長精度数値を扱うクラスがあるのではないかと想像しますが、今のところ知りません。すみません。<br />
<br />
<h4>
日付・日時列</h4>
日付、日時の列ならば colClasses に Date クラスが指定できます。ただフォーマットの指定は colClasses ではやりようが無いので、フォーマットが合わないときは、character で読み込んでおいて後から format.Date() で直せば良いでしょう。他の方法としては、既定のフォーマットを変更してから実行するとか、独自の Date クラスのサブクラスを作るといったことが考えられます。<br />
<br />
<h4>
ミリ秒まで持っている時刻データ列</h4>
ミリ秒精度だと、Date クラスでは保持できないので、POSIXct を指定します。ミリ秒も読み込まれます。表示上は見えませんが、options(digits.secs=3) とすればミリ秒まで見えてきます。ただ、精度の問題はあるようで、例えば 2020-05-05 00:00:00.123 を POSIXct に読み込ませて、options(digits.secs=3) の状態で表示させると、2020-05-05 00:00:00.122 となってしまいました。これはちょっと問題ですね。<br />
<br />
<h4>
colClasses に列名を指定する</h4>
colClasses に与えるベクトルを名前付きベクトルにすれば、列名で指定できます。<br />
<br />
<h3>
試行錯誤中は nrows で行数を限定する</h3>
read.csv() のパラメータの指定は割と試行錯誤になってしまうので、対象のファイルが大きいと、試す都度時間がかかって大変です。<br />
と、思っていたら、nrows パラメータで最大行数が指定できるので、ここを少な目の行数に設定して試せば早いですね。もっと早く気付きたかったです。<br />
<br />
<h2>
TSV ファイルの読み込みは read.delim()</h2>
TSV ファイルの読み込みは、read.delim() を用います。read.csv() と同様、read.table() の引数の既定値を変更しただけのものです。<br />
<br />
<h2>
データフレームの加工</h2>
とりあえずデータフレームに読み込めたとして、実際に分析を開始するまでには、通常はある程度加工が必要となるでしょう。<br />
<br />
<h3>
列の抽出は添字ベクトル</h3>
<div>
解析に必要な行、あるいは必要な列だけを抽出したいということはよくあります。</div>
<div>
列の抽出だけであれば添字ベクトルを使って、df <- df[select] とすることができます。select は添字ベクトルと呼ばれており、列名のベクトルです。</div>
<pre class="prettyprint lang-r">> df[c("x", "y")]
x y
1 1 4
2 2 5
3 3 6
</pre>
<br />
<h3>
行・列を一度に絞り込むなら subset()</h3>
必要な行、列を一度に絞り込むには subset(df, subset, select) が使えます。データフレーム df から、論理式 subset を満たす行を抽出し、文字列ベクトル select に含まれる列名の列を抽出する、という動作をします。<br />
<pre class="prettyprint lang-r">> df <- data.frame(x = c(1, 2, 3), y = c(4, 5, 6), z = c(7, 8, 9))
> subset(df, x %in% c(1, 3))
x y z
1 1 4 7
3 3 6 9
> subset(df, select = c("x", "z"))
x z
1 1 7
2 2 8
3 3 9
> subset(df, y %in% c(5, 6), c("y", "z"))
y z
2 5 8
3 6 9
</pre>
<br />
<h3>
ファイル読み込み時に行や列を除外する</h3>
ファイルを全部読み込んでから列や行を削るのは効率が悪いので、読み込み時点である程度除外できると、効率が良いです。<br />
行については、nrows (読み込む行数) と skip (指定した行数より前の行を読み飛ばす) を指定することで、ある程度対応できます。<br />
列については、colClasses パラメータで NULL を指定した列はスキップされます。<br />
<br />
<h3>
列のクラスを変更する</h3>
読み込んだデータフレームの列のクラスを後から変更したいこともよくあります。この記事でも、一旦 character で読み込んでおいて、後でどうにかしましょう、という話をしていました。<br />
実際には、例えば3列目を character から factor にしたければ以下のようにします。<br />
<pre class="prettyprint lang-r">> df[, 3] <- factor(df[, 3])
</pre>
<br />
複数列まとめて変更したければ、lapply() を使うことができます。3列目と74列目を character から factor に変更するなら、こうです。<br />
<pre class="prettyprint lang-r">> df[, c(3, 74)] <- lapply(df[, c(3, 74)], factor)
</pre>
<br />
<h3>
データフレームの結合は rbind()</h3>
データが複数ファイルに分割されている場合など、複数のデータフレームに読み込まれたデータを1つのデータフレームに結合したいこともよくあります。こういったときは、rbind() を使うことができます。<br />
例えば、filenames に読み込みたいファイル名のベクトルがあって、これらを全て読み込んで1つのデータフレームに結合したい場合は、次のようにすることができます。<br />
<pre class="prettyprint lang-r">> df <- data.frame()
> for (filename in filenames) {
+ dftmp <- read.delim(filename)
+ # 必要に応じて dftmp の subset() をとったりなんだり
+ df <- rbind(df, dftmp)
+ }</pre>
<br />
最初に空のデータフレーム df を用意しておき、そこにファイルから読み込んで、加工した上で df に結合していきます。<br />
加工が必要なければ以下のように、よりシンプルにできます。<br />
<pre class="prettyprint lang-r">> df <- Reduce(rbind, Map(read.delim, filenames))
</pre>
<br />
read.delim のかわりに加工も含めた関数を与えれば、前者と同じこともできますね。管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-43499499518009534002020-05-27T05:36:00.004+09:002020-05-29T09:46:28.893+09:00都心にカブトムシいた!<div class="separator" style="clear: both; text-align: left;">
先日、千鳥ヶ淵付近を散歩していたら、カブトムシを発見しました。</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-xJYndqcnFuc/Xs17u0rb7zI/AAAAAAAAAOA/uSbBwvmxQe0eQH93gWMJcX8NLHCs9ATIwCLcBGAsYHQ/s1600/IMG_1036.jpeg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img alt="カブトムシ" border="0" data-original-height="400" data-original-width="640" height="200" src="https://1.bp.blogspot.com/-xJYndqcnFuc/Xs17u0rb7zI/AAAAAAAAAOA/uSbBwvmxQe0eQH93gWMJcX8NLHCs9ATIwCLcBGAsYHQ/s320/IMG_1036.jpeg" title="カブトムシ" width="320" /></a></div>
<br />
こんな都心部にカブトムシがいるとは、と思いましたが、良く考えたら靖国神社や皇居にはいっぱい居るはずですね。そのあたりから飛んできたんでしょうか。それよりも、5月にこんなしっかりしたサイズのカブトムシが居ることの方が珍しいかもしれません。<br />
<br />
いずれにしても珍しいものを見ました。<br />
<br />管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-15260640543157691142020-05-27T04:36:00.001+09:002020-05-27T04:36:52.927+09:00R - 入門本編 7章 ファイルからデータを読み込む実践的には、プロンプトから data.frame() だの c() だのを駆使してデータフレームを構築することはまずありません。当然ながら、外部ファイルからデータを読み込んでデータフレームを構成することができますので、その方法を習得しておきましょう。<br />
<div>
<br /></div>
<a name='more'></a><h2>
read.table() 関数</h2>
read.table() は、特定の形をした外部ファイルを読み込んでデータフレームを生成する関数です。外部ファイルの形式としては、以下のようなことが要請されます。<br />
<ul>
<li>1行目が以下の条件を満たす場合、1行目はデータフレームの列名定義行になり、2行目以降がデータ行になる</li>
<ul>
<li>1行目の列数がデータ行の列数より1つ少ないこと</li>
</ul>
<li>そうでない場合は、1行目もデータ行とみなされ、列名は V1, V2, ..., と自動的に割り当てられる</li>
<li>各行はいわゆるホワイトスペース (1つ以上の半角スペース、タブ、CR, LF) で区切られた項目の集まりとして認識される</li>
<li>データ行の1列目は行ラベルとして認識される</li>
</ul>
<div>
例えば、houses.txt という外部ファイルがあるとしましょう。内容は以下のようなものです。</div>
<pre class="prettyprint lang-xml"> Price Floor Area Rooms Age Cent.heat
01 52.00 111.0 830 5 6.2 no
02 54.75 128.0 710 5 7.5 no
03 57.50 101.0 1000 5 4.2 no
04 57.50 131.0 690 6 8.8 no
05 59.75 93.0 900 5 1.9 yes</pre>
<div>
このファイルを読み込んでデータフレームを生成するには、以下のようにします。</div>
<pre class="prettyprint lang-r">> HousePrice <- read.table("houses.txt")
> HousePrice
Price Floor Area Rooms Age Cent.heat
01 52.00 111 830 5 6.2 no
02 54.75 128 710 5 7.5 no
03 57.50 101 1000 5 4.2 no
04 57.50 131 690 6 8.8 no
05 59.75 93 900 5 1.9 yes</pre>
デフォルトの振舞いでは、行ラベルを除く数値項目は数値変量として読み込まれます。そうではない項目 (この例では Cent.heat) は文字列として読み込まれます。
<br />
<h4>
(補足) バージョンによる違い</h4>
文字列として読み込まれるようになったのはバージョン4からです。過去のバージョンでは、因子として読み込まれていました。stringsAsFactors 引数を TRUE と指定すれば、従来通り因子になります。options(stringsAsFactors = TRUE) とする方法もありますが、非推奨です。<a href="https://www.cyberer.net/2020/05/r-version-4-released.html#toc_headline_6">こちらも合わせてご覧ください</a>。<br />
<br />
データフレームの実体はリストです。リストの各要素に class() を適用した結果と、全体を unclass() した結果を見てみましょう。<br />
<pre class="prettyprint lang-r">> lapply(HousePrice, class)
$Price
[1] "numeric"
$Floor
[1] "numeric"
$Area
[1] "integer"
$Rooms
[1] "integer"
$Age
[1] "numeric"
$Cent.heat
[1] "character"
> unclass(HousePrice)
$Price
[1] 52.00 54.75 57.50 57.50 59.75
$Floor
[1] 111 128 101 131 93
$Area
[1] 830 710 1000 690 900
$Rooms
[1] 5 5 5 6 5
$Age
[1] 6.2 7.5 4.2 8.8 1.9
$Cent.heat
[1] "no" "no" "no" "no" "yes"
attr(,"row.names")
[1] "01" "02" "03" "04" "05"</pre>
ファイルの1行目の列がリストのコンポーネント名になってています。lapply(HousePrice, class) の結果から、各列がどのクラスになったかが分かります。Cent.heat は character (文字列) になり、他の列は numeric (数値) になりました。データはベクトルで格納されています。行ラベルは row.names という属性に格納されていて、リストのコンポーネントにはなっていません。<br />
<br />
行ラベルを付けていないデータも扱うことができます。houses.txt から行ラベルを取り除いた houses2.txt を用意して、データフレームを生成してみます。houses2.txt の内容は次の通りです。
<br />
<pre class="prettyprint lang-xml">Price Floor Area Rooms Age Cent.heat
52.00 111.0 830 5 6.2 no
54.75 128.0 710 5 7.5 no
57.50 101.0 1000 5 4.2 no
57.50 131.0 690 6 8.8 no
59.75 93.0 900 5 1.9 yes
</pre>
結論から先に言うと、この形式のファイルを期待通りの形でデータフレームに読み込むには、header=TRUE 引数を指定する必要があります。
<br />
<pre class="prettyprint lang-r">> HousePrice2 <- read.table("houses2.txt", header=TRUE)
> HousePrice2
Price Floor Area Rooms Age Cent.heat
1 52.00 111 830 5 6.2 no
2 54.75 128 710 5 7.5 no
3 57.50 101 1000 5 4.2 no
4 57.50 131 690 6 8.8 no
5 59.75 93 900 5 1.9 yes
> lapply(HousePrice2, class)
$Price
[1] "numeric"
$Floor
[1] "numeric"
$Area
[1] "integer"
$Rooms
[1] "integer"
$Age
[1] "numeric"
$Cent.heat
[1] "character"</pre>
header=TRUE は、1行目が列名を表す行であることを明示するための引数です。
houses2.txt は、全行の列数が等しいため、デフォルトの振舞いでは1行目もデータ行と見做され、期待した構造では読み込まれません。試しにやってみると、次のようになります。<br />
<pre class="prettyprint lang-r">> HousesNG <- read.table("houses2.txt")
> HousesNG
V1 V2 V3 V4 V5 V6
1 Price Floor Area Rooms Age Cent.heat
2 52.00 111.0 830 5 6.2 no
3 54.75 128.0 710 5 7.5 no
4 57.50 101.0 1000 5 4.2 no
5 57.50 131.0 690 6 8.8 no
6 59.75 93.0 900 5 1.9 yes
> lapply(HousesNG, class)
$V1
[1] "character"
$V2
[1] "character"
$V3
[1] "character"
$V4
[1] "character"
$V5
[1] "character"
$V6
[1] "character"</pre>
1行目がデータ行になり、列名の指定が無いことになるので、列名は V1, V2, ..., V6 が割り当てられます。1行目がデータ行になったことによって、どの列も文字列要素が混じっていることになり、全列とも文字列扱いになってしまいました。<br />
<br />
<h2>
scan() 関数</h2>
scan() は、より低水準の外部ファイル読み込み関数です。データ行の項目数が均一で、各項目のモードも事前に明確であれば、scan() を用いて項目ごとのベクトルとして読み込むことができます。<br />
<br />
例えば上の houses.txt を scan() で読み込むとしたら、以下のようになります。<br />
<pre class="prettyprint lang-r">> HousePriceScan <- scan("houses.txt", list("", Price=0, Floor=0, Area=0, Rooms=0, Age=0, Cent.heat=""), skip=1)
Read 5 records
> HousePriceScan
[[1]]
[1] "01" "02" "03" "04" "05"
$Price
[1] 52.00 54.75 57.50 57.50 59.75
$Floor
[1] 111 128 101 131 93
$Area
[1] 830 710 1000 690 900
$Rooms
[1] 5 5 5 6 5
$Age
[1] 6.2 7.5 4.2 8.8 1.9
$Cent.heat
[1] "no" "no" "no" "no" "yes"</pre>
<br />
skip=1 は先頭から1行読み飛ばす指定です。scan() では列名定義行が扱えないので、読み飛ばしています。2つ目の引数に与えているリストは、各項目のモードを scan() に教えるためのダミーリストです。このリストのコンポーネント名が、生成されるリストのコンポーネント名に引き継がれます。<br />
リストではなくベクトルを与えることもでき、その場合は結果はベクトルで返ります (改行は関係なくなります)。ファイルに含まれる値の数の分、この引数のベクトルの要素が反復適用されて読み込まれます。例えばスカラーの 0 を指定すれば、全ての列が数値項目であることが期待され、その通りであれば数値ベクトルが返ります。<br />
<br />
<h2>
組込みデータセットにアクセスする</h2>
R には組み込みデータセットというものが用意されています。基本的なものは datasets パッケージに含まれており、他のパッケージも組み込みデータセットを持っていることがあります。datasets パッケージ内のデータセットの一覧は、data() で見ることができます。他のパッケージのデータセットの一覧であれば data(package = "boot") のようにパッケージを指定すれば確認できます。data(package = .packages(all.available = TRUE)) とすれば、全ての利用可能なパッケージからデータセットを探して一覧にしてくれます。<br />
データセットを読み込むには、data(infert) のようにデータセット名を指定します。datasets 以外のパッケージから読み込む場合は、data(amis, package="boot") のように package タグを指定してください。データセットを読み込むと、1つ以上の R のオブジェクトが読み込まれます。通常はデータセットと同名のデータフレーム1つであることがほとんどですが、そうでないこともあります。何が読み込まれるかはデータセット名でヘルプを引けばわかるはずです。なお、library() 関数を用いてパッケージを検索リストに入れれば、そのパッケージ内のデータセットは package タグで指定することなく読み込めます。<br />
<br />
さらに言えば、検索リスト上に存在するデータセットは、読み込まなくてもその名称で参照することは可能です。data() での読み込みは、データセットの複製をデフォルトの名前空間上に作っていることになります。変更を加えるような付値を行っても自動的に読み込まれます。<br />
<pre class="prettyprint linenums lang-r">> ls()
character(0)
> airmiles
Time Series:
Start = 1937
End = 1960
Frequency = 1
[1] 412 480 683 1052 1385 1418 1634 2178 3362 5948 6109 5981
[13] 6753 8003 10566 12528 14760 16769 19819 22362 25340 25343 29269 30514
> ls()
character(0)
> airmiles[1] <- 0
> airmiles
Time Series:
Start = 1937
End = 1960
Frequency = 1
[1] 0 480 683 1052 1385 1418 1634 2178 3362 5948 6109 5981
[13] 6753 8003 10566 12528 14760 16769 19819 22362 25340 25343 29269 30514
> ls()
[1] "airmiles"
> rm(airmiles)
> ls()
character(0)
> airmiles
Time Series:
Start = 1937
End = 1960
Frequency = 1
[1] 412 480 683 1052 1385 1418 1634 2178 3362 5948 6109 5981
[13] 6753 8003 10566 12528 14760 16769 19819 22362 25340 25343 29269 30514</pre>
<br />
最初の ls() で、デフォルトの名前空間には何も付値されていないことが分かります (1〜2行目)。airmiles という名前を参照すると、ts クラスのデータが参照されました (3〜9行目)。これは datasets パッケージに定義されています。datasets パッケージは検索リスト上にあるので、直接参照することができます。参照後に ls() をしても、やはり付値されていません (10〜11行目)。参照しただけでは複製は作られません。12行目で最初の要素を 0 に書き換えてみました。airmiles を参照してみると、書き変わっていることが分かります (13〜19行目)。ここで改めて ls() してみると、airmiles という名前で付値されていることが確認できます (20〜21行目)。rm(airmiles) として消してみると、デフォルトの名前空間から消えます (22〜24行目)。もう一度 airmiles を参照してみると、最初の要素が元に戻っています (25〜31行目)。これはデフォルト名前空間から airmiles が消えたことにより、また datasets パッケージの airmiles を参照するようになったためです。<br />
<br />
<h2>
データを編集する</h2>
edit() はテキストエディタを呼び出す関数ですが、データフレームや行列に対して呼び出されたときは、スプレッドシートが起動するようになっています。<br />
例えば Windows の場合は、以下のようになります。<br />
<pre class="prettyprint lang-r">> HousePrice4 <- edit(HousePrice)
</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-U5XaOZKJHyo/Xs1uAvMx9dI/AAAAAAAAANs/MFGeAjvUzwwJtJ0-d5h1MaF6lMxo54x5QCPcBGAYYCw/s1600/houseprice1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="166" data-original-width="387" src="https://1.bp.blogspot.com/-U5XaOZKJHyo/Xs1uAvMx9dI/AAAAAAAAANs/MFGeAjvUzwwJtJ0-d5h1MaF6lMxo54x5QCPcBGAYYCw/s1600/houseprice1.png" /></a></div>
<br />
edit() にデータフレームを渡すと、このようにスプレッドシートが起動します。このスプレッドシート上で Cent.heat[1] を yes に変更して、ウィンドウを閉じてみましょう。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-gtEI1XJLZq8/Xs1uAlJpLmI/AAAAAAAAANw/7ROFvHAI0jsXFwJdpB0S1EokUI-ziYyJQCLcBGAsYHQ/s1600/houseprice2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="166" data-original-width="387" src="https://1.bp.blogspot.com/-gtEI1XJLZq8/Xs1uAlJpLmI/AAAAAAAAANw/7ROFvHAI0jsXFwJdpB0S1EokUI-ziYyJQCLcBGAsYHQ/s1600/houseprice2.png" /></a></div>
<br />
edit() の結果を HousePrice4 に付値するようにしていたので、HousePrice4 を参照してみます。<br />
<pre class="prettyprint">> HousePrice4
Price Floor Area Rooms Age Cent.heat
01 52.00 111 830 5 6.2 yes
02 54.75 128 710 5 7.5 no
03 57.50 101 1000 5 4.2 no
04 57.50 131 690 6 8.8 no
05 59.75 93 900 5 1.9 yes</pre>
<br />
確かに変更されていることが確認できました。<br />
例えば、DF <- edit(data.frame()) のようにすれば、新しいデータフレームをスプレッドシートから定義することもできます。<br />
<br />管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-66189102854953061842020-05-24T08:20:00.000+09:002020-05-24T08:20:01.938+09:00アーキテクチャ設計のポイント業務要件を実現するだけであれば、アーキテクチャ設計は不要です。ではアーキテクチャ設計では何をすべきなのか? 私がアーキテクチャ設計を行うときのポイントを備忘録的にまとめました。<br />
<br />
<br />
<a name='more'></a><br />
<h2>
アーキテクチャ設計とは</h2>
アーキテクチャ設計は、システムのコンポーネント (構成要素) と、コンポーネント間のインターフェース (連携方法) を決めることに他なりません。アーキテクチャ設計により、システムはコンポーネントの集合体として再定義されます。コンポーネントには責務が与えられ、他コンポーネントとのインターフェースが指定されます。<br />
コンポーネントは、この責務とインターフェースを満たすように設計・実装されなければなりません。そのために各コンポーネントはまたそれぞれのアーキテクチャを持つことになります。アーキテクチャ設計は、それが自明、既知、実装済みなレベルに到達するまで、ブレイクダウンされていきます。<br />
<br />
<h2>
アーキテクチャのレベル</h2>
アーキテクチャのレベルは例えば次のような形になります。厳密なものではなく、状況に応じて適切なレベル感で記述します。<br />
<br />
<ul>
<li>エンタープライズレベル</li>
<ul>
<li>コンポーネント: 部門、部署、顧客、取引先 等</li>
<li>インターフェース: 手続、連絡手段、様式 等</li>
</ul>
<li>インターシステムレベル</li>
<ul>
<li>コンポーネント: 社内システム、外部システム 等</li>
<li>インターフェース: 接続/通信方式、介在業務 等</li>
</ul>
<li>インターサービスレベル:</li>
<ul>
<li>コンポーネント: サービス</li>
<li>インターフェース: プロトコル、ファイルフォーマット 等</li>
</ul>
<li>フレームワークレベル:</li>
<ul>
<li>コンポーネント: フレームワークコンポーネント、ライブラリ 等</li>
<li>インターフェース: 呼び出し、コールバック、キューイング、メッセージング、バッチ 等</li>
</ul>
<li>コードレベル:</li>
<ul>
<li>コンポーネント: メソッド、関数 等</li>
<li>インターフェース: 引数、戻り値、データ構造 等</li>
</ul>
</ul>
<br />
<br />
<h2>
アーキテクチャ設計と非機能要件</h2>
冒頭でも述べましたが、極端なことを言えば、業務要件の実現にややこしいアーキテクチャ設計は必要ありません。極端な話、多くの場合はコンピュータシステムが無くても、人間が紙とペンでがんばれば、出来ないことはないものです。しかし実際には、コンピュータシステムが必要であったり、既にシステム化されていたとしても、新しいシステムに変更する必要性があります。それは、業務要件とは違うビジネス上の要請、<b>非機能要件</b>のためです。非機能要件を実現するためにコンピュータシステムが必要で、そのためにはアーキテクチャ設計が必要なのです。<br />
非機能要件は、明確化するのがとても難しい領域です。ある程度高い水準で実現しなければならない一方、少し実現レベルを上げるだけで、<b>コストや期間が一気に増大する</b>ことが多いので、常に適正水準での妥協が必要です。<br />
非機能要件の観点は、インターサービスレベルから上の粒度であれば、<a href="https://www.ipa.go.jp/sec/softwareengineering/std/ent03-b.html">IPA 非機能要求グレード</a>などを参照すると良いと思います。例えばIPAでは以下の6カテゴリーに分類しています。<br />
<br />
<ul>
<li>可用性</li>
<li>性能・拡張性</li>
<li>運用・保守性</li>
<li>移行性</li>
<li>セキュリティ</li>
<li>システム環境・エコロジー</li>
</ul>
<br />
これら全てに対して、コスト・スケジュールとの兼ね合いを考えて実現レベルを決める必要があります。フレームワークレベル以下の粒度については、上位で決めた非機能要件を満たせるであろうように考慮して、決めていく必要があります。<br />
<br />
<h3>
特に重要なことは、未知の変更に対する弾力性・包容力</h3>
私は、アーキテクチャ設計において特に重要なことは、未知の変更に対する弾力性・包容力だと考えています。これは、コンポーネントが適度に適切に分割されていることと、コンポーネントが入れ替え可能であるように責務とインターフェースが決められていることに帰着します。<br />
<br />
<h2>
アーキテクチャ設計の観点</h2>
<div>
私がアーキテクチャ設計を行う際、以下のような観点を意識しています。</div>
<br />
<ul>
<li>非機能要件と、業務要件、コスト、スケジュールの間でコンフリクトが起きそうな事項は、重点的に必要性を確認する</li>
<li>非機能要件に以下のようなキーワードが出てきたら要注意</li>
<ul>
<li>リアルタイムで</li>
<li>直ちに</li>
<li>早く</li>
<li>同時に</li>
<li>必ず</li>
<li>確実に</li>
<li>絶対に</li>
<li>100%</li>
<li>永久に</li>
<li>無停止で</li>
<li>最優先で</li>
<li>全てのブラウザで</li>
<li>常に最新の環境で</li>
<li>BCP (事業継続計画)</li>
<li>DR (災害対策)</li>
</ul>
<li>キャパシティの余裕率は増強リードタイムに依存する</li>
<ul>
<li>高い余裕率は当然コストにはねる</li>
</ul>
<li>セキュリティ、コンプライアンスは関連部署のチェックを早い段階で入れる</li>
<ul>
<li>代替案の無い変更指示を受けることが多いので、後回しにするとコスト、スケジュールへの影響が出やすい</li>
</ul>
<li>保守性は「テスト可能性」と言い換えると分かりやすい</li>
<ul>
<li>テストを分割する境界の構築、ログ、リソース監視、エラー監視、検証環境の有無、……</li>
</ul>
<li>妥協したこと (=残存させたリスク) と、その理由はちゃんと残す</li>
<ul>
<li>複数案作って各案について全部残しておく</li>
</ul>
<li>全体的に「未知の変更に対する弾力性・包容力」のチェックを行うこと</li>
</ul>
管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-26939650395740173272020-05-23T05:31:00.002+09:002020-05-23T05:31:52.538+09:00時系列分析 (5) - 状態空間モデルAR, MA, ARMA, ARIMA, VAR, VARMA, VARIMA, ... などの時系列モデルは、観測値を直接モデル化するものでした。今回紹介する状態空間モデル (State-space model) の場合は、<b>状態</b>の時系列変化と、その状態から<b>観測される値</b>とに分けてモデル化する手法です。AR モデル、ARMA モデルなども状態空間表現を取ることができますし、状態空間モデルではより多様なシステムが記述できます。例えば、時変的な回帰係数を持つ回帰モデルなどが表現できます。<br />
また、状態空間モデルではカルマンフィルタという強力なアルゴリズムによって、条件付き同時分布が効率よく計算できるため、ARMA のパラメータ推定を行うに際して、ARMA の状態空間表現に対してのカルマンフィルタによる推定が良く行われています。<br />
<br />
<br />
<a name='more'></a><br />
<h2>
ガウス型状態空間モデル</h2>
<h3>
ガウス型状態空間モデルの定義</h3>
$Y_t$を$l$変量の時系列としたとき、これを以下の2つの方程式によって表現するモデルのことを状態空間モデルと呼びます。$$<br />
\begin{eqnarray*}<br />
X_t &=& F_tX_{t-1} + G_tv_t &\cdots システム方程式 \\<br />
Y_t &=& H_tX_t + w_t &\cdots 観測方程式<br />
\end{eqnarray*}<br />
$$<b>システム方程式は</b>、システムの<b>状態</b> ($X_t$) が1つ前の時点 ($X_{t-1}$) からどのように更新されるかを示しており、<b>観測方程式</b>は、システムの状態 ($X_t$) から<b>観測値</b> ($Y_t$) がどのように取り出されるかを示しています。<br />
$X_t$は<b>状態ベクトル</b>と呼ばれています。ここでは$k$変量のベクトルとします ($k$は状態空間モデルのパラメータになります)。$F_t$は<b>回帰係数行列</b>で$k$行$k$列の行列です。システムの過去の状態が現在にどのように影響を与えるかを示しています。$v_t$は<b>システムノイズ</b>と呼ばれ、システムの状態の変化の際に入り込む攪乱を表します。$m$変量の正規ホワイトノイズで、分散共分散行列$Q_t$に従うものとします。$G_t$はシステムノイズがシステムの状態に与える影響を表す$k$行$m$列の行列です。<br />
$H_t$はシステムの状態が観測値にどのように射影されるかを表す行列で、$l$行$k$列の行列です。$w_t$は<b>観測ノイズ</b>と呼ばれ、システムの状態を観測する際に入り込む攪乱を表します。$l$変量の正規ホワイトノイズで、分散共分散行列$R_t$に従うものとします。<br />
上記のような条件を置いた状態空間モデルは、ガウス型状態空間モデルとも呼ばれています。<br />
なお、状態空間モデルのシステム方程式、観測方程式は、真のモデルに対して一意ではありません。<br />
<br />
<h3>
AR モデルの状態空間表現</h3>
状態空間モデルの例として、$AR(p)$モデルの状態空間表現を見てみましょう。<br />
以下の式で定義される$y_t$は、$AR(p)$モデルです。$$<br />
y_t = \phi_1y_{t-1} + \phi_2y_{t-2} + \ldots + \phi_py_{t-p} + \epsilon_t, \epsilon_t \sim iidN(0, \sigma^2)<br />
$$定数項は$0$とし、攪乱項の分布は正規ホワイトノイズとします。$x_t = (y_t, y_{t-1}, \ldots, y_{t-p+1})'$とおくと、以下の関係式が成り立ちます。$$<br />
x_t = Fx_{t-1} + G\epsilon_t<br />
$$ただし、$$<br />
\begin{eqnarray*}<br />
F &=& \left(\begin{array}{ccccc}<br />
\phi_1 & \phi_2 & \ldots & \phi_p-1 & \phi_p \\<br />
1 & 0 & \ldots & 0 & 0 \\<br />
0 & 1 & \ldots & 0 & 0 \\<br />
\vdots & \vdots & \ddots & \vdots \\<br />
0 & 0 & \ldots & 1 & 0<br />
\end{array}\right) (p 行 p 列の行列) \\<br />
G &=& \left(\begin{array}{c}<br />
1 \\<br />
0 \\<br />
0 \\<br />
\vdots \\<br />
0<br />
\end{array}<br />
\right) (p 次元のベクトル)<br />
\end{eqnarray*}<br />
$$です。これが状態空間モデルにおけるシステム方程式にあたります。システムノイズの分散共分散行列$Q$は$Q = \sigma^2$になります。<br />
そして、$H = (1, 0, ..., 0)' (p 次元のベクトル)$と置けば、<br />
$y_t = HX_t$<br />
となり、これが観測方程式にあたります。観測ノイズの分散共分散行列$R$は$R = 0$となります。$F, G, H, Q, R$は、いずれも時不変なので、時不変の状態空間モデルです。<br />
<br />
<h3>
状態の推定</h3>
時系列データ$Y_t = (Y_1, Y_2, \ldots, Y_T)$に基づいて状態空間モデルを構築するにあたっては、状態$X_{t'} = (X_1, X_2, \ldots, X_T)$を推定する必要があります。また状態空間モデルから将来の観測値を予測するには、状態$X_{t'} = (X_{T+1}, X_{T+2}, \ldots)$の推定も必要です。<br />
$X_{t'} (1 \le t' < T)$の推定を<b>平滑化</b>、$X_{t'} (t' = T)$の推定を<b>フィルタ</b>、$X_{t'} (t' > T)$の推定を<b>予測</b>と呼んでいます。<br />
<br />
状態の推定は、$Y_t$が与えられた条件下での$X_{t'}$の条件付き分布$p(X_{t'}|Y_t)$を求める問題となります。状態空間モデルでは、ノイズが正規分布であるため、$p(X_{t'}|Y_t)$の分布も正規分布となります。つまり、$p(X_{t'}|Y_t)$の平均と分散共分散行列を求めれば、$X_{t'}$が推定できたことになります。<br />
<br />
<h3>
カルマンフィルタ</h3>
一般に$Y_t$が与えられた条件下での$X_{t'}$の条件付き同時分布を求める問題は、計算量が膨大になります。しかし、状態空間モデルにおいては逐次的に計算することが可能です。このアルゴリズムをカルマンフィルタと呼んでいます。<br />
$X_{t'}$の条件付き分布を求めるには、前述のようにその平均と分散共分散行列を求めれば良いことが分かっています。条件付き平均、条件付き分散共分散行列をそれぞれ以下のように表記することにします。$$<br />
\begin{eqnarray*}<br />
X_{t'|t} &=& E(X_{t'}|Y_t) \\<br />
V_{t'|t} &=& E((X_{t'} - X_{t'|t})(X_{t'} - X_{t'|t})')<br />
\end{eqnarray*}<br />
$$カルマンフィルタは、1期先予測 ($t = T - 1$の下での$t' = T$における$X_{t'}$の推定) と、フィルタ ($t = T$の下での$t' = T$における$X_{t'}$の推定) の2つの計算から成り立っています。<br />
1期先予測の計算は、以下のようになります。$$<br />
\begin{eqnarray*}<br />
X_{t'|t'-1} &=& F_{t'}X_{t'-1|t'-1} \\<br />
V_{t'|t'-1} &=& F_{t'}V_{t'-1|t'-1}F_{t'}' + G_{t'}Q_{t'}G_{t'}'<br />
\end{eqnarray*}<br />
$$1期先予測の平均は、現時点のフィルタに推移行列を掛けたものです。システムノイズの平均が$0$ですから、直感的にも分かります。分散共分散の第1項は推移行列の影響で、第2項はシステムノイズの影響を表しています。<br />
フィルタの計算は、以下のようになります。$$<br />
\begin{eqnarray*}<br />
K_{t'} &=& V_{t'|t'-1}H_{t'}'(H_{t'}V_{t'|t'-1}H_{t'}' + R_{t'})^{-1} \\<br />
X_{t'|t'} &=& X_{t'|t'-1} + K_{t'}(Y_t - H_{t'}X_{t'|t'-1}) = K_{t'}Y_t + (I - K_{t'}H_{t'})X_{t'|t'-1} \\<br />
V_{t'|t'} &=& V_{t'|t'-1} - K_{t'}H_{t'}V_{t'|t'-1} = (I - K_{t'}H_{t'})V_{t'|t'-1}<br />
\end{eqnarray*}<br />
$$ $K_{t'}$はカルマンゲインと呼ばれています。$Y_t - H_{t'}X_{t'|t'-1}$は$Y_t$の観測誤差、$H_{t'}V_{t'|t'-1}H_{t'}' + R_{t'}$はその分散共分散行列を示しており、カルマンゲインは観測誤差を状態の推定誤差に繰り込む働きをしています。<br />
フィルタの平均の後者の等式 ($X_{t'|t'} = K_{t'}Y_t + (I - K_{t'}H_{t'})X_{t'|t'-1}$) を見ると、フィルタが新しい観測値と1時点前までの情報による1期先予測との加重平均であると見ることもできます。<br />
分散共分散の前者の形 ($V_{t'|t'} = V_{t'|t'-1} - K_{t'}H_{t'}V_{t'|t'-1}$) を見ると、新しい観測値によって状態推定の精度が改善される、と解釈することもできます。<br />
<br />
<h2>
非ガウス型状態空間モデル</h2>
ガウス型状態空間モデルでは、ノイズが正規分布に従うため、平均値から離れた値を取る可能性が低く、状態や観測値の変位が滑らかになります。従って、観測値の変位が急激な場合は、あまり適合が良くありません。<br />
そこで、ノイズの分布を正規分布と仮定せず、一般のホワイトノイズであるとの仮定に緩めたものが、非ガウス型状態空間モデルです。ノイズにファットテールな分布を持ってくれば、観測値の急激な変化もモデルに織り込むことができる可能性が出てきます。<br />
<br />
<h3>
非ガウス型状態空間モデルの定義</h3>
定式化すると、以下のようになります。$$<br />
\begin{eqnarray*}<br />
X_t &=& F_tX_{t-1} + G_tv_t &\cdots システム方程式 \\<br />
Y_t &=& H_tX_t + w_t &\cdots 観測方程式<br />
\end{eqnarray*}<br />
$$方程式の形状はガウス型と同じで、各係数行列の条件も同じですが、$v_t, w_t$はそれぞれ、密度関数$q(v), r(w)$に従うホワイトノイズであるとします。<br />
<br />
<h3>
状態の推定</h3>
推定は、カルマンフィルタの拡張によって行うことができます。ガウス型では条件付き同時分布が正規分布であるという条件を置くことができましたが、非ガウス型ではそうとは限らないため、条件付き同時分布を直接求める必要があります。<br />
1期先予測は、以下のようになります。$$<br />
p(X_{t'}|Y_{t-1}) = \int p(X_{t'}|X_{t'-1}) p(X_{t'-1}|Y_1, \ldots, Y_{t-1}) dX_{t'-1}<br />
$$フィルタは、以下のようになります。$$<br />
p(X_{t'}|Y_t) = \frac{p(Y_{t'}|X_{t'}) p(X_{t'}|Y_1, \ldots, Y_{t-1})}{p(Y_{t'}|Y_1, \ldots, Y_{t-1})}<br />
$$<br />
<br />
<h2>
一般化状態空間モデル</h2>
非ガウス型状態空間モデルでは、状態の推移の線形性を仮定していましたが、その仮定を緩め、非線形な推移も許したものが一般化状態空間モデルです。<br />
例えば、確率的ボラティリティモデルのパラメータ推定を行うために、非線形変換を加えた上で、線形のガウス型状態空間モデルで近似して推定を行う方法がありますが、一般化状態空間モデルでは確率的ボラティリティモデルを近似無しに表現することができます。<br />
<br />
<h3>
一般化状態空間モデルの定義</h3>
一般化状態空間モデルは、以下のように定式化されます。$$<br />
\begin{eqnarray*}<br />
X_t &=& F(X_{t-1}, v_t) &\cdots システム方程式 \\<br />
Y_t &=& H(X_t, w_t) &\cdots 観測方程式<br />
\end{eqnarray*}<br />
$$ $v, w$は非ガウス型状態空間モデルと同様、密度関数$q(v), r(w)$に従うホワイトノイズであるとします。$F(X, v)$および$H(X, w)$は非線形の関数です。$H(X, w)$については、$Y = H(X, w)$に対して$w = G(Y, X)$を満たす$w$を一意に決定できる、$Y$について微分可能な関数$G(Y, X)$が存在することを要求します。また、初期状態$X_0$は密度関数$p_0(X)$に従うものとします。<br />
<br />
<h3>
モデルの推定</h3>
モデルの推定は解析的に行うことは一般にはできません。モンテカルロ・フィルタと呼ばれるアルゴリズムによって行います。モンテカルロ・フィルタでは、仮定した分布に従う「粒子」を乱数で生成し、解析的に表現できない同時分布をこの乱数に基づく経験分布関数で近似するという方法をとります。概要としては以下の通りです。<br />
<br />
<ol>
<li>用いる「粒子」の数を$m$とし、以下、$j = 1, \ldots, m$とする</li>
<li>$p_0(X)$に従う$m$個の乱数ベクトル$f_0^{(j)}$を生成する (初期状態を乱数で設定)</li>
<li>$t = 1, \ldots, T$について、以下の手順を実行する</li>
<ol>
<li>$q(v)$に従う$m$個の乱数ベクトル$v_t^{(j)}$を生成する (システムノイズを乱数で生成)</li>
<li>各$j$について、$p_n^{(j)} = F(f_{n-1}^{(j)}, v_n^{(j)})$を計算する (新しい状態の候補値の算出)</li>
<li>各$j$について、$\alpha_n^{(j)} = r(G(Y_n, p_n^{(j)})) |\partial G/\partial Y_n|$を計算する (状態$p_n^{(j)}$が観測値$Y_n$となる確率の算出)</li>
<li>$p_n^{(j)}$の中から$\alpha_n^{(j)}$に比例する確率で$m$回サンプリングを行い、その値を$f_n^{(1)}, \ldots, f_n^{(m)}$とする (新しい状態を決定)</li>
</ol>
</ol>
<br />
サンプリングは、 ランダムサンプリングでも良いのですが、計算量を減らすため階層化リサンプリング法といったものを用いることもできます。<br />
<div>
<br /></div>
管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-4636283455365120022020-05-23T01:33:00.004+09:002020-05-23T01:33:57.833+09:00時系列分析 (4) - 多変量の時系列分析複数の変量が相互作用を持って発展していく形の時系列は珍しくありません。こういったデータは、各変量を個別にモデル化しようとしてもうまくいきません。このような場合は、ベクトル自己回帰 (Vector Autoregressive / VAR) モデルが有効な場合があります。VAR モデルは、AR モデルを多変量に拡張したものです。今回は VAR モデルについて簡単に触れたいと思います。<br />
<br />
<br />
<a name='more'></a><br />
<h2>
ベクトル過程</h2>
VAR モデルの構築にあたっては、観測データをベクトルの時系列 (ベクトル過程) とみなすことから始めます。<br />
$n$個の変量からなる観測データの時点$t$における値を$n$行$1$列のベクトル$Y_t = (y_{1,t}, y_{2,t}, \ldots, y_{n,t})'$と考え、観測データ全体をこのベクトルからなるベクトル過程として捉えます。その上で、スカラー観測値を対象に構築された各種概念をベクトル観測値に拡張して適用します。<br />
<br />
<h3>
期待値、自己共分散のベクトル過程への拡張</h3>
期待値の拡張としての期待値ベクトル$E(Y_t)$は、各変量の期待値からなるベクトル、すなわち$E(Y_t) = (E(y_{1,t}), E(y_{2,t}), \ldots, E(y_{n,t}))'$と定めます。<br />
自己共分散の拡張は自己共分散行列となります。$k$次-自己共分散行列$Cov(Y_t, Y_{t-k})$は、$Cov(Y_t, Y_{t-k}) = (Cov(y_{i,t}, y_{j,t-k}))_{ij}$となります。これは$n$行$n$列の正方行列で、対角成分は各変量の自己共分散に等しくなっています。また、$0$次-自己共分散行列は、各変量を通常の確率変数として見たときの分散共分散行列にあたります。$k$の関数として見て、自己共分散関数と呼ぶこともできます。<br />
<br />
<h3>
ベクトル過程の弱定常性</h3>
期待値と自己共分散が定義できたので、これでベクトル過程の弱定常性が定義できます。ベクトル過程$Y_t$の期待値ベクトルと自己共分散関数が時点$t$に依存しないとき、$Y_t$は弱定常ベクトル過程と呼ばれます。弱定常過程においては、期待値ベクトルも自己共分散関数もともに$t$に依存しないので、以後それぞれ$\mu, \Gamma_k$で表すことにします。<br />
<br />
<h3>
自己相関行列</h3>
自己共分散と同様、自己共分散行列も変量の大きさに依存していますので、関係の強弱を比較しづらいものとなっています。そのため、自己相関のように基準化する必要があります。ベクトル過程に対する自己共分散行列を基準化したものを自己相関行列と呼びます。自己相関行列$\rho_k$は$\rho_k = (Corr(y_{i,t}, y_{j,t-k}))_{ij}$と定義されます。<br />
$\rho_k$は、以下のようにも書けます。$$<br />
\rho_k = D^{-\frac{1}{2}}\Gamma_kD^{-\frac{1}{2}}<br />
$$ただし、$D = diag(Var(y_1), \ldots, Var(y_n)))$です。$diag()$は対角行列の意味です。<br />
<br />
<h3>
ベクトルホワイトノイズ</h3>
ホワイトノイズの定義もベクトルに拡張しておきましょう。ベクトル過程 $\epsilon_t$が、$$<br />
\begin{eqnarray*}<br />
E(\epsilon_t) &=& 0 \\<br />
E(\epsilon_t\epsilon_{t-k}') &=& \begin{cases}<br />
\Sigma &(k = 0) (\Sigmaは正定値で\epsilon_tの分散共分散行列にあたる)\\<br />
0 &(k \ne 0)<br />
\end{cases}<br />
\end{eqnarray*}<br />
$$であるとき、$\epsilon_t$はベクトルホワイトノイズであると言い、$\epsilon_t \sim WN(\Sigma)$と表記します。$\Sigma$は対角行列でなくともかまいません。このことは、異なる時点間ではどの変量においても他の変量とも自分自身とも相関を持っていてはならないが、同一時点内において変量間に相関があることは許容されるということを意味します。<br />
<br />
<h2>
VAR モデル</h2>
VAR モデルは、AR モデルをベクトルに拡張したものです。$VAR(p)$は$p$時点まで遡った値で回帰したモデルということになります。数式で書くと以下のようになります。$$<br />
Y_t = c + \Phi_1Y_{t-1} + \ldots + \Phi_pY_{t-p} + \epsilon_t, \epsilon_t \sim WN(\Sigma)<br />
$$ $c$は$n$行$1$列の定数ベクトル、$\Phi_1$から$\Phi_p$は$n$行$n$列の行列となります。<br />
$n$変量の$VAR(p)$モデルは、実質的には$n$本の回帰式からなり、各回帰式が$np$個の回帰係数と$1$個の定数を持っています。また、攪乱項の分散共分散行列 (\Sigma) の中に$\frac{n(n+1)}{2}$個のパラメータがあるので、トータルで$(p + \frac{1}{2}) n^2 + \frac{3n}{2}$という膨大な数のパラメータを持っていることになります。最もパラメータ数の少ない$2$変量$VAR(1)$でもパラメータ数は$9$個あります。<br />
<br />
<h3>
VAR モデルの定常条件</h3>
VAR モデルの定常性は、AR 特性方程式をベクトルに拡張したもので判定できます。すなわち、$|I_n - \Phi_1z - \ldots - \Phi_pz^p| = 0$の全ての解の絶対値が1より大きいことが定常条件です ($I_n$は$n$行$n$列の単位行列です)。<br />
<br />
<h2>
VMA モデルと VARMA モデル</h2>
MA モデル、ARMA モデルも VAR と同様にベクトルに拡張することができます。VAR を $VMA(\infty)$に書き直すこともできます。<br />
VARMA モデルは VAR よりもさらに膨大なパラメータ数になるので、あまり使われません。<br />
<br />
<h2>
VAR モデルの推定</h2>
VAR モデルにおいては、各回帰式を独立に最小二乗法で推定したものが最良の線形不偏推定量となることが分かっており、攪乱項が多変量正規分布に従う場合は最尤推定量とも一致します。次数は基本的には情報量規準に基づいて選択すれば良いのですが、パラメータ数が多いことから、あまり有意でないパラメータがオーバーフィッティングを起こし、小さすぎる次数のモデルが選択される恐れがあります。予め定性的に高い次数が必要とされることが分かっている場合は、それを考慮にいれて推定を行うべきです。<br />
<br />管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-53356714729359256392020-05-22T19:48:00.000+09:002020-05-23T00:21:37.223+09:00時系列分析 (3) - モデルの推定と診断いよいよ、時系列データからモデルを作ることを考えます。モデルを推定し、そのモデルの善し悪しを診断する方法を紹介します。<br />
<br />
<a name='more'></a><h2>
モデルの推定</h2>
時系列モデルの推定とは、時系列データがどのような種類の確率過程に従うか、そしてそのパラメータは何かを推定することです。モデルが推定できれば、未来の時点の値の予測に用いることができます。<br />
モデルを推定する方法としては、最小二乗法 (Ordinary least squares/OLS) と、最尤推定 (Maximum likelihood estimation/MLE) が良く知られています。最小二乗法はシンプルかつ良い推定量を導出できるのでよく用いられますが、複雑なモデルでは適用できないことがあり、その場合は最尤推定を用いることになります。<br />
<br />
<h2>
最小二乗法</h2>
最小二乗法は、モデルによって説明できない項 (誤差項) を最小にするようにパラメータを推定する方法です。<br />
例えば、ある過程$y_t$を$AR(1)$モデル$y_t = c + \phi y_{t-1} + \epsilon_t, \epsilon_t \sim WN(\sigma^2)$に当てはめることを考えましょう。$c'$と$\phi'$がそれぞれ$c$と$\phi$の推定値だとすると、このモデルで説明できない値$e_t$は$e_t = y_t - c' - \phi'y_{t-1}$となります。これを残差 (residuals) と呼びます。この残差の絶対値を小さくすれば、良いモデルに近づくことになります。最小二乗法では、残差平方和 (Sum of squared residuals/SSR) を最小にするようなパラメータを選択します。このモデルでは、残差平方和は$SSR = \sum_{t=2}^T(y_t - c' - \phi'y_{t-1})^2$となります。この$SSR$を最小化するような$c', \phi'$ (これを OLS 推定量と呼び、ここでは$c'', \phi''$ で表すことにする) を求めれば良いことになります。そのため、$c', \phi'$それぞれで偏微分します。<br />
$$<br />
\begin{eqnarray*}<br />
\frac{\partial SSR}{\partial c'} &=& -2\sum_{t=2}^T (y_t - c' - \phi'y_{t-1}) \\<br />
\frac{\partial SSR}{\partial \phi'} &=& -2\sum_{t=2}^T y_{t-1}(y_t - c' - \phi'y_{t-1})<br />
\end{eqnarray*}<br />
$$<br />
この2式が$0$になるような$c', \phi'$が$c'', \phi''$となるので、$c'', \phi''$は以下の2式を満たします。<br />
$$<br />
\begin{eqnarray*}<br />
\sum_{t=2}^T (y_t - c'' - \phi''y_{t-1}) = 0 \\<br />
\sum_{t=2}^T y_{t-1}(y_t - c'' - \phi''y_{t-1}) = 0<br />
\end{eqnarray*}<br />
$$<br />
これを正規方程式と呼び、解くと以下のようになります。ただし、$\bar{y}_{a,b} = \frac{1}{b-a+1}\sum_{t=a}^b y_t$ ($y_a \sim y_b$の平均)です。<br />
$$<br />
\begin{eqnarray*}<br />
c'' &=& \bar{y}_{2,T} - \phi''\bar{y}_{1,T-1} \\<br />
φ'' &=& \frac{\sum_{t=2}^T (y_t - \bar{y}_{2,T})(y_{t-1} - \bar{y}_{1,T-1})}{\sum_{t=2}^T (y_{t-1} - \bar{y}_{1,T-1})^2}<br />
\end{eqnarray*}<br />
$$<br />
誤差項の分散$\sigma^2$の OLS 推定量$s^2$は、OLS 残差の不偏標本分散なので、以下のようになります。<br />
$$<br />
s^2 = \frac{1}{T - 2}\sum_{t=2}^T (y_t - c'' - \phi''y_{t-1})^2<br />
$$<br />
OLS 推定量は、不偏推定量 (平均的に真の値をとらえる)、一致推定量 (標本数が大きくなれば真の値に限りなく近づく) の性質を持ちます。また、線形不偏推定量の中で最小の分散共分散行列を持ち、$\epsilon_t \sim iidN(\sigma^2)$であれば任意の不偏推定量の中で最小の分散共分散行列を持ちます。つまり、OLS はこういった条件の下で<b>最適性</b>を持っています。OLS 推定量を基準化すると、正規分布に漸近するため、仮説検定ではこの性質も活用します。<br />
<br />
<h2>
最尤法</h2>
最尤法は、モデルが観測値を実現しやすくなるようにパラメータを選択する方法です。もう少し具体的に言うと、観測値が実現される確率を<b>尤度関数</b> (または単に尤度) というパラメータの関数として見て、この尤度関数を最大化するという方法です。尤度関数を最大化するパラメータを<b>最尤推定値</b>と呼びます。最尤推定値は解析的に求めることができる場合もありますが、解析解が無く数値的に求めなければならない場合もあります。解析的に求める場合は尤度関数を微分することになりますが、尤度関数の最大化と尤度関数の対数を取ったもの (<b>対数尤度</b>と言う) の最大化は同じ結果となるので、微分しやすい対数尤度を用いることが多いです。<br />
<br />
<h3>
例1 簡単な確率の推定</h3>
例えば簡単な例として、「サイコロを10回振ったら、1が2回出た」という観測値があったとします。ここから、1が出る確率$p$を最尤法で推定してみましょう。<br />
<br />
尤度関数$L(p)$は、「1が出る確率が$p$のときに、10回中2回1が出る確率」ですから、$L(p) = {}_{10}\rm{C}_2p^2(1-p)^8$となります。このままでは微分しづらいので対数をとると、$\log L(p) = \log({}_{10}\rm{C}_2p^2(1-p)^8) = \log {}_{10}\rm{C}_2 + 2 \log p + 8 \log(1-p)$となります。<br />
微分すると、$$\frac{d}{dp} \log L(p) = \frac{2}{p} - \frac{8}{1-p} = \frac{2 - 10p}{p(1-p)}$$なので、これが$0$となるような$p$が最尤推定値となります。$2 - 10p = 0$すなわち$p = 0.2$となります。<br />
<br />
<h3>
例2 AR(1) のパラメータ推定</h3>
最小二乗法で例示したものと同様の$AR(1)$のパラメータ推定を最尤法で行ってみましょう。<br />
今回は$y_t = c + \phi y_{t-1} + \epsilon_t, \epsilon_t \sim iidN(0, \sigma^2)$とします。最小二乗法の時は$\epsilon_t \sim WN(\sigma^2)$としていましたが、$\epsilon_t$の分布が分からないと尤度関数が定まらないので、一般の$WN$ではなく、$iidN$としています。<br />
<br />
尤度関数は、同時密度を条件付き確率に分解する$P(X, Y) = \frac{P(X | Y)}{P(Y)}$を用いると、確率密度関数$f()$を使って$$L(c, \phi, \sigma^2) = f_{Y_T, Y_{T-1}, \ldots, Y_2|Y_1}(y_T, y_{T-1}, \ldots, y_2|y_1; c, \phi, \sigma^2) = \prod_{t=2}^T f_{Y_t|Y_{t-1}}(y_t|y_{t-1}; c, \phi, \sigma^2)$$と書けます。対数をとると、$$\log L(c, \phi, \sigma^2) = \sum_{t=2}^T \log f_{Y_t|Y_{t-1}}(y_t|y_{t-1}; c, \phi, \sigma^2)$$となります。ここで$\epsilon_t \sim iidN(0, \sigma^2)$を用いると、$y_t|y_{t-1} \sim N(c + \phi y_{t-1},\sigma^2)$より$$<br />
f_{Y_t|Y_{t-1}}(y_t|y_{t-1}; c, \phi, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}exp\left(-\frac{(y_t-(c+\phi y_{t-1}))^2}{2\sigma^2}\right)<br />
$$なので、対数尤度は$$<br />
\log L(c, \phi, \sigma^2) = -\frac{T}{2}log 2\pi - \frac{T}{2}\log\sigma^2 - \sum_{t=2}^T\left(\frac{(y_t-c-\phi y_{t-1})^2}{2\sigma^2}\right)<br />
$$となります。<br />
この対数尤度の最大化は、$$\sum_{t=2}^T (y_t - c - \phi y_{t-1})^2$$の最小化になりますが、これは最小二乗法で見た残差平方和そのものです。つまり、最尤推定量は OLS 推定量に一致することになります。<br />
<br />
なお、最尤推定量は不偏推定量とは限りません。例えば、$\sigma^2$の最尤推定量は、$\sigma^2 = \frac{1}{T}\sum_{t=2}^T(y_t - c - \phi y_{t-1})^2$ですが、これは OLS 推定量とは一致しません。OLS 推定量の方が不偏推定量で、最尤推定量は不偏推定量ではないのです。最尤推定量は一致推定量ではあります。また基準化すると正規分布に漸近します。最尤推定量は、一致推定量の中で漸近的に最小の分散共分散行列を持ちます。<br />
<br />
<h2>
偏自己相関</h2>
パラメータ推定以前の問題として、そもろもモデルとして MA, AR, ARMA のどれが適切であるかを選択する必要があります。これらのモデルの違いは、自己相関の組み込み方の違いですから、選択にあたってはまずは標本自己相関関数を見て判断することになります。ある次数$q$から先が$0$に近いようであれば、$MA(q)$で当てはめられる可能性が高いわけです。しかしそうではない場合、AR 過程とみるか、ARMA 過程とみるかは、標本自己相関関数からの判断は難しくなってきます。この判断にあたっては、<b>偏自己相関</b>を用いるのが適切です。<br />
$k$次-自己相関は、$y_t$と$y_{t-k}$の相関でしたが、$k$次-偏自己相関は、$y_t$と$y_{t-k}$の相関からその間にある$y_{t-1}, \ldots, y_{t-k+1}$の影響を除去したものです。自己相関関数と同様、$k$の関数として見て、偏自己相関関数とも言います。<br />
<br />
自己相関関数は、$AR(p), ARMA(p, q)$では減衰、$MA(q)$では$q$次を越えると$0$になる性質がありましたが、これに対して偏自己相関関数は、自己回帰項の影響が取り除かれることによって、$AR(p)$では$p$次を越えると$0$になる性質があります。$MA(q), ARMA(p, q)$では減衰します (MA は$AR(\infty)$で書き直せることからわかります)。<br />
<br />
従って、標本自己相関関数と標本偏自己相関関数を見ることによって、MA, AR, ARMA のいずれが適しているかを見定めることができます。MA, AR であれば切断箇所から次数も見積もることができます。ARMA の次数は自己相関、偏自己相関からは分かりません。またもちろん、標本自己相関関数、標本偏自己相関関数には真の値との間に誤差があることも忘れてはなりません。<br />
<br />
<h2>
情報量規準</h2>
あるモデル1つを仮定すれば、そのパラメータの推定量を最尤法などで求めることはできますが、そもそもどのモデルが最適なのかを決める術が必要です。その道具となるのが<b>情報量規準</b>です。<br />
情報量規準は、モデルの当てはまりの良さと、モデルの複雑さに対するペナルティから構成されています。情報量規準$IC$は以下のような一般形で表されます。$$<br />
IC = -2\log L(\Theta) + p(T)k<br />
$$<br />
$\Theta$は最尤推定値ベクトルで、$L(\Theta)$は最大尤度、$k$はパラメータ数、$p(T)$はペナルティ関数です。つまり第1項がモデルの当てはまりの良さを表現しており、第2項が複雑さに対するペナルティを表しています。$IC$の小さい方が、良いモデルであると考えます。<br />
<br />
<h3>
AICとSIC</h3>
ペナルティ関数に選び方によって$IC$の値は変わってきます。良く使われるのが赤池情報量規準 ($AIC$) と Schwarts 情報量規準 ($SIC$) です。$SIC$はベイズ情報量規準 ($BIC$) とも呼ばれています。$AIC$は$p(T) = 2$、$SIC$は$p(T) = \log T$としたもので、全体としては以下のようになります。$$<br />
\begin{eqnarray*}<br />
AIC &=& -2\log L(\Theta) + 2k \\<br />
SIC &=& -2\log L(\Theta) + k\log T<br />
\end{eqnarray*}<br />
$$<br />
$T \ge 8$では$SIC > AIC$です。通常は$T \ge 8$でしょうから、$SIC$は$AIC$よりもペナルティが大きいことになります。すなわち、$SIC$に基づく方がより小さいモデル (パラメータの少ないモデル) を選択する傾向があるということになります。また、$AIC$は$T$を大きくしても、パラメータの多すぎるモデルを選択してしまう確率が$0$になりませんが、$SIC$は$T$を大きくするほど最適なモデルを選択する確率が上がっていきます。ただし、$AIC$が過大なパラメータ数を持つモデルを選択しても、多すぎる部分のパラメータの最尤推定量は$0$に漸近するので、$AIC$も十分に有用です。<br />
<br />
<h2>
モデルの診断</h2>
最適とされるモデルが選択されたら (あるいは候補が絞り込まれたら)、そのモデルが適切かどうかの診断を行います。この診断で適切と判断されなければ、これまで候補に挙がったモデルの中に適切なモデルが無かったということになるので、おそらくはモデルの候補の範囲を広げて推定からやり直すことになるでしょう。<br />
<br />
例えば真のモデルが MA, AR, ARMA 等だとすると、$\epsilon_t \sim WN(\sigma^2)$ですから、$\epsilon_t$は自己相関を持たないはずです。そこで、最適とされたモデルにおける残差が自己相関を持たないことをかばん検定を用いて検定するという方法が使えます。ただし、パラメータ数の分だけ自由度が下がることを考慮し、$Q(m)$と$\chi^2(m-k)$ ($k$はパラメータ数) とで検定を行います。管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-89876289080154074302020-05-21T02:56:00.001+09:002020-05-28T09:55:32.179+09:00時系列分析 (2) - 自己相関のモデルある時系列データが、自己相関検定を経て自己相関があると分かったら、その自己相関のモデル化に取り組む価値があります。自己相関のモデル化にあたっては、移動平均過程 (MA過程) と自己回帰過程 (AR過程) という2つの過程が基本となります。この2つの過程とその組み合わせである自己回帰移動平均過程 (ARMA過程) について見ていきましょう。<br />
<br />
<a name='more'></a><h2>
自己相関のモデル化</h2>
前述の通り、自己相関のモデル化は、MA過程とAR過程という2つの過程が基礎となります。時系列データに$m$次-自己相関があることが想定されるとしましょう。このとき、MA 過程では、$y_t$と$y_{t-m}$のモデルに共通の成分を含めることで自己相関をモデル化します。共通の成分があるので相関が発生する、ということです。AR 過程では、$y_t$のモデルに$y_{t-m}$を含めることでモデル化します。$y_t$に$y_{t-m}$が含まれているので相関が発生する、ということです。<br />
<br />
<h2>
MA 過程 (移動平均過程)</h2>
MA 過程 (移動平均過程/Moving average process) は、のことである。$q$次-MA 過程$MA(q)$とは、平均$\mu$に攪乱項としてホワイトノイズ$\epsilon_t \sim WN(\sigma^2)$を加え、さらに$\epsilon_{t-1}$から$\epsilon_{t-q}$までの重み付け合計値 (重み$\theta_{t-1} \sim \theta_{t-q}$)を加えることにより、$q$次までの自己相関を組み込んだものです。ホワイトノイズを$\mu$だけ底上げし、自己相関項を組み込んで拡張したものと捉えることができます。<br />
定式化すると $y_t = \mu + \epsilon_t + \theta_1\epsilon_{t-1} + \theta_2\epsilon_{t-2} + \ldots + \theta_q\epsilon_{t-q}, \epsilon_t \sim WN(\sigma^2)$となります。<br />
<br />
<h3>
MA(1)</h3>
MA 過程の最も単純な形である、$MA(1)$ について詳しく見てみましょう。<br />
<br />
<h4>
MA(1)の構造</h4>
$MA(1)$は、前述の$MA(q)$の式について、$q=1$としたものですから、$y_t = \mu + \epsilon_t + \theta_1\epsilon_{t-1}, \epsilon_t \sim WN(\sigma^2)$と定義されます。$y_{t-1} = \mu + \epsilon_{t-1} + \theta_1\epsilon_{t-2}$ですから、$y_t$と$y_{t-1}$は共通項$\epsilon_{t-1}$を持っています。それによって1次-自己相関を生じています。<br />
<br />
<h4>
MA(1)の基本統計量</h4>
$MA(1)$の平均は、$E(y_t) = E(\mu + \epsilon_t + \theta_1\epsilon_{t-1}) = E(\mu) + E(\epsilon_t) + \theta_1E(\epsilon_{t-1}) = \mu$ となります ($\epsilon_t, \epsilon_{t-1}$はホワイトノイズなので平均$0$です)。<br />
$MA(1)$の分散$\gamma_0 = Var(y_t) = Cov(y_t, y_t)$ (分散は$0$次-自己共分散です) は、$\gamma_0 = Var(\mu + \epsilon_t + \theta_1\epsilon_{t-1}) = Var(\epsilon_t) + 2\theta_1Cov(\epsilon_t, \epsilon_{t-1}) + \theta_1^2Var(\epsilon_{t-1}) = \sigma^2 + 0 + \theta_1^2\sigma^2 = (1 + \theta_1^2)\sigma^2$ となります。攪乱項であるホワイトノイズの分散が$\sigma^2$ですから、それより$\theta_1^2\sigma^2$だけ大きくなっています。また、$\theta_1$が正の場合は$\theta_1\epsilon_{t-1}$と$\epsilon_{t-1}$が同じ符号であるため、$y_t$と$y_{t-1}$が同じ方向に動く傾向が出ます (<b>正の自己相関</b>)。負であれば、逆方向に動く傾向が出ることになります (<b>負の自己相関</b>)。また$\theta_1$の絶対値が大きいほど自己相関は強くなります。グラフ的に言うと、正の自己相関が強いほどグラフは滑らかになり、負の自己相関が強いほどグラフは尖ります。<br />
<br />
$MA(1)$の自己相関係数を求めて、どのように自己相関が折り込まれているか見てみましょう。$1$次-自己共分散$\gamma_1$は、$$<br />
\begin{eqnarray*}<br />
\gamma_1 &=& Cov(y_t, y_{t-1}) = Cov(\mu + \epsilon_t + \theta_1\epsilon_{t-1}, \mu + \epsilon_{t-1} + \theta_1\epsilon_{t-2}) \\<br />
&=& Cov(\epsilon_t, \epsilon_{t-1}) + Cov(\epsilon_t, \theta_1\epsilon_{t-2}) + Cov(\theta_1\epsilon_{t-1}, \epsilon_{t-1}) + Cov(\theta_1\epsilon_{t-1}, \theta_1\epsilon_{t-2}) \\<br />
&=& 0 + 0 + \theta_1Cov(\epsilon_{t-1}, \epsilon_{t-1}) + 0 \\<br />
&=& \theta_1\sigma^2<br />
\end{eqnarray*}<br />
$$となります。<br />
<br />
$1$次-自己相関係数$\rho_1$は、$\rho_1 = \frac{\gamma_1}{\gamma_0} = \frac{\theta_1}{1+\theta_1^2}$となります。$|\theta_1| = 1$のとき、$|\rho_1|$は最大値$1/2$をとります。つまり、$1$次-自己相関係数の絶対値が$1/2$を越えるような過程は$MA(1)$ではモデル化できないことがわかります。また、$2$次より大きい自己相関係数は全て$0$となります。<br />
上記から、$MA(1)$は平均も自己共分散も時点に依存しないことがわかります。したがって$MA(1)$は弱定常過程です。<br />
<br />
<h3>
MA(q)</h3>
$MA(q)$に一般化すると、導出は省略しますが、平均は$\mu$、$k$次-自己相関係数は<br />
$$<br />
\rho_k = \begin{cases}<br />
\frac{\theta_k + \theta_1\theta_k-1 + ... + \theta_q-k\theta_q}{1 + \theta_1^2 + \theta_2^2 + \ldots + \theta_q^2} & (1 \le k \le q) \\<br />
0 & (k \ge q + 1)<br />
\end{cases}<br />
$$<br />
となります。$MA(q)$もやはり弱定常過程です。<br />
<br />
<h3>
MA過程の反転可能性</h3>
MA 過程では、任意の MA 過程に対して、それと同一の平均、同一の自己相関関数を持つ MA 過程が複数 ($2q$個) 存在するという問題があります。大雑把に言えば、移動平均項の係数とホワイトノイズの分散の組み合わせ方が複数存在します。そこで、扱いやすい MA 過程を選ぶ基準として、<b>反転可能性</b>というものがあります。反転可能とは、「$AR(\infty)$に書き換え可能」という意味です。反転可能である場合、$y_t$を過去の$y$から自己回帰的に予測したときに、$\epsilon_t$の予測誤差を生じると解釈できます。このときの$\epsilon_t$を本源的攪乱項と呼びます。書き換え可能な$2q$個の同一平均、同一自己相関関数の$MA(q)$過程に対して、反転可能なものは1つしかありません。<br />
<br />
$MA(q)$の反転可能条件は、$1 + \theta_1z + \theta_2z^2 + \ldots + \theta_qz^q = 0$(MA 特性方程式) の全ての解の絶対値が$1$より大きいことです。$MA(1)$であれば、$|\theta| < 1$が反転可能条件です。<br />
<br />
<h2>
AR 過程</h2>
AR 過程とは、自己回帰過程 (Autoregressive process) のことです。$p$次 AR 過程$(AR(p))$とは、平均$\mu$に攪乱項$\epsilon_t \sim WN(\sigma^2)$を加え、さらに$y_{t-1}$から$y_{t-p}$までの重み付け合計値を加えることにより、自己相関を組み込んだものです。後述しますが、$p$次を越える自己相関が生じます。$y_t = c + \phi_1y_{t-1} + \phi_2y_{t-2} + ... + \phi_py_{t-p} + \epsilon_t, \epsilon_t ~ WN(\sigma^2)$と書けます。<br />
<br />
<h3>
AR(1)</h3>
AR 過程の性質を知るため、やはり$AR(1)$について詳しく見てみましょう。<br />
<br />
<h4>
AR(1)の構造</h4>
$AR(1)$は、$y_t = c + \phi_1y_{t-1} + \epsilon_t, \epsilon_t \sim WN(\sigma^2)$と定義されます。前項と攪乱項から値が決まります。前項が直接組み込まれていますので、当然自己相関が生じます。前項を参照するので、初項が必要となりますが、条件無し分布が分かっていればその分布に従う確率変数としたり、そうでなければ適当な定数を置いたりします。定常過程であれば初項の影響は時間と共に限りなく小さくなっていくので、あまり気にする必要はありません。ただし、一般に AR 過程は定常過程になるとは限りません。<br />
$AR(1)$が定常過程となる条件は、$|\phi_1| < 1$です。$|\phi_1| = 1$かつ$\epsilon_t \sim iid(\sigma^2)$の場合は単位根過程 (ランダムウォーク) となります。<a href="https://www.cyberer.net/2020/05/timeseries-analysis.html">時系列分析 (1)</a>で触れましたが、ランダムウォークは非定常過程です。$|\phi_1| > 1$のときは爆発的 (explosive) な過程と呼ばれます。どんどん$|y_t|$が大きくなっていきます。<br />
<br />
<h4>
定常なAR(1)の基本統計量</h4>
<b>定常な</b>$AR(1)$の平均は、$E(y_t) = E(c + \phi_1y_{t-1} + \epsilon_t) = c + \phi_1E(y_{t-1})$であり、定常ですから$E(y_t) = E(y_{t-1}) = \mu$と置くことができますので、これを$\mu$について解くと$\mu = \frac{c}{1 - \phi_1}$となります。<br />
分散は、$Var(y_t) = Var(c + \phi_1y_{t-1} + \epsilon_t) = \phi_12Var(y_{t-1}) + 2Cov(y_{t-1}, \epsilon_t) + Var(\epsilon_t) = \phi_12Var(y_{t-1}) + 0 + \sigma^2$であり、やはり定常性から$\gamma_0 = Var(y_t) = Var(y_{t-1})$と置いて$\gamma_0$について解くと、$\gamma_0 = \frac{\sigma^2}{1 - \phi_1^2}$となります。<br />
<br />
定常な$AR(1)$の自己相関係数について調べましょう。$k$次-自己共分散$\gamma_k$を求めると、$\gamma_k = Cov(y_t, y_{t-k}) = Cov(\phi_1y_{t-1} + \epsilon_t, y_{t-k}) = Cov(\phi_1y_{t-1}, y_{t-k}) + Cov(\epsilon_t, y_{t-k}) = \phi_1\gamma_{k-1}$となります。従って、$k$次-自己相関係数は、$\rho_k = \phi_1\rho_{k-1}$となります。これは Yule-Walker 方程式と呼ばれるものです。これを解けば、$\rho_k = \phi_1^k$であることがわかります。つまり、定常な$AR(1)$の自己相関係数の絶対値は、次数が上がるにつれて指数的に減少していきます。$\phi_1$が正であれば単調に、負であれば振動しながら減少していくことになります。<br />
<br />
<h3>
AR(p)</h3>
$AR(1)$が定常過程となる条件は$|\phi_1| < 1$でしたが、$AR(p)$に一般化した場合、定常となる条件は少し複雑です。$AR(p)$過程$y_t = c + \phi_1y_{t-1} + \phi_2y_{t-2} + \ldots + \phi_py_{t-p} + \epsilon_t, \epsilon_t \sim WN(\sigma^2)$が定常となる条件は、方程式$1 - \phi_1z - \phi_2z^2 - \ldots - \phi_pz^p = 0$の<b>全ての解の絶対値が1より大きい</b>ことです。この方程式を AR 特性方程式と呼びます。例えば、$AR(1)$であれば、AR 特性方程式は$1 - \phi_1z = 0$ となり、解は$z = \frac{1}{\phi_1}$です。$|z| > 1$となるのは$|\phi_1| < 1$のときですから、前述の通りの定常条件となりました。$AR(2)$ならば、$\phi_2 > -1, \phi_2 < 1 + \phi_1, \phi_2 < 1 - \phi_1$となります。<br />
<br />
定常な$AR(p)$の平均は、$\mu = E(y_t) = \frac{c}{1 - \phi_1 - \phi_2 - \ldots - \phi_p}$、分散は$\gamma_0 = Var(y_t) = \frac{\sigma^2}{1 - \phi_1\rho_1 - \phi_2\rho_2 - \ldots - \phi_p\rho_p}$です。$k$次-自己共分散は$\gamma_k = \phi_1\gamma_{k-1} + \phi_2\gamma_{k-2} + \ldots + \phi_p\gamma_{k-p} (k \ge 1)$、$k$次-自己相関係数は$\rho_k = \phi_1\rho_{k-1} + \phi_2\rho_{k-2} + \ldots + \phi_p\rho_{k-p} (k \ge 1)$となります。$AR(p)$においても$k$次-自己相関係数の式は Yule-Walker 方程式と呼ばれます。<br />
<br />
ちなみに、定常な AR 過程は、係数が指数的に減少するような$MA(\infty)$で書き直すことができます。例えば$AR(1)$なら、$y_t = \phi_1^0\epsilon_t + \phi_1^1\epsilon_{t-1} + \phi_1^2\epsilon_{t-2} + \ldots$と書き直すと、$\theta_k = \phi_1^k$とおいた$MA(\infty)$過程であることがわかります。<br />
<br />
<h2>
MA 過程と AR 過程の性質</h2>
MA 過程は定常なので扱いやすいのですが、$MA(q)$には$q$次までの自己相関しか組み込まれていませんから、逆に言うと、$q$次の自己相関をモデル化するには$MA(q)$が必要となります。これはモデルのパラメータが$q+2$個必要 ($\theta_1 \sim \theta_q と \mu と \sigma$) となることを意味しています。長期間の自己相関を扱うにあたっては、期間分のパラメータが必要となることは不都合があります。また、攪乱項の和で表現されていることから、モデルに意味付けをすることが困難であり、そのモデルから導かれる予測の扱いが難しくなります。<br />
一方、AR 過程は少ないパラメータで長期の自己相関を表現できます。特に$AR(p)$は周期$p / 2$の周期的な自己相関を織り込むことができ、これは季節変動があるデータのモデル化には好都合です。しかし、AR 過程では自己相関の強さが指数的に減少していくような形でしかモデル化できません。<br />
<br />
<h2>
ARMA 過程</h2>
ARMA 過程とは、自己回帰移動平均過程 (Autoregressive moving average process) のことです。AR 過程に移動平均項を追加したもの、つまり AR 過程と MA 過程を足し合わせたモデルです。$(p, q)$次のARMA 過程 ($ARMA(p, q)$) は、以下のように定式化されます。<br />
$$<br />
\begin{eqnarray*}<br />
y_t &=& c &+& \phi_1y_{t-1} &+& \phi_2y_{t-2} &+& \ldots &+& \phi_py_{t-p} \\<br />
&+& \epsilon_t &+& \theta_1\epsilon_{t-1} &+& \theta_2\epsilon_{t-2} &+& \ldots &+& \theta_q\epsilon_{t-q}, \epsilon_t \sim WN(\sigma^2)<br />
\end{eqnarray*}<br />
$$<br />
$ARMA(p, q)$は AR 過程の性質を持つため、常に定常になるとは限りませんが、定常な$ARMA(p, q)$は、以下の性質を持ちます。<br />
$$<br />
\begin{eqnarray*}<br />
E(y_t) &=& \frac{c}{1 - \phi_1 - \phi_2 - \ldots - \phi_p} \\<br />
\gamma_k &=& \phi_1\gamma_{k-1} + \phi_2\gamma_{k-2} + \ldots + \phi_p\gamma_{k-p} (k \ge q + 1) \\<br />
\rho_k &=& \phi_1\rho_{k-1} + \phi_2\rho_{k-2} + \ldots + \phi_p\rho_{k-p} (k \ge q + 1)<br />
\end{eqnarray*}<br />
$$<br />
なお、定常なAR過程同様、自己相関は指数的に減少します。<br />
<br />
$q$次までの自己相関は、移動平均項の影響で Yule-Walker 方程式が成り立ちませんが、$p>q$ならば$q+1$次以降の自己相関は自己回帰項の影響のみとなるので、$AR(p)$同様 Yule-Walker 方程式が成立します。<br />
<br />
ARMA 過程の定常性は、自己回帰項部分の定常性にのみ着目すれば十分です。ARMA 過程は AR 過程と MA 過程を足し合わせたものであり、MA 過程は常に定常であるからです。したがって、自己回帰項から AR 特性方程式を作り、その全ての解の絶対値が1より大きいことが定常性を持つ条件となります。<br />
<br />
一方、ARMA 過程の反転可能性は、自己回帰項の部分がすでに AR 過程で表現できているので$AR(\infty)$で表現できることは分かっていますから、残る移動平均項が AR 過程で表現できるかどうかによって決まります。したがって、移動平均項から MA 特性方程式を作り、その全ての解の絶対値が1より大きいことが反転可能性を持つ条件となります。<br />
<br />管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-85786233142810617222020-05-20T01:27:00.002+09:002020-05-28T09:55:51.501+09:00時系列分析 (1) - 時系列の性質の把握為替レートなどの価格時系列データの分析をよく行うので、統計学的に分析する場合の一般論について書いてみました。<br />
理系大学の学部生が使う教科書に書いてあるレベルより少し丁寧に説明しているつもりです。<br />
<br />
<a name='more'></a><h2>
時系列分析の流れ</h2>
<div>
時系列分析は、時系列データを分析して、何らかの回帰モデルを導出し、それに基づいて予測を行うことが主な目的になります。</div>
<div>
そのための手順は概ね以下のような流れになります。</div>
<ul>
<li>分析対象となる原系列を用意する</li>
<li>必要に応じて、対数系列、差分系列、対数差分系列などに変換 (主に定常性の仮定を満たすため) したり、季節変動要因 (seasonality) の調整を行ったりする</li>
<li>基本統計量で要約する </li>
<ul>
<li>期待値 (mean)</li>
<li>分散 (variance) = $0$次の自己共分散 (autocovariance)</li>
<li>$k$次-自己共分散 ($k$時点離れたデータ間の共分散) </li>
<li>$k$次-自己相関係数 (自己共分散の基準化) → $k$の関数として自己相関関数 (autocorrelation function/ACF) と呼ばれる </li>
<li>コレログラム (自己相関関数のグラフ) </li>
</ul>
<li>系列に対応する確率過程またはデータ生成過程 (DGP) を仮定する</li>
<li>過程を検証する </li>
<ul>
<li>時系列間の回帰分析を行う場合は、単位根検定によりその時系列が単位根過程でないことを確認しておく必要がある (単位根過程同士には見せかけの相関が発生することがある) </li>
<li>時系列$X$と$Y$がそれぞれ$I(d)$である ($d$階の階差が定常である) とき、$X$と$Y$の線形結合$(X - \alpha Y)$が$I(d-b)$であれば、$X$と$Y$は$(d, b)$次の共和分 (cointegration) 関係である</li>
</ul>
</ul>
<br />
<h2>
確率変数の基本統計量</h2>
<div>
まずは、一般の確率変数における基本統計量の導出方法を確認しておきましょう。</div>
<div>
<br /></div>
<div>
<h3>
期待値と分散、標準偏差</h3>
<div>
$E(x), Var(x) = E[(x-E(x))^2] = E(x^2) - (E(x))^2$は、確率変数$x$の期待値および分散です。分散は、値の散らばり具合を示しており、大きければ広い範囲に散らばっていることを意味します。また、$\sqrt{Var(x)}$を標準偏差といいます。</div>
<div>
<br /></div>
<h3>
共分散、相関係数</h3>
<div>
$Cov(x, y) = E[(x - E(x))(y - E(y))]$ を確率変数$x, y$の共分散と言います。$Var(x) = Cov(x, x)$ です。</div>
<div>
$Corr(x, y) = \frac{Cov(x, y)}{\sqrt{Var(x)}\sqrt{Var(y)}}$を2つの確率変数$x, y$の相関係数と言います。共分散の値を各確率変数の標準偏差の積で割ったもので、$-1 \le Corr(x, y) \le 1$の範囲の値になります。共分散は確率変数の分布によって異なる水準の値を取るため、そのままでは他の確率変数と比較できませんが、相関係数にすることで、$-1$は完全な負の相関、$1$は完全な同期、$0$は完全に独立といった意味を読み取ることができます。</div>
<div>
<br /></div>
</div>
<h2>
定常性 (stationarity)</h2>
<div>
時系列分析を行う上で、とても重要な性質が、定常性です。基本的に非定常時系列は時系列分析の対象とすることができないといっても過言ではありません。非定常時系列に対しては、見せかけの相関によって合理的な回帰モデルが構成できず、予測に使用することができません。定常性には、弱定常性と強定常性の2種類があります。<br />
<br />
以下の説明では、時点$t$という添字を持つ確率変数$y_t$からなる確率過程$Y = {y_t}$を対象とします。<br />
<br /></div>
<h3>
弱定常性 (weak stationarity)</h3>
過程$y_t$が任意の$t$に対して$E(y_t) = \mu$かつ、任意の$t, k$に対して$Cov(y_t, y_{t-k}) = E[(y_t - \mu)(y_{t-k} - \mu)] = \gamma_k$であるとき、$y_t$は弱定常過程であるといいます。$\gamma_k$ と書いているのは、$t$に依存せず、次数$k$にのみ依存した一定の値になる、ということを意味しています。<br />
<br />
<h3>
強定常性 (strict stationarity)</h3>
過程$y_t$が任意の$t$と$k$に対して、$(y_t, y_{t+1}, \ldots, y_{t+k})'$ ($(\cdots)$は転置ベクトル) の同時分布が同一となる場合、$y_t$は強定常過程であるといいます。分散が有限であれば、強定常過程は弱定常過程でもあります。<br />
<br />
<h3>
定常過程の性質</h3>
弱定常過程では、異なる時点のデータ間の線形依存関係が時点に依存しない (時間差のみに依存する) ことを必要としますが、強定常過程では、全ての形の依存構造が時点に依存しないことを必要とします。<br />
正規過程 (任意の$t$と$k$に対して、$(y_t, y_{t+1}, \ldots, y_{t+k})'$の同時分布が多変量正規分布となるような過程) においては、弱定常正規過程は強定常過程となりますが、一般には弱定常過程は強定常過程となるとは限りません。<br />
また、定常過程においても、<b>条件付き</b>期待値・分散が時点に依存しないということは言えません。<br />
<br />
<h3>
定常過程の探し方</h3>
非定常過程も差分過程を取ると定常過程になることがあります。対象の時系列が定常過程ではないときは、差分過程を確認しましょう。また、非定常となる要因が判明している場合は、まずその部分をモデル化し、残差を取り出して、それが定常過程であることを期待するという方法もあります。<br />
<br />
<h2>
定常過程の基本統計量</h2>
定常過程という前提を置くと、様々な有益な基本統計量を導出することができます。<br />
<br />
<h3>
平均と分散、標準偏差</h3>
$y_t$はそれぞれ確率変数ですから、平均$E(y_t)$、分散$Var(y_t)$、標準偏差$\sqrt{Var(y_t)}$を考えることができます。確率過程においては、平均を$\mu_t$、標準偏差を$\sigma_t$と表記することが多いです。<br />
ここで、$Y$が定常過程であるという前提を当てはめてみます。<br />
弱定常性の定義に$E(y_t) = \mu$とあった通り、定常過程においてはどの時点においても一定の期待値となります。<br />
分散についても弱定常性の定義に従うと、$Var(y_t) = Cov(y_t, y_t) = \gamma_0$ でやはり時点によらず一定です。<br />
<br />
<h3>
自己共分散、自己相関係数</h3>
確率過程においては、異なる2時点の確率変数における共分散、相関係数を考えることができます。それぞれ自己共分散、自己相関係数と呼ばれます。<br />
弱定常性の定義にある $Cov(y_t, y_{t-k}) = E[(y_t - \mu)(y_{t-k} - \mu)] = \gamma_k$ がまさに自己共分散です。$k$時点離れたデータ間の共分散を意味しており、$k$次-自己共分散と呼びます。どの時点$t$を選んでも共分散が一定であるというのが、定常過程の性質です。<br />
相関係数についても同様に$Corr(y_t, y_{t-k})$を考えることができます。$k$次-自己相関係数と言います。$Corr(y_t,y_{t-k}) = \frac{Cov(y_t, y_{t-k})}{\sqrt{Var(y_t)}\sqrt{Var(y_{t-k})}} = \frac{\gamma_k}{\gamma_0} = \rho_k$となり、やはり時点によらず一定です。<br />
<br />
噛み砕いて言えば、定常過程においては、期待値も自己共分散、自己相関係数も時点に依存しないということになります。<br />
<div>
<br /></div>
<h3>
自己相関関数、コレログラム</h3>
<div>
$k$次-自己相関係数$\rho_k$の系列を$k$の関数として見たものを自己相関関数 (autocorrelation function/ACF) と呼びます。また ACF をグラフに描画したものをコレログラムと呼んでいます。</div>
<div>
コレログラムを見ると確率過程の自己相関の特徴を捉えやすいので、便利です。コレログラムの例を2つほど見ていただきましょう。</div>
<div>
<br /></div>
<div>
次のグラフは、ある株価指数の前日比 (つまり株価指数の1階の階差) のコレログラムです。横軸の Lag が$k$にあたります。$0$次の自己相関係数は、同一の確率変数との相関で、完全に一致するわけですから、当然$1$です。$1$次から先は非常に小さい自己相関係数となっていることが分かります。</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/-DsJruhciSc8/XsP4Ydlcn0I/AAAAAAAAANQ/Fkp58RQjRAIMhcDimi4yymQ5V7-LVxAjgCLcBGAsYHQ/s1600/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-05-20%2B0.16.25.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="928" data-original-width="976" height="608" src="https://1.bp.blogspot.com/-DsJruhciSc8/XsP4Ydlcn0I/AAAAAAAAANQ/Fkp58RQjRAIMhcDimi4yymQ5V7-LVxAjgCLcBGAsYHQ/s640/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-05-20%2B0.16.25.png" width="640" /></a></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
次のグラフは、$\bf{CO}_2$濃度の前月比 (これも1階の階差) のコレログラムです。見ての通り自己相関係数が有意に$0$とは異なるように見えます。ですが、実は1年周期の季節変動要素を示しているだけで、$\bf{CO}_2$濃度の変化が前月の変化幅の影響を受けているとは言えません。このようなコレログラムが見えたときは、季節変動をモデル化して差し引き、残差に対して分析を続ける必要があります。</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://1.bp.blogspot.com/--64vXr2Zw3M/XsP5ZCV32SI/AAAAAAAAANc/PayhB-NJi9Q98tO62aTQs3A7u9D7V3WTQCLcBGAsYHQ/s1600/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-05-20%2B0.20.50.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="928" data-original-width="976" height="608" src="https://1.bp.blogspot.com/--64vXr2Zw3M/XsP5ZCV32SI/AAAAAAAAANc/PayhB-NJi9Q98tO62aTQs3A7u9D7V3WTQCLcBGAsYHQ/s640/%25E3%2582%25B9%25E3%2582%25AF%25E3%2583%25AA%25E3%2583%25BC%25E3%2583%25B3%25E3%2582%25B7%25E3%2583%25A7%25E3%2583%2583%25E3%2583%2588%2B2020-05-20%2B0.20.50.png" width="640" /></a></div>
<div>
<br /></div>
<h2>
ホワイトノイズ、iid過程</h2>
最も単純な弱定常過程として、時点によらず平均$0$、分散$\sigma^2$であるような過程があります。これをホワイトノイズと呼びます。$y_t$がホワイトノイズであるとき、$y_t \sim WN(\sigma^2)$と書きます。<br />
<br />
ホワイトノイズより強い仮定を置いた過程として、$iid$過程があります。これは、時点によらず平均$\mu$、分散$\sigma^2$の独立で<b>同一の分布</b>に従う過程です。$\mu = 0$の$iid$過程はホワイトノイズの性質を満たします。$y_t$が$iid$過程であるとき、$y_t \sim iid(\mu, \sigma^2)$と書きます。<br />
<br />
さらに強い仮定として、正規過程に従うホワイトノイズである、正規ホワイトノイズがあります。正規ホワイトノイズは正規分布に従うため、同一分布ですから、$iid$過程でもあります。$y_t$が正規ホワイトノイズ過程であるとき、$y_t \sim iidN(0, \sigma^2)$と書きます。<br />
<br />
<h2>
自己相関の検定</h2>
ある時系列データが定常過程に従うと仮定すると、標本統計量を用いて真の期待値、自己共分散、自己相関係数を推定することができます。期待値は標本平均、自己共分散は標本自己共分散、自己相関係数は標本自己相関係数がその推定値となります。<br />
ところで、時系列データに予測可能性があるかどうかは、その時系列と相関があり、先行している他の時系列があるか、あるいは自己相関があるかで決まってきます。そのため、標本から自己相関の有無を判別することは重要な関心事の1つです。<br />
<br />
自己相関の有無を検定するには、自己相関係数が$0$であるという帰無仮説、$0$でないという対立仮説をおいて、標本自己相関係数を用いて検定します。<br />
<br />
$iid$過程 (独立同一分布の過程) であれば、その自己相関係数は平均$0$、分散$1/T$ ($T$は時点数) の正規分布に漸近します。$k$次-標本自己相関係数が$-1.96/\sqrt{T} \sim 1.96/\sqrt{T}$の範囲を越えると、有意水準5%で帰無仮説は棄却され、$k$次-自己相関を持つと言えます ($1.96$は標準正規分布の両側95%点です)。前掲のコレログラムは、R の acf() 関数を用いてプロットしたものですが、この棄却点が青線でプロットされています。ci パラメータで有意水準を変更することができます。<br />
<br />
また、$m$次までの自己相関係数が全て$0$であるという帰無仮説を検定する手法があり、「かばん検定」と呼ばれています。Ljung and Box の方法が知られており、R では Box.test() を用いて行うことができます。この帰無仮説が棄却されるということは、$m$次までのどこかに自己相関があるということを意味します。具体的には$Q(m) = T(T+2) \sum_{k=1}^{m}(\rho_k^2 / (T-k))$が自由度$m$のカイ二乗分布に漸近することを利用して検定します。<br />
<br />
<h2>
単位根検定</h2>
過程$y_t$ が $y_t = y_{t-1} + e_t, e_t \sim iid(0, \sigma^2)$であるとき、この過程は「ランダムウォークである」、または「単位根を持つ」と言います。$E(y_t) = E(y_{t-1}) + E(e_t) = E(y_{t-1}) $ ですから($e_t \sim iid(0, \sigma^2)$なので$E(e_t)=0$です)、ランダムウォークは期待値が1つ前の時点に依存しています。すなわち非定常過程です。しかし、1階の階差を取ると$e_t$だけが残りますから、$iid$過程となります。<br />
これを調べるのが単位根検定です。Phillips-Perron 検定と Augmented Dickey-Fuller 検定が良く知られています。管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.comtag:blogger.com,1999:blog-8741021964326180750.post-37156022200565921722020-05-19T01:03:00.007+09:002020-05-19T01:03:52.826+09:00R - 研究編 5 メモリ管理に関する関数R のメモリ管理に関する関数をまとめて紹介します。<br />
<br />
<a name='more'></a><h2>
object.size(x)</h2>
オブジェクト x のサイズを object_size クラスで返します。実体はバイト数を表す整数です。<br />
別の記事で紹介した <a href="https://www.cyberer.net/2020/05/r-pkginfo.html">pkginfo()</a> でも使用しています。<br />
object_size クラス用の print() のオーバーライドがあり、units 属性で単位を指定できます。units に指定できる値は、"b", "Kb", "Mb", "Gb", "B", "KB", "MB", "GB", "auto" です。<br />
<br />
<h2>
gc(verbose=getOption("verbose"), reset=FALSE)</h2>
ガーベジコレクションを実行します。verbose を TRUE にすると、cons cell, vector の各領域についての統計情報を出力します。reset を TRUE にすると、メモリの最大使用量をリセットします。<br />
<br />
<h2>
gcinfo(verbose)</h2>
verbose を TRUE にして実行すると、以後ガーベジコレクションが実行される度に統計情報を出力するようになります。FALSE にすれば出力されなくなります。<br />
<br />
<h2>
gctorture(on=TRUE)</h2>
on を TRUE にして実行すると、以後メモリを確保しようとする度にガーベジコレクションを実行するようになります。超遅くなります。<br />
<br />
<h2>
gctorture2(step, wait=step, inhibit_release=FALSE)</h2>
step 回メモリ確保を行う度にガーベジコレクションを実行するように設定します。wait を指定すると、その回数分メモリ確保が行われるまでガーベジコレクションの強制実行を開始するのを待ちます。<br />
inhibit_release を TRUE にすると、メモリを解放せず、再利用に備えるようになります。<br />
<br />
<h2>
memory.size(max=FALSE)</h2>
Windows 専用です。現在の使用メモリ量を取得します (max=FALSE)。max を TRUE にすると、OS から報告される割り当て済の最大メモリ量になります。max を NA にすると、メモリ上限量になります。<br />
<br />
<h2>
memory.limit(size=NA)</h2>
Windows 専用です。メモリ上限量を取得します (size=NA)。size に数値を指定すると、その値を MB として解釈し、メモリ上限量を変更します。<br />
<br />
<h2>
memory.profile()</h2>
cons cell の使用状況を SEXPREC のタイプごとに取得します。<br />
<br />
<h2>
Rprofmem(filename="Rprofmem.out", append=FALSE, threshold=0)</h2>
filename で指定されたファイルにメモリ割り当て状況を出力します。NULL, 空文字列にすると出力されなくなります。append を TRUE にすると既存のファイルに追記します。<br />
threshold に数値を指定すると、その値をバイト数として解釈し、そのサイズ以上のメモリ割り当てのみが記録されるようになります。<br />
<br />
<h2>
tracemem(x)</h2>
オブジェクト x を監視し、コピーが発生したときにその旨を表示するようにします。監視を特定する文字列を返します。<br />
<br />
<h2>
untracemem(x)</h2>
tracemem(x) によって設定したオブジェクト x の監視を外します。<br />
<br />
<h2>
retracemem(x, previous=NULL)</h2>
オブジェクト x の監視を行うにあたって、x を生成する元となっているオブジェクトとの紐付けを行います。previous には元となっているオブジェクトの監視特定文字列を指定します。<br />
一般的には、tracemem(a); b <- foo(a); retracemem(b, retracemem(a)) のように使用します。<br />
<br />
管理人http://www.blogger.com/profile/11928044228828325974noreply@blogger.com