霞と側杖を食らう

ほしいものです。なにかいただけるとしあわせです。[https://www.amazon.jp/hz/wishlist/ls/2EIEFTV4IKSIJ?ref_=wl_share]

怪異・妖怪伝承との衝動遊戯:第1ゲーム

【ABSTRACT】

怪異の専門家に手っ取り早くなるために、怪異・妖怪伝承データベースのデータに対して、Rとテキストマイニング(自然言語処理データマイニング)の手法を用いて、怪異について学習する。

【OBJECTIVE】

モラトリアムが終わろうとしている。このまま就職して本当に良いのだろうか。そんな疑問が頭を過ぎる。就活をしているとき、自分の過去と向き合って、やりたいことを探していた。面接でしゃべるためだけの志望理由を探すために。
そんな過去への省察をもう一度やってみた。すると、一人の男が思い浮かんだ。忍野メメである。忍野メメは、西尾維新の『化物語』と、それに続くシリーズの登場人物で、怪異・妖怪変化の類の専門家である(化物語は、高校生の頃、初めて見た深夜アニメだったように思う。表現の仕方が独特で、アニメが面白いものだと気付くきっかけだった。このままでは脱線してしまうので、これ以上は語らない。)。 忍野メメは一見、胡散臭いが、実際は頼りになる人物で、とても魅力的に感じていた。そうだ、怪異の専門家になってみよう。それが今回の遊びの目的である。
ところで、怪異とは何なのか。忍野メメによると、
「人間がそこにあると思うからそこにあるもの、それが怪異」
とのことである。

【METHODS】

データは、国際日本文化研究センターが作成した『怪異・妖怪伝承データベース』を利用する。 国際日本文化研究センター | 怪異・妖怪伝承データベース
データを検索すると、それぞれのカードに、番号、呼称、出典、話者(引用文献)、地域、要約、貼り付け事例がある。要約の項の説明に、「貼り付け事例からピックアップした民俗語彙や特徴、行為、怪異・妖怪現象の本質などを、100字程度の文章形式にまとめたもの。なお、データには現在の観点からみて不適切と判断される表現も含まれているが、データの資料性を重視して、貼り付け事例に従って表記した。」とある。要約のところだけをスクレイピングで抽出して集めることにする。

データベースには、全部で35701件の事例データがあるが、さすがに全部扱うのは厳しい。代表的な怪異である「オニ」を取り扱う。「オニ」には「ウシオニ」、「オニババ」、「オニビ」などあるが、「オニ」として登録されているものだけを取り扱う。全部で496件である。

そのうちの1件を具体的に見てみることにする。これはオニの検索結果で一番上に表示されたデータである。 http://www.nichibun.ac.jp/YoukaiCard/0640010.shtml
要約の項目には、「昔、鬼が毎年現れて田畑を荒らした。村人は困り、もし一夜のうちに100段の石段を作れなければ、以後、姿を現さないと鬼に約束させた。ある夜、鬼が石段をつくり、あと一つで100段というところで、一番鶏が鳴いたので、鬼は姿を隠してしまった。」とある。

これらの要約部分を抽出して集めていく。スクレイピングはrvestパッケージで、read_html()とhtml_table()で、できると思ったのだが、496件中105件、謎の文字化けが発生してしまった(文字化けって怪異では?)。解決方法が分からなかったので、別の方法を探すことに。XMLパッケージのreadHTMLTable()も使えるみたいなので、やってみたら、文字化けを起こさずに上手くいった。無事、要約を抽出した496件のtxtデータが入手できた。このへんのコードは第2ゲームに載せる予定。

さて、データは得られた。35701件を496件に絞り込んだと言えども、全部を読んで頭に入れている時間は無いので、テキストマイニングの力を借りて、オニという怪異の特徴を把握する。頻度集計、ワードクラウドでオニ全体の特徴をつかみ、トピックモデルを用いて、各オニを分類してみて、傾向をつかんでみる。

【RESULTS】

頻度集計

オニ全体の特徴をざっくりと把握するため、496件のデータをひとまとめにして頻度集計を行う。 まず、自立語(動詞、形容詞、形容動詞)についての頻度集計の結果がこちら。
f:id:moratoriamuo271:20181024033654p:plain
食べる、死ぬ、殺す 、さらうなどのおどろおどろしい言葉が見られる。
次に、名詞についての頻度集計の結果がこちら。
f:id:moratoriamuo271:20181024034149p:plain
鬼が多いのは当然として、退治するケースが多そう。他に、節分、豆、正月、菖蒲といった言葉が目立ち、鬼の来る季節というのがありそう。月が気になってデータベースを検索してみたら、月と関係する話もあるが、大半は何月の話かという月が多かった。夜という言葉も多い。オニは夜行性のようだ。目や首は、オニの話なのか、人の話なのか、ここからは分からない。

ワードクラウド

はてなブログ映えを意識して、ワードクラウドによる可視化も行ってみた。これもひとまとめにしたデータで行った。
f:id:moratoriamuo271:20181024035104p:plain
オニのイメージはこんな感じということだ。

トピックモデル

496件のオニに関する怪異・妖怪伝承データをトピックモデルによる分類を行った。
トピック数は30にして、よく出る言葉を並べたものを用意した。30にしたのは、3,4,5,10,20,30と色々やってみて、解釈しやすかったから。 f:id:moratoriamuo271:20181024050434p:plain
トピック1は夫婦ものの話だろう。金と血が出ているのが少し怖い。トピック2は季節は秋で兄弟ものだ。
トピック3は鬼の形相を書いているものだろうか。赤や青、角、歯と鬼の描写かと思われる。トピック4は晩に山からやってきて妹か息子を鬼が殺してしまった話か。 トピック5は自然と関係する。トピック6は池周辺の家の話か。 トピック7は正月の話で、ショウブが関係する。菖蒲といえば、端午の節句では、と思ってググってみたら、食わず女房という話が日本にはあるらしい。これだ。
食わず女房 - Wikipedia
トピック8は山と狐と集落と仕業が気になる。トピック9はトイトイが気になったが、麻雀は関係なく、鳥のことらしい。
トピック10は正月に門松で厄除けという話だろう。トピック11は部落の神社と弓の話か。
トピック12は退治の話だろうか。トピック13は節分の話だ。トピック14はお地蔵さんに団子でもお供えして鬼から身を守ったのだろうか。トピック15は鬼は外、福は内の話だ。
トピック16は僧が鬼を追い払うかと。トピック17は旅人が泊まった家の山姥の正体が鬼だったとかではないか。天井から隠れて見ていた話かもしれない。
トピック18は綱が気になってデータベースを見てみたら、綱と名がつく人の話がいくつか見られた。渡辺綱という名前がよく出てきて調べてみると、
渡辺綱 - Wikipedia
平安時代の武将で、源綱ともいう。羅生門の話と関係があり、一条戻橋が舞台の話の書き換えなんだとか。
トピック19は霧のかかった川を船で渡っていたら子供が連れ去られたとかだろうか。トピック20は節分にイワシの頭と柊を戸口にって話だ。トピック21は剣で目玉を刺したとかか。
トピック22は端午の節句だ。菖蒲。ヨモギ。トピック23は罪人が地獄に落ちて、鬼にいじめられているところを神が救おうとしたとかではないか。
トピック24は大晦日が気になった。データベースを覗いてみると、大晦日は鬼が家に入らないように各地でいろんな工夫をしているみたいだった。
トピック25は相撲が気になった。鬼と相撲するケースがいくつかあるようだ。 朝鮮の話では、鬼と箒が関係あるらしい。中でも笑ったのが、この話で、
http://www.nichibun.ac.jp/YoukaiCard/1231055.shtml
「ある人が鬼と相撲を取った。倒してみると、正体は便所の箒に鼻血がついたものだった。」というものだ。
トピック26は庄屋が気になる。庄屋とは村長のような存在らしい。村に鬼がやってきて困ったとかって話か。
トピック27は勅命。退治する命令でも下ったのだろう。トピック28は山鳥が気になる。夢のお告げや山鳥の尾で作った矢が関係するらしい。
トピック29は酒呑みが関係しそうである。興味深い。トピック30は鶏の鳴きと博打が気になる。博打を調べてみると鬼が博打してる話が何件も出た。
トピックの解釈はこんな感じだろう。分類してみると、新しい発見もあるものだ。

【CONCLUSIONS】

スクレイピングでは、文字化けという名の怪異に憑りつかれて、なかなか眠りにつけなかった。周囲の友人からも、最近元気ないねと心配される次第であった。
分析上の課題としては、盛り沢山である。
word2vecは試してみたい。時間の制約上、今回は、できなかった。
トピックモデルについては、LDAでトピック数を雰囲気で決めてしまったが、ldatuningというパッケージがあるらしいので、これを使って、トピック数を決定してみるべきだった。
ディリクレ過程を使ったノンパラベイズによるトピックモデルならば、トピック数を事前に決めなくても良いらしいのだが、まだキャッチアップできていない。というか、そもそもトピックモデルを理論的に把握していないので、学習に努めたいところ。
解釈には、やはりドメイン知識が不可欠だ。柳田國男あたりを読んで伝承に詳しくなるべきかもしれない。
今回はオニにのみ注目したが、他の知らない怪異を見てみるのもよいかもしれない。さらに、最近発見したのが、 まんが日本昔ばなしのデータベースだ。
まんが日本昔ばなし〜データベース〜 - 1975年〜1994年、TBS系列で放送されたTVアニメ「まんが日本昔ばなし」の全話総まとめ、全話あらすじデータベースです。
これを使ってスクレイピングして、分析するのも面白そうだ。

コードや分析の詳細は別記事に追記する予定で、第2ゲームで、RMarkdownでHTML変換して貼り付ける予定。
【2018/10/25に貼り付けました。】
怪異・妖怪伝承との衝動遊戯:第2ゲーム - 霞と側杖を食らう

スクレイピングテキストマイニングの入門には多少成功したようには思うし、スクレイピングテキストマイニングの力を借りて、怪異の専門家に少しは近付けたように思う。怪異の専門家の力を借りたくなったら、ご連絡ください。力は貸します。助けはしません。一人で勝手に助かってもらうだけですが。

ところで、
「人間がそこにあると思うからそこにあるもの、それが怪異」
という怪異の説明を最初に挙げた。
これは統計モデリングと似ていないだろうか。
何か不思議な現象があったとする。怪異の専門家ならば、現象を、見えない何かによるものとして、その何かを怪異と呼ぶだろう。 統計家ならば、現象を、確率分布(統計モデル)によるデータ発生と捉えるだろう。怪異の専門家も統計家もやっていることの本質は同じなのであった。

【REFERENCES】

・『R Advent Calendar 2017 rvestを用いてポケモンデータをスクレイピング&分析してみた』
R Advent Calendar 2017 rvestを用いてポケモンデータをスクレイピング&分析してみた | かものはしの分析ブログ
スクレイピングの際に、参考にしました。

・『readHTMLTable:R から HTML の表を読み込む』
readHTMLTable:R から HTML の表を読み込む - 廿TT
XMLパッケージのreadHTMLTable()で、tableの読み取りができることに気付かせていただきました。

・小林『Rによるやさしいテキストマイニング: 機械学習編』
Amazon CAPTCHA
これを読んでテキストマイニングの仕方を学びました。読んだときの記録は昔書きました。
Rによるテキストマイニング(機械学周編)の学習記録 - 霞と側杖を食らう

・石田・市川・瓜生・湯谷『Rによるスクレイピング入門』
http://amzn.asia/d/clK9Tjm
スクレイピングの学習に使っています。

国際日本文化研究センター『怪異・妖怪伝承データベース』
国際日本文化研究センター | 怪異・妖怪伝承データベース
データベースを利用させていただきました。

動学的最適化の基礎の学習記録

【学習動機】

変分ベイズを学ぶ前に、変分法のおさらいをしようと思ったことと、マクロ経済学で、動学的最適化問題に帰着させて、最適化を行うのだけれど、変分法だとかハミルトニアンだとか動的計画法(DP)だとか、出てくるけど、いったいぜんたいなんなんだってなって整理したかった。とくに、中級マクロをやると、二期間モデルで最適化解きますまでは良いのだが、唐突に「こいつがオイラー方程式って呼ばれる」って書かれる。
たとえば、
二神・堀『マクロ経済学
http://amzn.asia/d/bC2ZtZE
では、たしかそのようだったと記憶している。オイラー方程式ってなんなんだよってなるのに、それを教えてはくれない。とても気持ち悪い(二神・堀『マクロ経済学』は、それ以外は読みやすく、中級マクロ本としては良著だと思います。)。上級マクロの講義を受けると、最適化を、オイラー方程式だったり、ハミルトニアンだったり、ベルマン方程式だったりの、どれかで、突然解かれていて、何がなんだか分からないうちにどんどん前に進んでいってしまうことがある(とくに一番最初は面食らって意味不明の沼の底に落とされた。)。一応、解き方を真似してその場しのぎをしてみるけれど、やはり、気持ち悪い。
動学的最適化の諸々の違いは何なのか。各手法の旨味は何なのか。そのへんが分からなくて気持ち悪い。
そんなときに、日本語でいいものがないか探していて出会ったのが、
A.C.チャン『動学的最適化の基礎』
http://amzn.asia/d/5FvrgpB
だった。そんなわけで、これをざっくりと読むことにした。

【学習記録】

<総括>

動学的最適化問題には、3つの主要なアプローチがあり、変分法と最適制御理論(ハミルトニアン)と動的計画法がある。
変分法は古典的アプローチで、汎関数積分可能で、登場するすべての関数が連続で連続的に微分可能であることが仮定される。
最適制御理論は、時間変数 tと状態変数 y(t)に加えて、制御変数 u(t)を扱う。制御変数に焦点を当てるということは、状態変数が二次的な地位になる。これは、状態変数について初期条件が与えられると、制御経路の決定が、その副産物として状態変数の経路を必ず決める場合にのみ受け入れられる。このため、最適制御問題は状態変数yと制御変数uを結ぶ式である \frac{dy}{dt}=f(t,y(t),u(t))がないといけない。これは運動方程式(遷移式、状態方程式)と呼ばれ、uの選択によるyへの影響を表す。 動的計画法の特徴は、与えられた制御問題を制御問題の一群に包含し、その結果それはある与えられた問題を解くことが結果的に問題の一軍を解くことになることと、この一連の問題の焦点は汎関数の最適値 V^{ * }に向けられ、変分法における最適状態経路 y^{ * } や最適制御理論における最適制御経路  u^{  *  }  の特徴に関心がない点である。

この本は、第一部が序論、第二部が変分法、第三部が最適制御理論で構成されている。動的計画法については、連続時間だと偏微分方程式の知識が必要になるため、深入りせず、変分法と最適制御理論に焦点が当てられている。この記事では、各章について記述するのではなく、気になっていたことについてだけ書いていくことにする。

結局、 オイラー方程式とは、 \max{\int{F(t, y(t), y'(t))}dt} s.t.  y(0)=A, y(T)=Z変分法で解いたときに、変分法で設定された \epsilon p(t)が消え、必要条件として現れる F_y - \frac{d}{dt}  F_{y'} = 0 \forall t \in [0,T ]のことだった。(このブログの記事では場合分けして書くのが面倒なので、積分区間を曖昧に濁している。本ではしっかり書かれている。)このオイラー方程式は、 F_{y' y'} y''(t) + F_{y y'} y'(t)+F_{t y'}-F_{y} = 0 \forall t \in [0,T ]とも表現でき、これはオイラー方程式が一般に、2階の非線形微分方程式であることを意味している。これは。一般解は2つの任意の定数をもつから、初期点と終点があれば確定した解が得られるってことにつながる。終点が可動的な場合を考えると、横断性条件が導ける。
制約条件付き静学的最適化問題同様、制約条件付き動学的最適化問題でも、ラグランジュの未定乗数法が役に立つ。ラグランジュ乗数を、それぞれが1つのオイラー方程式(オイラーラグランジュ方程式)における付加的状態変数として扱えばよい。動学的最適化問題に帰着して、ラグランジュ関数を立てて解いてるのは、変分法だったわけだ。

最適制御理論の特徴の一つは、変分法は内点解しか扱えない一方で、最適制御理論ではコーナー解も扱える点である。制御経路はジャンプする非連続性が認められている。他方で、状態経路は時間期間を通して連続でなければならないが、有限個の鋭点、角を持つことが認められていて、区分的に微分可能であればよい。最適制御理論のもう一つの特徴は、制御変数 uに対する制約を直接扱えること。もうひとつ、変分法との違いは、もっとも簡単な例が、変分法では固定終点のケースであったが、最適制御理論では自由終点(垂直的終点直線)のケースである点。
最大値原理と呼ばれる1階の必要条件が重要で、ハミルトン関数(ハミルトニアン) Hと共役状態変数(costate variable、あるいは補助変数)  \lambda の概念によって記述される。共役状態変数は、ラグランジュ乗数と似ていて、性質上、評価変数であり、関連する状態変数の潜在価格(シャドープライス)を測るもの。ハミルトニアン
 H(t,y,u,\lambda )=F(t,y,u) + \lambda (t) f(t,y,u)
と定義される。状態変数 yに関する1つの2階微分方程式であったのがオイラー方程式だったが、最大値原理は、状態変数 yと共役状態変数 \lambdaに関する2つの1階微分方程式を含んでいる。さらに、あらゆる時点において、制御変数 uに関してハミルトニアンが最大化されることが必要となる。
たとえば、最適制御の問として最も簡単なものを
 \max{V=\int{F(t,y,u)}dt}
s.t.  \dot{y}=f(t,y,u),y(0)=A,y(T)=Z Aは所与でZは自由。および、 \forall t, u(t) \in U
と記述でき、ハミルトニアン
 H(t,y,u,\lambda )=F(t,y,u) + \lambda (t) f(t,y,u)
となり、最大値原理の条件は、
 \max_{u}{H(t,y,u,\lambda)} \forall t
 \dot{y}=\frac{\partial H}{\partial \lambda}
 \dot{\lambda}=-\frac{\partial H}{\partial y}
 \lambda(T)=0
となる。
最大値原理の詳細な証明はこの本では扱っていないが、制御問題の変分法的な見方で以て、理論的根拠を示している。変分法と最適制御理論の比較をし、両者の同等性についても述べられている。その他、経済学的な解釈の仕方まで学べる。 たとえば、利潤を最大化する企業を考え、状態変数は資本ストック K、制御変数 uは経営上の意思決定として、最適制御問題を設定すると、
 \max{\Pi} = \int{\pi (t,K,u)}dt
s.t.  \dot{K} = f(t,K,u),  K(0)=K_0, K(T)=K_T ( K_0, Tは所与、 K_Tは自由。)
このとき、
 \frac{\partial \Pi^{ * }}{\partial K_0} =\lambda ^{ * } (0)
 \frac{\partial \Pi^{ * }}{\partial K^{ * } (T)} = - \lambda ^{ * } (T) となる。最適な共役状態変数 \lambda ^{ * } (0)は、所与の初期資本に対する最適総利潤 \Pi^{ * }の感度の尺度となっており、初期資本を1単位増やしたら総利潤が \lambda ^{ * } (0)大きかったことを意味する。したがって、初期資本の帰属価値、潜在価格であると理解できる。逆に、   \lambda ^{ * } (T)は計画期間の最後であと1単位保存したいとすると総利潤を  \lambda ^{ * } (T)だけ犠牲にしなくてはいけない。つまり、終点資本ストックの潜在価格を意味する。
同様に、一般的に、 \lambda ^{ * } (t)はその特定時点における資本ストックの潜在価格となる。

紛れもなく、読みやすい良著。上級マクロで初見殺しを食らう前に知っておきたかった。

<捕捉>

マクロの消費と投資については、 阿部修人先生講義ノート
http://www.ier.hit-u.ac.jp/~nabe/2007consumption.pdf

http://www.ier.hit-u.ac.jp/~nabe/investment2013.pdf
が流れを追いやすく参考にした。

橘永久先生の最大値原理の資料
http://www.le.chiba-u.ac.jp/~ttachi/doc/maximumprinciple-061124.pdf
もとても参考になった。

物理のかぎしっぽの変分法のページは今回読まなかったが、過去に変分法で躓いたときに参考した。
http://hooktail.sub.jp/mathInPhys/variations1/

経済セミナーの最適化特集回も昔読んだので、記録しておく。
経済セミナー2011年10・11月号|日本評論社

後から見つけたので、読めてはいないが、
須賀晃一先生の経済数学の資料
http://www.f.waseda.jp/ksuga/matheco05II.pdf
が、とても有用そうな気がした。
他に、細矢祐誉先生の経済学のための数学の講義資料がすごいので、 http://stairlimit.html.xdomain.jp/text16-8.pdf
つよい人は読んでみるとよいかと。

最適化について、もっと初めの方から始めたい方は
金谷『これなら分かる最適化数学―基礎原理から計算手法まで』
http://amzn.asia/d/8l4kbdI
がオススメ。

<追記>

ルジャンドル変換が気になって少し調べて雰囲気理解しただけだけど、参考になった資料のリンクを追加する。
ルジャンドル変換とはなにか?(動画 バージョン)
http://irobutsu.a.la9.jp/mybook/ykwkrAM/sim/LegendreTR.html
ルジャンドル変換とは何か(Legendre transformation) http://fnorio.com/0146Legendre_transformation/Legendre_transformation.html
EMANの物理学 ルジャンドル変換
https://eman-physics.net/analytic/legendre.html

【学習予定】

マクロのお気持ちは少しずつ理解できてきた気がする。 変分ベイズにいい加減、取り組みたい。

Rによるテキストマイニング(機械学周編)の学習記録

【学習動機】

これまでベイズ統計を学んできたけど、ノンパラベイズをまだ知らない。ノンパラベイズについて調べてみると、ディリクレ過程が出てきて、ディリクレについて調べてみると、トピックモデルが出てきた。トピックモデルで何かしたいというわけじゃないが、理論を追ってトピックモデルは行けたらいいなと思った。
トピックモデルは自然言語処理に関連する。機械学習関連分野に、自然言語処理と画像認識があるが、どちらかできるようにしておくといいなと思っていた。でも、画像認識は出来れば面白そうだけどマシンパワーないと辛そうだし、自然言語処理文字コードとか面倒だからいいやと放置していた。今回、自然言語処理と触れ合えるチャンスなので、入門してみようという気持ちが起こった。できるようになれば、趣味としても面白いような気もするし、実力ついたら、それはそれで将来どこかで仕事に使えるかもしれない。そんなわけで、とりあえずRで基本的なことをできるようにしようと思い立ったのである。

【学習記録】

<総括>

小林『Rによるやさしいテキストマイニング[機械学習編]』を読むことにした。

https://www.amazon.co.jp/R%E3%81%AB%E3%82%88%E3%82%8B%E3%82%84%E3%81%95%E3%81%97%E3%81%84%E3%83%86%E3%82%AD%E3%82%B9%E3%83%88%E3%83%9E%E3%82%A4%E3%83%8B%E3%83%B3%E3%82%B0-%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E7%B7%A8-%E5%B0%8F%E6%9E%97-%E9%9B%84%E4%B8%80%E9%83%8E/dp/4274221008

200ページほどで、自然言語処理、前処理、スクレイピング機械学習のRを用いたやり方が分かるのは良い。駆け抜けていくような本。ざっくり手法を知るのには適していると思われる。コードの重複するような部分はサポートページに任せて、関数の説明をもう少し詳しくすればベターだったと感じた。この『Rによるやさしいテキストマイニング』はシリーズ化されており、無印版と活用事例編版がある。活用事例編は気になるところ。

サポートページに、付属データがないんだけど、、、
https://www.ohmsha.co.jp/book/list.htm?keyword=r%E3%81%AB%E3%82%88%E3%82%8B%E3%82%84%E3%81%95%E3%81%97%E3%81%84%E3%83%86%E3%82%AD%E3%82%B9%E3%83%88%E3%83%9E%E3%82%A4%E3%83%8B%E3%83%B3%E3%82%B0&cmid=752

と思ったら、ここにあった。

ML - Rによるやさしいテキストマイニング

第1章から第4章までがテキストマイニングの基本を扱っており、第5章から第7章が機械学習の手法を扱っている。

c1:自然言語処理

形態素解析は、文を単語単位で分割。各単語の品詞などを同定する処理。
日本語の形態素解析ツールとして、MeCabChaSen、JUMAN、KyTeaが挙げられている。
規則や定義がツールごとに違うので、1つのプロジェクトでは同一のものを使うべし。

【MeCab】を【R】で使えるようにするための【RMeCab】を導入する方法・流れ【Windows編】|さぎのみや家の分析録

これに従って、MeCabとRMeCabのインストールをした。日本語はこれでなんとかなりそう。
英文解析では、パッケージのopenNLPとNLPを使っている。
やってみると、library(openNLP)がうまくいかない。
error: JAVA_HOME cannot be determined from the Registry
とエラーを吐いてくる。windowsのデフォルトのjavaが32bitなのが悪いらしい(?)
rJavaがロードできない - 盆栽日記 これで解決できたぽい。

構文解析(係り受け解析)は、文を区単位で分割し、区の係り受け関係(修飾・非修飾関係)を分析する技術。
日本語の構文解析ツールとして、CaBoChaとKNPが挙げられている。 英語では、openNLPmodels.enパッケージとopenNLPパッケージのParse_Annotator()で解析可能。

・意味解析は、いくつか種類があって、固有表現抽出、術後項構造解析、語義曖昧性解消、評判分析(感情分析)がある。ここでは、評判分析の例が紹介されている。ポジ・ネガの単語が入ってる評価表現辞書が必要。 英文の評価分析をするならば、tidytextが有用。パッケージに3つの辞書(afinn, bing, nrc)があり、言葉とそれに対応する感情(nrcならば、anger, fear, positive, negativeなど10種類)のデータフレームとして入っている。

head(get_sentiments("afinn"))
summary(get_sentiments("afinn")$score)

と確認したところ、afinn辞書はwordに対応するscoreがmax5のmin-5で与えられている。

head(get_sentiments("bing"))
unique(get_sentiments("bing")$sentiment)

と確認したところ、bing辞書はwordに対応するposiとnegaのみ。
tidytextは単語単位の評判分析ができたが、文単位の評判分析はSentimentAnalysisパッケージで可能。本書の第一版第一刷時点では、CRANに公開されていなかったみたいだが、2018/10/11現在では公開されていて、普通にインストールできた。
日本語の評価表現辞書としては、単語感情極性対応表、日本語評価極性辞書が挙げられている。

R+RMeCabで感情分析

このへんを参考にすれば日本語はできそう。

・言語判定は、テキストが何の言語か判定する技術。textcatパッケージで可能。

・文書要約は、重要な情報を自動で抽出する技術。LSAfunパッケージのgenericSummary()で可能。

c2:テキスト処理

テキストの読み込みと整形と頻度集計と用例検索を扱っている。

テキストの読み込みはscan()で、整形はstringrやtmで小文字変換、句読点、数字削除。tmのstopwordsも。日本語なら、stringrのstr_replace_all()で。頻度集計は日本語ならRMeCabTxet()で得た結果をmap_chrの第二引数をextract(1)と指定して、table()で集計して、sort()で並べ替え、データフレーム化してあげればいい。もしくはRMeCabFreq()を使ってもできる。
複数文書の頻度集計すると文書ターム行列が得られる。(NLPの理論本ぱらぱら見てて、潜在意味解析(LSA)の特異値分解のとこで出てきたやつか!)
英語なら、tmのCourpus()を整形して、DocumentTermMatrix()で文書ターム行列作成して、as.matrix()でゲット。
日本語なら、RMeCabのdocDFでできる。 n-gramの頻度集計、日本語なら、RMeCabのNgram()で可能。共起語の頻度集計は、日本語なら、RMecabのcollocate()で。 用例検索はforループで。

c3:スクレイピング

rvestを使用。これは、もう少し詳しく掘り下げたいため、今回は簡単にコード回して別の機会に。

『Rによるスクレイピング入門 』 https://www.amazon.co.jp/R%E3%81%AB%E3%82%88%E3%82%8B%E3%82%B9%E3%82%AF%E3%83%AC%E3%82%A4%E3%83%94%E3%83%B3%E3%82%B0%E5%85%A5%E9%96%80-%E7%9F%B3%E7%94%B0-%E5%9F%BA%E5%BA%83/dp/486354216X

で学びたいところ。

c4:データハンドリング

これも、宇宙本で過去に学習しているため、流し見しただけ。テキストに対してどのように使うかは追加的に学べる。

Rstudioとtidyverseの学習記録 - 霞と側杖を食らう

c5:回帰

この章から機械学習の手法のパートになるが、だいたいの手法は知っているため、コードを書いて流し見。気になったところはメモする形になる。 p119のcolnames(BNCbiber)の結果が違う。たとえば、"f_01_past_tense"は"f01_past_tense"。 p120のselect(BNCbiber, ~~~)は、エラーが出る。

Error in UseMethod("select_") : 
  no applicable method for 'select_' applied to an object of class "c('matrix', 'double', 'numeric')"

classがmatrixだとダメなのでは?仕方がないので、as.data.frame(BNCbiber)でデータフレームに変換したら解決。
平滑化回帰の交差検証のところで、bisoregパッケージのloess.wrapperをやっても、返ってこない、なぜか不明。
正則化回帰はglmnetパッケージ使ってやってる。問題ない。

c6:分類

k近傍法をclassパッケージのknn()でやっている。knn(trainの説明変数, testの説明変数, trainの目的変数, kの数) 書かれていないが、caretパッケージのconfusionMatrix()を使うにはe1071パッケージも必要みたい。
ナイーブベイズをe1071パッケージのnaiveBayes()でやっている。
ニューラルネットワークをnnetパッケージで。SVMをe1071パッケージのsvm()で。
ちなみに、e1071パッケージの名前は、" Misc Functions of the Department of Statistics, Probability Theory Group (Formerly: E1071), TU Wien"とあり、このパッケージを作ったウィーン工科の研究チームの識別子がe1071だったことに由来するみたい。

https://cran.r-project.org/web/packages/e1071/e1071.pdf

決定木のCARTをrpartパッケージとrpart.plotパッケージを使っている。
バギングはadabagパッケージのbagging()、ブースティングはadabagパッケージのboosting()、勾配ブースティングはgbmパッケージのgbm()を使用。ランダムフォレストはrandomForestパッケージのrandomForest()。

このへんのアンサンブル学習は、なんかよく使われているらしいけど、何も理論的背景を知らないので、入門してみたいところ。

c7:教師なし学習

階層型クラスタリングはdist()で距離を得て、hclustに、得た距離をつっこんで、methodを指定すればオーケー。
階層型クラスタリングは可視化性能は良いが、計算量が多く、大きいデータには使いにくい。 kmeans()でk-means法。
自己組織化マップはkohonenパッケージのsomgrid()とsom()を使用。
word2vecをwordVectorsパッケージで実行。やってみたけど、たぶん、メモリ不足で上手くいかなかった。 トピックモデル 1つの文書なら、ldaパッケージで可能。日本語でやるなら、RMeCabで形態素解析済みのデータを作ってやる必要がある。
複数文書なら、topicmodelsパッケージを使うとよい。

【学習予定】

トピックモデルが気になるので、理論含めて学んでいきたい。

佐藤『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』

https://www.amazon.co.jp/%E3%83%88%E3%83%94%E3%83%83%E3%82%AF%E3%83%A2%E3%83%87%E3%83%AB%E3%81%AB%E3%82%88%E3%82%8B%E7%B5%B1%E8%A8%88%E7%9A%84%E6%BD%9C%E5%9C%A8%E6%84%8F%E5%91%B3%E8%A7%A3%E6%9E%90-%E8%87%AA%E7%84%B6%E8%A8%80%E8%AA%9E%E5%87%A6%E7%90%86%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA-%E4%BD%90%E8%97%A4%E4%B8%80%E8%AA%A0/dp/4339027588/ref=pd_bxgy_14_img_2?_encoding=UTF8&pd_rd_i=4339027588&pd_rd_r=1c761cb3-d025-11e8-a088-e38b52a8e10c&pd_rd_w=cjY9N&pd_rd_wg=hAy82&pf_rd_i=desktop-dp-sims&pf_rd_m=AN1VRQENFRJN5&pf_rd_p=a4de75e6-d8f7-4a34-bd69-503ea4866e6c&pf_rd_r=YMC0NPW68TZJ139H80GE&pf_rd_s=desktop-dp-sims&pf_rd_t=40701&psc=1&refRID=YMC0NPW68TZJ139H80GE

岩田『トピックモデル (機械学習プロフェッショナルシリーズ)』

https://www.amazon.co.jp/%E3%83%88%E3%83%94%E3%83%83%E3%82%AF%E3%83%A2%E3%83%87%E3%83%AB-%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%83%97%E3%83%AD%E3%83%95%E3%82%A7%E3%83%83%E3%82%B7%E3%83%A7%E3%83%8A%E3%83%AB%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA-%E5%B2%A9%E7%94%B0-%E5%85%B7%E6%B2%BB/dp/4061529048/ref=pd_rhf_dp_s_cp_0_6?_encoding=UTF8&pd_rd_i=4061529048&pd_rd_r=32a053ef-9690-4827-b938-0ee907c2dc6d&pd_rd_w=KXxPv&pd_rd_wg=wwj3v&psc=1&refRID=1S99Q54MSNH93XTCV1QR

あたりで学習していきつつ、実際に何かしらのデータを用いて試していきたい所存。

VSCodeやGIt・Github、Markdownの整備の覚書

【用途】 

エディタの選択、EmacsとかAtomとかあるけど、VSCode使いやすそうではということでVSCode整備。今のところ、pythonmarkdownとかが書ければいいかなくらい。RはRStudioで良いし。このへんを整備するが、 また整備する機会のために個人的メモとして、記録しとく。あと、整備で詰まったところとか書き残してる。

 

【内容】 

VSCode

元からsurfaceに入ってたので、それを使用。

コマンドパレットは Ctrl + Shift + P

 

python

python anaconda3

www.messyer813.com

 

分かりやすい。

 

Jupyter をVSCode上で展開したくて、後程、拡張でそれに関係するものを入れる。

 

<Git・Github

Gitが入ってなかったのと、Githubのアカウントは昔作ったまま放置してた。

これに従って入れて、整備していくことにした。

employment.en-japan.com

 

vimのcondfigの終わらせ方が分からなくて一瞬戸惑った。

qiita.com

 

hacknote.jp

 

qiita.com

 

リポジトリのクローンで

Permission denied (publickey)ってエラーが出て、困った。

hacknote.jp

基本的にこれに従えばよいのだが、「`ssh-agent`」を「"ssh-agent"」と打っていたせいで、上手くいかなかった。"と`は別物なのね。。。

 

コミットしようとしたら、who are you?って聞かれた

qiita.com

これでなんとかなった。

コミットでコメント残したと思ったら、aborting commit due to empty commit messageと言われる。

qiita.com

どうやらこういうことらしい。あと#の後にスペースが要る。

  

<拡張>

このへんから欲しい拡張機能を選んでインストールした

ics.media

qiita.com

 

vscode-icons

アイコンついて見やすい。

・Clipboard Ring

コピーを複数記憶できる。

 ・Bracket Pair Colorizer

カッコが色付きになって判別しやすくなる。

・Path Autocomplete

パスの入力が楽に。

・Rainbow CSV

csvファイルを見るときに列で色別になって見やすくなる。

・Partial Diff

選択範囲の差分を調べられる。写経のときにスペルミスとかで溶かしてた時間を削減できるはず。

・Trailing Spaces

スペースを強調してくれる。変なスペースで溶かしてきた時間ともこれでお別れ。

 

Markdown

 マークダウン関連の拡張

 ・Markdown All in One

mdファイルを開いた際に自動的にプレビューを起動(デフォルトでは無効)。

ショートカットキーでMarkdownのタグを入力可能(ctrl+bで太字、ctrl+shift+]でインデントの追加など)。

・Auto-Open Markdown Preview

自動でプレビュー起動。

Markdown PDF

mdからpdfにエクスポート。mdファイル保存して開いて右クリックで選択。

これの問題点として、数式入っているmdファイルをpdfにしても、数式を変換してくれない。なので、数式使うときは下のMarkdown Preview Enhancedを使用する。

Markdown Preview Enhanced

プレビューがつよくなってる。数式入りMarkdownもpdfにできる。直接吐き出す方法は分からないけれど、ブラウザで開く(open in browser)で、chromeで開き、印刷でpdfで保存とすれば出来るみたい。

Markdown + Math

数式入れるため。$$内に数式を書いてctrl+shift+.と打って、数式対応のプレビューが表示。

qiita.com

次に、数式を複数行書くとき、=を揃えて書こうとしてつまづいた。

\begin{aligned}

a &=b\\

&=c

\end{aligned}

で、できた。KaTeXだと、{align}じゃなくて、{aligned}だった。

tera-kun.hatenadiary.com

・Paste Images

 画像をコピーしたら、Ctrl+Alt+V でMarkdownファイル内に貼り付けられ、同一ディレクトリに保存。(保存先は設定ファイルで変更可能。)

・Jupyter

Jupyterを使うため。

VS Code Jupyter Notebook Previewer

プレビューに使うため。

・ Tables Generator

excelの表を変換してくれるサービス。markdownで表作るのってめんどくさいけど、これがあれば楽になるはず。

www.procrasist.com

 

【記憶の検索キーワード】

VSCode, GIt, GIthub, Markdown, 拡張, 整備, エディタ

 

Rによる多重代入法の学習記録

【学習動機】

Kaggleで欠測値を埋めなきゃということがある。他にも埋めなきゃいけない状況がある。
欠測値の理論は、因果推論の一般的な枠組みとして捉えるのに役立つ。
あと、ベイズと関わってくるので、勉強しておきたいと思った。実際に、Rで実行できるように。

【学習記録】

高橋・渡辺『欠測データ処理 Rによる単一代入法と多重代入法』で学んでいくこととする。

<総括>

Rで実際にどうやって欠測値に多重代入法を用いて対処するかが分かったのは、とても良い。
今出ている日本語書籍では、理論を学ぶことはできるが、どれをどうやって使うのかがすぐに分からなくて困っていた。
データセットは書籍HPに用意すべき。データは書面に記載しているってのはモダンじゃない。
norm2とmiceとAmelia。代入・分析・統合。

c1:Rによるデータ解析

Rの使い方と欠測あるときの反応を少しチェックしている。

c2:不完全データの統計解析

多くの統計ソフトウェアでは、リストワイズ除去で、欠測値を含む観測対象の行をすべて削除し、擬似的「完全」な状態にして分析することが多い。
欠測パターンと欠測メカニズムで分類される。
欠測パターンには、単変量欠測パターン、単調欠測パターン、計画的欠測パターン、一般的な欠測パターン。多重代入法では、いずれの欠測パターンにも有効で、欠測パターンはそんなに重視されない。
欠測メカニズムには、完全に無作為な欠測(MCAR:Missing Completely At Random)と条件付きで無作為な欠測(MAR:Missing At Random)と無作為ではない欠測(NMAR:Not Missing At Random)がある。 Dをn×pのデータセット D = {  D _ {obs}, D _ {miss} }。Kは回答指示行列でn×p行列。 K _ {i j}=1で観測、  K _ {ij}=0で欠測
・MCAR
 Pr(K|D) = Pr(K) であり、欠測する確率が、データと独立。
・MAR
 Pr(K|D) = Pr(K|D _ {obs}) で、観測データを条件とし、欠測確率の分布が非観測データから独立。欠測データがMARならば、欠測を除去する分析は偏っている。この偏りは共変量(covariate)または補助変数(auxiliary variable)による代入で是正可能。 ・NMAR
 Pr(K|D) \neq Pr(K|D _ {obs}) で、欠測する確率が欠測した変数の値自体に依存している。代入法で偏りを是正できるとは限らず、選択モデルやパターン混合モデルを用いた分析が必要で、感度分析が有用。

分析の目的が母集団の推定なら、多重代入法を使用すべき。

c3:単一代入法

従来、公的統計では、調査データの集計が主目的で、確定的単一代入法を用いることが通例。方法として、よく用いられるのが、回帰代入法、比率代入法、平均値代入法、ホットデック法。
・確定的単一代入法はガウス・マルコフの仮定あるとよい。
・比率代入法は不均一分散のとき使える。
・平均値代入法は百害あって一利なし。
・ホットデック法(hot deck imputation)は、ノンパラな代入法。ホットデックでは、欠測している人の補助変数の情報を利用して、回答者から似通った値のものを選び出し、その回答者の値を代入値として採用。何を基準にするかは議論あり。最近隣法が多い。
HotDeckImputationパッケージのimpute.NN_HD()で実行できる。
・確率的回帰代入法
各々の代入値に誤差項を加える。誤差項は、平均が0、分散が回帰モデルの残差分散の正規分布から無作為抽出するなどのやり方があり、miceのmice()で引数methをnorm.nobとすれば実行できる。

c4:多重代入法の概要

ガウス・マルコフの仮定とMARの過程がすべて満たされた場合、多重代入法によるパラメータ推定量はBLUEになる。
代入、分析、統合の3ステージ
統合方法。  \tilde{\theta_m}をパラメータ \thetaのm番目の代入済みデータセットに基づいた推定値とする。
統合した点推定値は
 \bar{\theta_M} = \frac{1}{M} \Sigma \tilde{\theta_m}
多重代入法の代入値の分散は代入間分散(between)と代入内分散(within)に分解でき、
代入内分散は、
 \bar {W _ m}  = \frac{1}{M}  \Sigma var(\tilde{\theta_m})
代入間分散は、
 \bar{B _ M} = \frac{1}{M-1} \Sigma (\tilde {\theta _ m} -  \bar{\theta_M}) ^2
 \bar{\theta_M}の分散 T_Mは、これらを合算して
 T_M = \bar{W_m} + (1 + \frac{1}{M} ) \bar{B_M}

この方法をRubin(1987)のルールという。

代入間の繰り返しをどれだけ行えば収束するかは、様々な要素が絡んでくるが、とりわけ欠測情報の比率が重要。
適切な多重代入法(proper multiple imputation)でないと、推測の結果が妥当でなくなるおそれあり。
代入モデルと分析モデルが同一の変数を持ち、同じ数のパラメータを推定するとき、2つのモデルには適合性(congeniality)があるという。適合していない場合、多重代入法のパラメータ推定量の一致性は保証されない。
代入モデルよりも大きな分析モデルを使うと、リストワイズ除去より悪い結果になりうる。
また、多重代入済みデータ数は多い方が望ましいとされる。

c5:多重代入法のアルゴリズム

3つのアルゴリズム
・DAアルゴリズム。norm2でできる。しかし、norm2が消えてる。
FCSアルゴリズム。miceでできる。これが主流か?
・EMBアルゴリズム。Ameliaでできる。 (EMアルゴリズムについてはまた別の機会で詳しく学ぶ)

norm2パッケージがCRANから消えてるので

url <- "https://cran.r-project.org/src/contrib/Archive/norm2/norm2_2.0.1.tar.gz"
pkgFile <- "norm2_2.0.1.tar.gz"
download.file(url = url, destfile = pkgFile)
install.packages(pkgs=pkgFile, type="source", repos=NULL)
unlink(pkgFile)

あと、ホットデックのコードを動かした後、Ameliaパッケージの多重代入するとエラー起きるかもしれないらしい。
Amelia、norm2、miceと見比べて、すべての変数が量的な連続変数なら、どれも変わらず、miceが時間的に効率悪いか。
質的データや時系列データになると、アルゴリズムの使い分けが必要。

c6:多重代入モデルの診断

MARが正しいと仮定したとき、代入の結果がMARの仮定と整合性があるかを診断する。
MCARの場合を除くと、観測データの分布と代入データの分布は一致せず、むしろ、異なっていることが期待される。診断で、分布の違いを発見し、MARの仮定と整合性があるか確認する。分布が極端に違う場合も、どの変数に問題があるか見つけるのに診断ツールが役立つ。 Ameliaとmiceの診断ツールを用いる。

c7,8:量的データの多重代入法

・多重代入法を用いた平均値のt検定 AmeliaとMKmiscを使うか、miceのpool.scalar()の返り値から直接計算。
欠測データによって、自由度の計算が大変。分散の比率λと代入の不確実性の分散の相対的な増加rを計算する。

このへん、本のコードが変。本が書かれた後にアップデートされて関数が変わったのか知らないが。 miceのpool.scalar()は引数にmethodはとらないし、返り値にlambdaはない。
rは返ってくるので、 \lambda = r/(1+r) で自分で計算しないとかと。

・重回帰分析
回帰診断で、ジャック・ベラの正規性検定(JB検定)で誤差項の正規性を検証、ブルーシュ・ペイガン検定(BP検定)で不均一分散の検証、分散拡大要因(VIF:Variance Inflation Factor)で多重共線性の検証、ボンフェロー二外れ値検定で外れ値の検証。

c9,10:質的データの多重代入法

・ダミー変数のある重回帰分析 質的データの代入法に関して、決定打はないようだ。ここでは、応用研究で評判が高いとされるmiceによる方法とhot.deckの多重ホットデックによる方法が紹介されている。

・ロジスティック回帰
miceでいつも通り代入して、glm.mids()の引数familiyをbinomialとすればよさげ。

c11:時系列データの多重代入法

一般的に、時系列データにおける欠測値は、補完(interpolation)やカルマンフィルタを用いて処理されることが多い。
一方で、多重代入法を用いて欠測値を処理することで不確実性を反映した分析を可能にする。

c12:パネルデータの多重代入法

Ameliaで実行可能。amelia()の引数csに横断面の単位、引数tsに時系列の単位、引数polytimeに3以内の多項式を指定、引数lagsにラグ付き従属変数を指定。未来の値から過去の値を予測することで代入モデルの精度を向上させたいなら引数leadsを指定してもよい。
Ameliaに入っているafricaのデータを用いて、plmでパネルデータ分析をする。
p151のコードを回してみたが、変量効果モデルの推定がうまくいかない。特異行列が出て、solve()で逆行列を求めるところがうまくいってないみたいだが。。。

c13:感度分析:NMARの統計分析

SensMice(SensMiceDA)パッケージを使えば、感度分析は可能。感度分析の目的は、欠測データに関する設定を変えて、分析結果にどう影響するかを調べること。デルタ修正のパラメータデルタをどう設定するかが重要。デルタの設定方法として、変数の標準偏差を基準とする方法がある。標準偏差の-1倍, -0.5倍, 0倍, 0.5倍, 1倍した5つのデルタ値を用意して、どう影響するかを見て、ドメイン知識と照合して検討。

c14:事前分布の導入

Ameliaならamelia()の引数boundsに(変数の位置、最小値、最大値)を行にした行列を指定することができる。リッジ事前分布を使ってモデルの安定性を実現したいなら、amelia()の引数empriで指定。
norm2ならば、emNorm()のprior="ridge",prior.df=0.5などとする。miceはsqueeze()を用いる。

【学習予定】

高井・星野・野間『欠測データの統計科学』

www.amazon.co.jp

阿部『欠測データの統計解析』

www.amazon.co.jp

あたりで理論を補強していきたい。

制限付き最尤法(REML)の学習記録

【学習動機】

  統計数理研究所公開講座で、経時的反復測定データの取り扱い方で、MIxed models for repeated measures(MMRM法)について
習った。そのとき、パラメータ推定の仕方として、REML法が出てきたが、いまいちよく分かってなかったため。
制限付き最尤法(REML: REstricted Maximum Likelihood)(REMLは残差最尤法(REsidual Maximum Likelihood)ともいうみたい。罰則付き最尤法とは別物。)を、土居正明先生のpdf資料を用いて、学んでいく。
直交行列、Helmert変換、多変量正規分布、正準形あたりが前提知識となっており、それらも同様にして、準備していく。
その後、まだ理解が不足しているので、追加的に随時、補足する。

【学習記録】

土居正明『直交行列と Helmert 変換』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/orthogonal_mtx.pdf

直交行列の定義。 P’P=I_n
命題。n×n行列Pが直交行列であることはPが \mathbb{R}^nのある正規直交基底であることは同値。
性質。 O(n)={P: n×n行列|P’P|=I_n}として、
 P \in O(n)に対して、群の性質 (a)-(c)に加えて、(d)と(e)がある。

 (a) I_n \in O(n)
 (b) P,Q \in O(n) \Rightarrow PQ \in O(n)
 (c)  P^{-1} \in O(n)
 (d)  P’ = P^{-1}
 (e) |P| = \pm1

土居正明『一次独立と正規直交基底』講義資料
http://www012.upp.so-net.ne.jp/doi/math/anova/linear_indep.pdf
Parsevalの等式は射影と相性が良い。正準形の理論でも役立つ。
 \mathbb{R}^n において正規直交基底\{ e_1, ... , e_n \} をとると、]
 \mathbb{R}^n中の任意のベクトルaに対して、以下の(i)(ii)が成り立つ。
 (i) a = (a \cdot e_2)e_1 + (a \cdot e_2)e_2 + ... + (a \cdot e_n)e_n
 (ii) || a ||^2 =  {(a \cdot e_1)}^2 + {(a \cdot e_2)}^2 + ...  + {(a \cdot e_n)}^2

Parsevalの等式は射影と相性が良い。正準形の理論でも役立つ。

特殊な直交行列として、Helmert行列の作り方。
(i) ある手順で \mathbb{R}^nの直交基底を作る。
(ii) (i)で作った直交基底を正規直交基底に正規化する。
(iii) (ii)で作った正規直交基底を横に並べて直交行列にする。
(iv) (iii)で作った行列を転置する。

こうして作られたHelmert行列 H_nを用いた変換をHelmert変換という。
(エルミート行列(Hermitian matrix)のことか?と思ったけど、違うね。)
https://en.wikipedia.org/wiki/Helmert_transformation

これは、 1) \bar{X}  {s}^2 は互いに独立である。 2) n{s}^2/{\sigma}^2 ~{\chi}^2(n-1) を同時に示すのに使える。

c.f. 竹村数理統計p67,68

土居正明『多変量正規分布』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/m_normal.pdf

多変量正規分布は多変量標準正規分布を一次変換して、定数ベクトルを加えたもの。
「定理4:標準正規分布の標本平均と標本分散の独立性」
「系2:一般の正規分布の標本平均と標本分散の独立性」
上の1)と2)のこと。

土居正明『正準形1』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/canonical_form1.pdf

群の数を Iとし、 n_i (i=1,2,...,I)とおき、 \Sigma{n_i} = Nとする。
 y= X \beta + \epsilon,    \epsilon ~  N(0, \sigma ^2 I_N)
ただし、 \beta = (\mu_1, \mu_2, ..., \mu_I )

このとき、適当な直交行列G(Helmert行列で考えたのと同様)でyをzに変換すると
 z_1 ~  N(\sqrt{n_1}\mu_1, \sigma ^2),  z_I ~  N(\sqrt{n_I}\mu_I, \sigma ^2)
 z_{I+1}, ..., z_N ~  N(0, \sigma ^2)
で、すべて独立で、
 \hat{\mu_1} = z_1 / \sqrt{n_1},..., \hat{\mu_I} = z_I / \sqrt{n_I},
 \hat{\sigma ^2} =  \frac{1}{N-I}  \Sigma _ {i=I+1} ^{N}  z_i ^2
となり、 \hat{\mu_1}, \hat{\mu_2}, ..., \hat{\mu_I}, \hat{\sigma ^2}はすべて独立。
 \hat{\mu_i} ~  N( \mu_i, \frac{1}{n_i}  \sigma ^2)
 (N-I) \frac{\hat{\sigma ^2}} {\sigma ^2}  = \Sigma _ {i=I+1} ^{N} \frac{z_i ^2}{\sigma ^2} ~  \chi ^2 (N-I)
となる。

 \beta と \sigma ^2 のUMVUである\hat{\beta}と \hat{\sigma ^2}が互いに独立であり、それぞれ従う分布が確認できた。

(UMVUとは)
 \theta の不偏推定量 \hat{\theta} が一様最小分散不偏推定量(UMVU : Uniformly Minimum Variance Unbiased Estimator)であるとは  \theta の任意の不偏推定量 \hat{\theta}に対して、
 V_{\theta}  (\hat{\theta})
 \leq V_{\theta} (\hat{\theta}),  \forall \theta
が成り立つときをいう。
(土居正明『不偏推定量・UMVU・大数の法則』資料2 http://www012.upp.so-net.ne.jp/doi/math/anova/UMVU.pdf より)

土居正明『正準形2』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/canonical_form2.pdf
分散分析のF検定の分子と分母の定数倍が帰無仮説のもとで独立なchi ^2分布に従うことを示して、
F統計量が帰無仮説のもとでF分布に従うことを示している。

土居正明『制限付き最尤法(REML)』分散分析資料
http://www012.upp.so-net.ne.jp/doi/math/anova/REML.pdf REMLの問題意識は、正規分布から得られたデータからの、分散の最尤推定量は不偏推定量にならない。
正準形を用いてデータを2つに分け、平均と分散の別々の尤度関数を立てて推定する方法を制限付き最尤法という。
これは、xについての尤度関数を、正準形にしてzについての尤度関数の積に分解したものと考えられる。

島健悟『線形推測論』
http://nshi.jp/contents/doc/glm/glm2016_10.pdf

ここまでの説明ではランダム効果の説明が入っていない。

長島先生がREMLの導出のメモをしているものを見つけました。
https://nshi.jp/contents/stat/reml/

REML以外にも良い資料がたくさんありました。
https://nshi.jp/contents/

南美穂子『制限付き最尤推定法(REML推定法)』 https://www.jstage.jst.go.jp/article/jappstat1971/25/2/25_2_73/_article/-char/ja
REMLに関する覚書き。この方は岩波DSで野球のバントの因果推論をしていたような気がします。

田豊『線形混合モデルにおける分散成分の推定』分散成分の推定と歴史
http://www.obihiro.ac.jp/~suzukim/masuda2/stat/review_variance.html
帯広畜産大学のwebサイトのリニューアルで消えたらしいのですが、増田豊先生がミラーではてなブログに転載しているのを見つけました。

yutakamasuda.hatenablog.com

上の記事は2007年に書かれたもので、2019年の状況を踏まえた注釈の記事も書かれていました。
yutakamasuda.hatenablog.com

田豊『制限付き最尤法』
https://slidesplayer.net/slide/11509221/
混合モデルのことを分散成分モデル(variance componets model)と呼んでいた。 REMLの推定のアルゴリズム
導関数を使用しない方法(Derivative Free REML)
・1次導関数を使用する方法(EM REML)
・2次導関数を使用する方法(AI REML)
・AI REML以外にも,いくつかの準ニュートン法による分散成分の推定法が提案されている
・REML以外のアプローチとして、ベイズ

黒木玄さんの不偏分散の分母がn-1の理由に関するマストドンのpdf https://genkuroki.github.io/documents/mathtodon/2017-05-16%20%E4%B8%8D%E5%81%8F%E5%88%86%E6%95%A3%E3%82%84%E4%B8%8D%E5%81%8F%E5%85%B1%E5%88%86%E6%95%A3%E3%81%AE%E5%AE%9A%E7%BE%A9%E3%81%AB%E3%81%8A%E3%81%84%E3%81%A6%E3%81%A9%E3%81%86%E3%81%97%E3%81%A6n%E3%81%A7%E5%89%B2%E3%82%89%E3%81%9A%E3%81%ABn%E2%88%921%E3%81%A7%E5%89%B2%E3%82%8B%E3%81%8B.pdf
n-1になるのは線形代数の直交射影で空間の次元が落ちるからという話。統計学の教科書はなかなか教えてくれない、自由度が何たるかを理解する助けになる。REMLの話も、上で見てきたように、同様の線形代数なので、この話もとても有用。

浅野・中村『計量経済学
https://www.amazon.co.jp/%E8%A8%88%E9%87%8F%E7%B5%8C%E6%B8%88%E5%AD%A6-y21-%E6%B5%85%E9%87%8E-%E7%9A%99/dp/4641163367
の12章の12.2に「分散要素モデル」あり。あとで読む。

【学習予定】

REMLがまだよくわかっていない。ランダム効果モデル、混合モデルについて、もう少し詳しくなるべきで、
みどりぼんやアヒル本、『一般化線形モデル入門』や『統計モデル入門』を読もうと思う。

Rstudioとtidyverseの学習記録

【学習動機】

moratoriamuo.hatenablog.com

これのときに、気付いたのが、Hadley教徒になりたいけど、ぐぐって断片的に集めた知識じゃ、すぐ抜け落ちるってこと。というわけで、本で体系的に学ぼうということで、通称、宇宙本をやることにした。

 

松村 優哉『RユーザのためのRStudio[実践]入門−tidyverseによるモダンな分析フローの世界−』

http://amzn.asia/d/639ydJK

 

【学習記録】

<総括>

RStudioやtidyverseの入門として、とてもよかった。次に、各章の専門的に書かれた本に飛べばよいはず。

 

<忘れやすいショートカットめも>

Ctrl + Shift + M ;パイプ演算子%>%

Ctrl + Shift + C ;選択した行をコメント/非コメント

Ctrl + Alt + K :レンダリング実行(Knit)

Ctrl + Alt + I :Rチャンク挿入

 

c1:RStudioの基礎

・ファイルの読み込みは標準関数のread.csv()ではなく、readrのread_csv()を使うべし。

他に、data.tableのfread()も使えるが、ここでは扱わない。

標準関数のread.csv()は遅いし、読み込む型が不親切。

・読み込む列の型を指定したいなら、col_types引数を使う。

文字列[c]、整数[i]、実数[d]、TRUEFALSE[l]、日付[D]、時間[t]、日付時間[T]、数字以外の文字が含まれても数字として返す[n]、推測する[?]、列を読まない[_]

エンコーディングの指定。UTF-8で読み込まれるが、Shift-JIS(CP932)のファイルを読むと文字化けするので、local引数で指定。文字コードが分からないときは、guess_encoding()を使うといい。

・書き出しにはwrite_csv()を使う。標準関数のwrite.csv()より高速で、デフォルトでrow.names = FALSEなので、行名が入らなくてよい。

Excelファイルの読み込みにはreadxlのread_excel()を使えばよく、sheet引数でシートの指定ができる。

SASSPSS、STATAファイルの読み込みはhavenパッケージが使える。

 

c2:スクレイピングによるデータ収集

スクレイピングにはHTMLとCSSの知識必要。

・URLのHTML構造を抽出し、要素にアクセスし、その中身を取り出すか、属性にアクセスし、その値を取り出す。CSSセレクタXMLパスのような構文を利用したスクレイピングによってもドキュメントの要素や値を取得できる。

スクレイピングを実行するためにはrvestパッケージを使用。

まず、URLをread_htmlで読み込む。読み込むと、DOM(Document Object Model)という形で保存される。これはHTMLの要素を階層構造に変換したもの。HTMLの要素やクラスにあたる部分、DOMでノードという。DOMはブラウザでチェック可。

ページソースを読み解きやCSSセレクタXpathの文法暗記は困難なので、スクレイピングではデベロッパーツールというchromeの開発ツールを利用する。WindowsではFn+F12で使え、ノードが確認できる。ここから取得したいところにカーソルあてて、Copy XPathでコピー。クオーテーションマークはシングルとダブルを区別することに注意。

 

・ページを複数取りたい場合は、繰り返しを使う。URLが規則的な構造ならforループでなんとかなるが、そうでない場合はRSeleniumパッケージを使う。

(RSeleniumがCRANから消えている。2018/09/04現在。)

API(Application Programming Interface)は、つまり、ソフトウェアやサービスを外部から使用できるようにしたもの。無制限には使えず、単位時間での使用量決まっていること多い。

Twitterには、REST APIとSTREAMING APIがある。前者はタイムラインやDM、FF、リストなど取得できる。検索も可能。後者は流れるデータの1%を受信できる。

https://togetter.com/li/1250205

仕様変更していて、異なるところがあるみたい。

 

『Rによるスクレイピング入門』

 

で補強していくべき。

 

c3:dplyr/tidyrによるデータ前処理

 

Hadleyの提唱するtidy dataとは、

・1つの列が1つの変数を表す。

・1つの行が1つの観測を表す

・1つのテーブルが1つのデータセットだけを含む

 

(横長のワイド(wide)形式のデータを縦長のロング(long)形式のデータに変えるなど。パネルデータ分析の前処理で出てきたりする操作。)

ロング形式の方が、ベクトル演算しやすい。

データをtidy dataに変形するのに、tidyrのgather()とspread()が有用。

ワイドをロングに変形するのはgather()。

ロングをワイドに変形するのはspread()。

引数keyとvalueを指定、gather()なら並べる列名のkeyも指定、spread()はkeyが列名になるので指定不要。

 

列や行の絞り込みといったデータ操作はdplyrパッケージが有用、dplyrの主な関数は3つに大別でき、

・1つのデータフレームを操作する関数  行を絞るfilter()や列を絞るselect()など。

・2つのデータフレームを結合する関数  行同士で結合のjoin()系関数、集合演算intersect()やnion()など。

・ベクトルを操作する関数  ウィンドウ関数など。

 

・select()は引数に絞りたい列名を列挙。ヘルパ関数として、starts_with(), ends_with(), contains(), matches(), num_range(), one_of(),everything()が使える。

・filter()は引数に条件式を書いて絞り込み。論理式はいつも通り。

・mutate()は新しい列名を追加、transmute()は追加でなく別のtibble作成。

・arrange()が行を並べ替える。選択した列の昇順がデフォ。数値列なら-をつければ降順。そうでなければ、desc()でくくるといい。

・summarise()で集計計算。グループ化と組み合わせると強力。データを1行に要約する関数。

・group_by()を適用すると、データはグループ化されたデータフレーム(group_df)になり、mutate()やfilter()を使うと違う挙動を見せる。

 

セマンティクスに関しては、これ

https://speakerdeck.com/yutannihilation/dplyrfalseselecttomutatefalsesemanteikusufalsewei-i

 

・tidy dataに変換できないデータなら、複数の列を同時に処理するscoped functionが有用

・scoped functionは以上の関数の複数列バージョンで、サフィックスは_all, _if, _at。

これ、なかなか便利では!?

 

・2つのデータセットの結合と絞り込みで、_join()の関数を使う。

 

c4:ggplot2によるデータ可視化

カラーページで描かれていてよい。

パッケージggplot2の利点は、統一された文法、グラフの再現性、豊富な拡張パッケージ。

水準ごとの識別

・グループ化(mapping引数のaes()で指定。)

・ファセット(facet_wrap()とfacet_grid())

・統計的処理(stat())

水準ごとに自動で配色を割り当てる機能は可視化の効率化に大きく貢献する一方で、配布資料や印刷物でグラフを示す場合、モノクロやグレースケールが望ましいこともある。

背景真っ白なら、theme_classic()、目盛り線を残したいなら、theme_bw()、文字サイズやフォントの変更はtheme_()の中で変更。

 

 

c5:R Markdownによるレポート生成

ここのパートの著者が書いてるものが、

R Markdown入門

https://kazutan.github.io/kazutanR/Rmd_intro.html

 あと、これも役立つ。

www.amazon.co.jp

 

 

【学習予定】

 データハンドリングとスクレイピングをもう少し習得したい。

www.amazon.co.jp

 か

www.amazon.co.jp

 をやっていこうと思う。