第59期王位戦挑戦者決定リーグ、紅組は羽生善治竜王が優勝 [将棋]
昨年失陥してしまった羽生竜王ですが、王位戦挑戦者決定リーグの最終局で松尾歩八段に勝利しました。
紅組での優勝なので、この後、白組の優勝者と勝負をして勝つことができれば、王位戦の挑戦者になることができます。
さて、最終局は、全て同時対局となるわけですが、羽生竜王はなんと55手という最短手順で勝利することができました。
横歩取りの将棋から、お互い大技を掛け合って、さて、これからどういう構想で相手玉を攻めるか、、、という状況の中で、いきなり、香車を捨て、馬も捨てて、バッサリと松尾玉を仕留めました。
将棋ソフトも最初、何が起こっているのか分からないまま、評価値が下がったのですが、その後、羽生さんの手順が進んでいくと、あっさりと手のひらを返して羽生さんが優勢の評価となりました。
将棋ソフト(AI)の読み筋をも超えた鬼神の将棋でした。
この将棋の内容なら、今の将棋ソフト(AI)にも勝つことができるんじゃないか?って思ってしまった羽生先生の将棋でした。
紅組での優勝なので、この後、白組の優勝者と勝負をして勝つことができれば、王位戦の挑戦者になることができます。
さて、最終局は、全て同時対局となるわけですが、羽生竜王はなんと55手という最短手順で勝利することができました。
横歩取りの将棋から、お互い大技を掛け合って、さて、これからどういう構想で相手玉を攻めるか、、、という状況の中で、いきなり、香車を捨て、馬も捨てて、バッサリと松尾玉を仕留めました。
将棋ソフトも最初、何が起こっているのか分からないまま、評価値が下がったのですが、その後、羽生さんの手順が進んでいくと、あっさりと手のひらを返して羽生さんが優勢の評価となりました。
将棋ソフト(AI)の読み筋をも超えた鬼神の将棋でした。
この将棋の内容なら、今の将棋ソフト(AI)にも勝つことができるんじゃないか?って思ってしまった羽生先生の将棋でした。
タイヤに亀裂が・・・ [自動車 / バイク]
タイヤ(後部)に亀裂がありました。
アウディではなく、スバルの方ですが。
かれこれ3年くらい乗っているので、そろそろ消耗する箇所も出てきてもおかしくない時期ですかね。
スバルに持って行ったところ、すぐに破裂する心配はないが、交換したほうが良いとのことでした。
しかも、バランスが悪くなるので、後ろは二本とも変えましょうと。
前はまだまだ大丈夫そうです。
後ろ二本の取り換えで7万ちょっと。
それなりの出費です。(涙
ただ、スバルの車は、高速道路を使って割と遠くまで旅行することが多いので、致し方なしというところです。
アウディではなく、スバルの方ですが。
かれこれ3年くらい乗っているので、そろそろ消耗する箇所も出てきてもおかしくない時期ですかね。
スバルに持って行ったところ、すぐに破裂する心配はないが、交換したほうが良いとのことでした。
しかも、バランスが悪くなるので、後ろは二本とも変えましょうと。
前はまだまだ大丈夫そうです。
後ろ二本の取り換えで7万ちょっと。
それなりの出費です。(涙
ただ、スバルの車は、高速道路を使って割と遠くまで旅行することが多いので、致し方なしというところです。
紫キャベツの実験 [ファミリー]
塾の理科の授業で「紫キャベツ 実験」をしました。
リトマス試験紙とかBTB溶液と同じ実験ですね。
酸・アルカリの指示薬として、ムラサキキャベツの色素(アントシアン)を使うことができます。
具体的には、どんな物質の酸性、アルカリ性を調べたかは不明ですが、本人は楽しかったようです。
「ムラサキキャベツを切るのが大変だった。」と言っていました。w
リトマス試験紙とかBTB溶液と同じ実験ですね。
酸・アルカリの指示薬として、ムラサキキャベツの色素(アントシアン)を使うことができます。
具体的には、どんな物質の酸性、アルカリ性を調べたかは不明ですが、本人は楽しかったようです。
「ムラサキキャベツを切るのが大変だった。」と言っていました。w
東京ディズニーシー ディズニー・イースター 2018 [Disney / ディズニー]
ディズニーのイベント「ディズニー・イースター」に行って来ました。
意外と混んでいませんでした。
当日の流れは、
ニモ&フレンズ・シーライダー
↓
ミラコスタでFTWランチ
↓
海底2万マイル
↓
スカットルのスクーター
↓
ワールプール
↓
フランダーのフライングフィッシュコースター
↓
インディ・ジョーンズ・アドベンチャー: クリスタルスカルの魔宮
でした。
映画『インクレディブル・ファミリー』が公開されるので、それに向けて、絶賛宣伝中!
アメリカンウォーターフロントにある、クラシックなヴィークル。
大きくなった今でも大好きです。
ミラコスタのレストランからショーを観戦しました。
一望できるけど、キャラとの距離が遠いのが少し残念です。
夕食は、行きつけの「ザンビーニ・ブラザーズ・リストランテ」へ。
コスパ良く美味しい料理(パスタ)を食べることができます。
意外と混んでいませんでした。
当日の流れは、
ニモ&フレンズ・シーライダー
↓
ミラコスタでFTWランチ
↓
海底2万マイル
↓
スカットルのスクーター
↓
ワールプール
↓
フランダーのフライングフィッシュコースター
↓
インディ・ジョーンズ・アドベンチャー: クリスタルスカルの魔宮
でした。
映画『インクレディブル・ファミリー』が公開されるので、それに向けて、絶賛宣伝中!
アメリカンウォーターフロントにある、クラシックなヴィークル。
大きくなった今でも大好きです。
ミラコスタのレストランからショーを観戦しました。
一望できるけど、キャラとの距離が遠いのが少し残念です。
夕食は、行きつけの「ザンビーニ・ブラザーズ・リストランテ」へ。
コスパ良く美味しい料理(パスタ)を食べることができます。
【FTWミラコスタ】ウェディング・アニバーサリー・ランチ 2018 その1 [Disney / ディズニー]
9月の挙式ですが、2017年のアニバーサリーが出来ていなかったので、遅ればせながら行って来ました。
前回は「地中海料理レストラン オチェーアノ」だったので、今回は「中国料理レストラン シルクロードガーデン」を利用しました。
最初席に着いたときは、ドナルドが置いていました。
ショーを観て、席に戻ってくると、ミッキーとミニーが置いていました。
右がミッキーで、左がミニー(リボン付き)だそうです。
キッズのコースメニュー。
コースメニューだと食べきれない場合は、ミッキー・ミニーチャーハンがあります。
ちなみに、下の写真は、リボンが付いているので、ミニーチャーハンになります。
前回は「地中海料理レストラン オチェーアノ」だったので、今回は「中国料理レストラン シルクロードガーデン」を利用しました。
最初席に着いたときは、ドナルドが置いていました。
ショーを観て、席に戻ってくると、ミッキーとミニーが置いていました。
右がミッキーで、左がミニー(リボン付き)だそうです。
キッズのコースメニュー。
コースメニューだと食べきれない場合は、ミッキー・ミニーチャーハンがあります。
ちなみに、下の写真は、リボンが付いているので、ミニーチャーハンになります。
【FTWミラコスタ】ウェディング・アニバーサリー・ランチ 2018 その2 [Disney / ディズニー]
FTWは9月に挙式なので、毎年秋にやっています。
そのため、秋っぽいコース料理なのですが、今回は春。
春らしいコース料理となっていました。
【大人用のメニュー】
【季節の前菜盛り合わせ】
一口サイズなので、とても食べやすいです。
お腹が空いているので、少々物足りない感じです。w
【蟹肉とホワイトアスパラガスのスープ】
何が入っているんだろうと思って開けると、、、
スープが入っていました。
【帆立貝とポテトのフリッター オニオンソース】
香ばしいです。
オニオンソースとの相性もばっちり。
【海老とイカと季節の野菜の炒め 西洋わさび風味】
西洋わさび、ですがツーンとする辛さではなく、子供でも食べれそう。
【皮付き豚バラ肉の紹興酒煮込み】
とても柔らかいお肉。
紹興酒の香りが口に広がるので、大人向け。
人によって好き嫌いが分かれる味でした。
【白身魚と筍の炒飯 アサリ入りクリームソース】
白身魚がとても美味しかった。
クリームソースがかかっていましたが、個人的にはソースがない方が良かったです。
少しチャレンジ感がしました。
【苺とヨーグルトのプリン シナモン風味
揚げ胡麻団子 フレッシュフルーツ添え】
この辺りになってくると、お腹がいっぱいですが、頑張って別腹と言いたい。。。
【FTW記念のアニバーサリーケーキ】
ダメ押しのケーキ。
ここでKOとなりました。。。
そのため、秋っぽいコース料理なのですが、今回は春。
春らしいコース料理となっていました。
【大人用のメニュー】
【季節の前菜盛り合わせ】
一口サイズなので、とても食べやすいです。
お腹が空いているので、少々物足りない感じです。w
【蟹肉とホワイトアスパラガスのスープ】
何が入っているんだろうと思って開けると、、、
スープが入っていました。
【帆立貝とポテトのフリッター オニオンソース】
香ばしいです。
オニオンソースとの相性もばっちり。
【海老とイカと季節の野菜の炒め 西洋わさび風味】
西洋わさび、ですがツーンとする辛さではなく、子供でも食べれそう。
【皮付き豚バラ肉の紹興酒煮込み】
とても柔らかいお肉。
紹興酒の香りが口に広がるので、大人向け。
人によって好き嫌いが分かれる味でした。
【白身魚と筍の炒飯 アサリ入りクリームソース】
白身魚がとても美味しかった。
クリームソースがかかっていましたが、個人的にはソースがない方が良かったです。
少しチャレンジ感がしました。
【苺とヨーグルトのプリン シナモン風味
揚げ胡麻団子 フレッシュフルーツ添え】
この辺りになってくると、お腹がいっぱいですが、頑張って別腹と言いたい。。。
【FTW記念のアニバーサリーケーキ】
ダメ押しのケーキ。
ここでKOとなりました。。。
最新の消費者行動から考えるマーケティングの新潮流 [マーケティング / 仕事]
清水先生のセミナーに行って来ました。
最新の消費者行動から考えるマーケティングの新潮流
https://smart-analytics.jp/events_and_seminars/mes-1/
理論的な話や、事例紹介をがっつり1日するというセミナーでした。
消費者行動論は、どちらかといえば、アンケートなどの小規模のデータが中心となった事例が多いかと思われます。
最近、ブームとなっているデータは、人の行動履歴(閲覧履歴や購買履歴)が中心となっており、これらのデータを使ったデータ解析(主に予測)が盛んにおこなわれています。
実際は、テキトーにそれっぽいロジックを作ってしまっている場合が多いのですが、
本来は、消費者行動論的なアプローチをロジックを組み込み
丁寧に統計モデルを作って行くべきだと思います。
今回のセミナーでいくつか面白そうなアイデアをもらったので、早速実践して行きたいです。
最新の消費者行動から考えるマーケティングの新潮流
https://smart-analytics.jp/events_and_seminars/mes-1/
理論的な話や、事例紹介をがっつり1日するというセミナーでした。
消費者行動論は、どちらかといえば、アンケートなどの小規模のデータが中心となった事例が多いかと思われます。
最近、ブームとなっているデータは、人の行動履歴(閲覧履歴や購買履歴)が中心となっており、これらのデータを使ったデータ解析(主に予測)が盛んにおこなわれています。
実際は、テキトーにそれっぽいロジックを作ってしまっている場合が多いのですが、
本来は、消費者行動論的なアプローチをロジックを組み込み
丁寧に統計モデルを作って行くべきだと思います。
今回のセミナーでいくつか面白そうなアイデアをもらったので、早速実践して行きたいです。
【超・獣神祭】天下の大泥棒 石川五右衛門 [ゲーム]
名前から想像すると、♂を想像していましたが、思いっきり♀キャラですか・・・。w
「ふれた敵を加速状態にする」という新SSが楽しみです。
超・獣神祭といえば、マナを2体持っているので、最近は、引かなかったんですが、
久々に引いてみたいです。
満員電車にベビーカー [時事 / ニュース]
【子連れでおでかけ】満員電車にベビーカーを使って乗るのはあり? なし?
http://ure.pia.co.jp/articles/-/73957
田舎だと、車を持っている家庭がほとんどだし、そもそも満員電車ってのもあまりない気がします。
朝の通勤時間帯のギスギスした時間帯でベビーカーを使うのは、赤ちゃんにとって危険も多いため、できるなら乗らない方が良い。
ただ、自分が自動車を買ったきっかけでもあるのだが、どの時間帯でも混んでいる電車ってあります。
朝も混んでいるけど、その他の時間もめちゃくちゃ混んでいる!
そして、お出かけするときに、ラッシュ時間を外したとしても、
逆に帰りのラッシュ時間にぶち当たってしまい、満員電車に乗る羽目になってしまうこともありました。
1人だけだと、折りたたんで抱っこする、という事もできるのですが、
2人目とかになってくると、赤ちゃんに必要な小物(オムツやら着替えやら)も増えてくるし、そういう荷物の持つ手段としてのベビーカーを使っている人も多いかと思われます。
電車での移動に限界を感じたために、自分は2人目の子供のタイミングで自動車を買いました。
ただ、東京の場合、自動車を持っていない人も多いだろうし、ラッシュ時にどうしても移動しなければならない場合もあるかと思います。
そういう場合は、温かい気持ちで見守っていきたいものです。
http://ure.pia.co.jp/articles/-/73957
田舎だと、車を持っている家庭がほとんどだし、そもそも満員電車ってのもあまりない気がします。
朝の通勤時間帯のギスギスした時間帯でベビーカーを使うのは、赤ちゃんにとって危険も多いため、できるなら乗らない方が良い。
ただ、自分が自動車を買ったきっかけでもあるのだが、どの時間帯でも混んでいる電車ってあります。
朝も混んでいるけど、その他の時間もめちゃくちゃ混んでいる!
そして、お出かけするときに、ラッシュ時間を外したとしても、
逆に帰りのラッシュ時間にぶち当たってしまい、満員電車に乗る羽目になってしまうこともありました。
1人だけだと、折りたたんで抱っこする、という事もできるのですが、
2人目とかになってくると、赤ちゃんに必要な小物(オムツやら着替えやら)も増えてくるし、そういう荷物の持つ手段としてのベビーカーを使っている人も多いかと思われます。
電車での移動に限界を感じたために、自分は2人目の子供のタイミングで自動車を買いました。
ただ、東京の場合、自動車を持っていない人も多いだろうし、ラッシュ時にどうしても移動しなければならない場合もあるかと思います。
そういう場合は、温かい気持ちで見守っていきたいものです。
輸入車は壊れやすい? [自動車 / バイク]
スバル(国産車)とAudi(輸入車)を乗っています。
Audiを購入するときに、営業の人に「Audi(輸入車)は、壊れやすいか」確認したところ、割と素直に語ってくれました。
輸入車は壊れやすい?
→ YES!
理由は2つあるようです。
理由 1
100個パーツを作ったとして、5個不良品が混じっているとする。
このあたりの製造工程(不良品の発生率)は、基本、同じ。
国産車は、100個検品して不良品を除去するようにしている。
輸入車は、あらかじめ不良品が入っているという認識の元、105個納品する。
理由 2
国産車は、パワステやクーラーなど本来100の能力がある場合、80くらいで制限をかけており、目いっぱいの能力を使わないようにしている。
輸入車は、100の能力がある場合、リミットをかけずに100の機能を出せるようにしている。
ということで、壊れやすいのは本当っぽい。
そのため、壊れた時にの保証も含めて、ちゃんとメンテナンスに入っていく必要がありそうです。
とうぜん、メンテナンス費用も有料なので、輸入車は金食い虫ではありますが。。。
Audiを購入するときに、営業の人に「Audi(輸入車)は、壊れやすいか」確認したところ、割と素直に語ってくれました。
輸入車は壊れやすい?
→ YES!
理由は2つあるようです。
理由 1
100個パーツを作ったとして、5個不良品が混じっているとする。
このあたりの製造工程(不良品の発生率)は、基本、同じ。
国産車は、100個検品して不良品を除去するようにしている。
輸入車は、あらかじめ不良品が入っているという認識の元、105個納品する。
理由 2
国産車は、パワステやクーラーなど本来100の能力がある場合、80くらいで制限をかけており、目いっぱいの能力を使わないようにしている。
輸入車は、100の能力がある場合、リミットをかけずに100の機能を出せるようにしている。
ということで、壊れやすいのは本当っぽい。
そのため、壊れた時にの保証も含めて、ちゃんとメンテナンスに入っていく必要がありそうです。
とうぜん、メンテナンス費用も有料なので、輸入車は金食い虫ではありますが。。。
小学校&幼稚園の運動会 2018 [ファミリー]
小学校と幼稚園の合同運動会がありました。
幼稚園だけの運動会は秋にあります。
快晴ということで、かなり暑い運動会でした。
3年生になり、1,2年生に比べて、直線の競争からトラックの競争に変わったり、
踊る演技もキレキレのダンスをしたりと、
1年前に比べて随分と変わったなぁという印象です。
娘の方も年長になり、しっかりと年少さんのサポートが出来ていました。
幼稚園だけの運動会は秋にあります。
快晴ということで、かなり暑い運動会でした。
3年生になり、1,2年生に比べて、直線の競争からトラックの競争に変わったり、
踊る演技もキレキレのダンスをしたりと、
1年前に比べて随分と変わったなぁという印象です。
娘の方も年長になり、しっかりと年少さんのサポートが出来ていました。
修士課程で取得した単位を博士課程で認定 [大学院(MBA)]
修士課程で取得した単位を博士課程の単位として認定してもらうことができました。
基本的なルールは、以下の通り。
1. 博士課程に繰り越しできる単位は4単位まで。
2. 修士の卒業単位は30単位必要ですが、4単位持っていくには34単位必要。
仮に、修士の単位を28、博士の単位を4単位、合計32取得していた場合、
博士課程に持っていける単位は、卒業に必要な30単位を除いた2単位のみ。
3. 修士課程の時の成績(A、A+)などは、全て、「認」という単位になる。
個人的には、持っていけるなら、認定しもらっておくと心の余裕が出てくる。
博士課程でも面白い授業がたくさんあるので、研究中心の学生生活を送りつつ、
負担にならない範囲で、博士の授業も受けておく、というスタンスが良いかと思われます。
基本的なルールは、以下の通り。
1. 博士課程に繰り越しできる単位は4単位まで。
2. 修士の卒業単位は30単位必要ですが、4単位持っていくには34単位必要。
仮に、修士の単位を28、博士の単位を4単位、合計32取得していた場合、
博士課程に持っていける単位は、卒業に必要な30単位を除いた2単位のみ。
3. 修士課程の時の成績(A、A+)などは、全て、「認」という単位になる。
個人的には、持っていけるなら、認定しもらっておくと心の余裕が出てくる。
博士課程でも面白い授業がたくさんあるので、研究中心の学生生活を送りつつ、
負担にならない範囲で、博士の授業も受けておく、というスタンスが良いかと思われます。
リメンバー・ミー [Disney / ディズニー]
カーズに続いて、ピクサーの長編映画『リメンバー・ミー』のDVDが発売。
実は、映画館で観たかったのですが、ガイコツが出てくるってことで息子が怖がってしまい観に行けませんでした。。。
実際は、怖いホラー映画というよりかは、感動系の傑作らしいです。
2018年7月18日に発売ということで、少し先になりますが、予約注文しておきました♪
実は、映画館で観たかったのですが、ガイコツが出てくるってことで息子が怖がってしまい観に行けませんでした。。。
実際は、怖いホラー映画というよりかは、感動系の傑作らしいです。
2018年7月18日に発売ということで、少し先になりますが、予約注文しておきました♪
リメンバー・ミー MovieNEX [ブルーレイ+DVD+デジタルコピー(クラウド対応)+MovieNEXワールド] [Blu-ray]
- 出版社/メーカー: ウォルト・ディズニー・ジャパン株式会社
- メディア: Blu-ray
MACにCommon Lispをインストール [大学院(MBA)]
大学院の授業でLISPの授業がありました。
LISPは、名前を聞いたことがある程度の言語なので、全くの初心者。
まずは、自分のMACにインストールします。
インストールは簡単で、
でインストールすることができます。
起動は
です。
簡単な演算。
分散の計算。
答えば、
2/5
と出てきます。
なんと、分数の計算もできるんですね。
2/5 + 1/5 = 3/5ということで、出力結果も
3/5
と出てきます。
LISPは、名前を聞いたことがある程度の言語なので、全くの初心者。
まずは、自分のMACにインストールします。
インストールは簡単で、
brew install clisp
でインストールすることができます。
起動は
clisp
です。
簡単な演算。
(+ 1 2)
分散の計算。
(/ 2 5 )
答えば、
2/5
と出てきます。
なんと、分数の計算もできるんですね。
(+ 2/5 1/5)
2/5 + 1/5 = 3/5ということで、出力結果も
3/5
と出てきます。
階層ベイズで作成したパラメータのシミュレーションについて その1 [階層ベイズ]
里村先生の『マーケティング・モデル』をベースに、階層ベイズのモデルを作り、シミュレーションを行いました。
コードは、教科書に載っていますし、共立出版のホームページからもダウンロードできます。
以下、里村先生のコードをベースに修正したものを載せていますが、
一部、教科書そのままだと、上手く動かない部分があるので修正しています。
(オリジナルのコードは、上記を参照してください。)
まずは、初期設定から。
階層部分のパラメータ(Delta)は、
> Delta
[,1] [,2]
[1,] -2 1
[2,] -1 2
で与えています。
続いて、シミュレーションデータの発生とロジットモデルの尤度関数を定義しています。
単純な二項ロジットなので、尤度関数も簡単です。
推定方法ですが、メトロポリス・ヘイスティングで推定しています。
シミュレーション回数は1万回。
1万回のシミュレーションなので、数分待てば、計算が終了します。
続けて、出てきた結果を確認。
ここまでは、オリジナルコードとほとんど同じです。
続けて、階層ベイズで作ったモデルを予測で使う場合、どうなるかです。
マーケティング・モデル 第2版 (Rで学ぶデータサイエンス 13)
- 作者: 里村 卓也
- 出版社/メーカー: 共立出版
- 発売日: 2015/04/09
- メディア: 単行本
コードは、教科書に載っていますし、共立出版のホームページからもダウンロードできます。
以下、里村先生のコードをベースに修正したものを載せていますが、
一部、教科書そのままだと、上手く動かない部分があるので修正しています。
(オリジナルのコードは、上記を参照してください。)
まずは、初期設定から。
## 階層ベイズ2項ロジットモデルの推定
library(bayesm)
library(psych) ## パラメータの散布図で使用
set.seed(1234)
nvar <- 2 ## ロジットモデルの説明変数の数
hh <- 1000 ## 個人数
nobs <- 10 ## 個人あたりの選択回数
nz <- 2 ## 個人属性の説明変数の数
Z <- matrix(c(rep(1,hh),runif(hh,min=-1,max=1)),hh,nz)
Delta <- matrix(c(-2,-1,1,2),nz,nvar)
iota <- matrix(1,nvar,1)
Vbeta <- diag(nvar)+.5*iota%*%t(iota)
階層部分のパラメータ(Delta)は、
> Delta
[,1] [,2]
[1,] -2 1
[2,] -1 2
で与えています。
続いて、シミュレーションデータの発生とロジットモデルの尤度関数を定義しています。
単純な二項ロジットなので、尤度関数も簡単です。
## シミュレーションデータの発生
hhdata=NULL
for (i in 1:hh) {
beta <- t(Delta)%*%Z[i,]+as.vector(t(chol(Vbeta))%*%rnorm(nvar))
X <- matrix(runif(nobs*nvar),nobs,nvar)
prob <- exp(X%*%beta)/(1+exp(X%*%beta))
unif <- runif(nobs,0,1)
y <- ifelse(unif < prob,1,0)
hhdata[[i]] <- list(y=y,X=X,beta=beta)
}
## 2項ロジットモデルの対数尤度関数の定義
loglike <- function(y, X, beta) {
p1 <- exp(X %*% beta)/(1 + exp(X %*% beta))
ll <- y * log(p1) + (1 - y) * log(1 - p1)
sum(ll)
}
推定方法ですが、メトロポリス・ヘイスティングで推定しています。
シミュレーション回数は1万回。
## ベイズ推定のための設定
R<- 10000
sbeta<- 1.5
keep<- 1
nhh <- length(hhdata)
nz <- ncol(Z)
nvar <- ncol(X)
## 事前分布のパラメータ
nu<- nvar+3
V<- nu * diag(rep(1,nvar))
ADelta <- 0.01 * diag(nz)
Deltabar <- matrix(rep(0, nz * nvar), nz, nvar)
## サンプリング結果の保存スペースの作成
Vbetadraw <- matrix(double(floor(R/keep) * nvar * nvar), floor(R/keep), nvar * nvar)
betadraw <- array(double(floor(R/keep) * nhh * nvar), dim = c(nhh, nvar, floor(R/keep)))
Deltadraw <- matrix(double(floor(R/keep) * nvar * nz), floor(R/keep), nvar * nz)
## 棄却率と対数尤度
reject <- array(0, dim = c(R/keep))
llike <- array(0, dim = c(R/keep))
## 初期値の設定
oldbetas <- matrix(double(nhh * nvar), nhh, nvar)
oldVbeta <- diag(nvar)
oldVbetai <- diag(nvar)
oldDelta <- matrix(double(nvar * nz), nz, nvar)
betad <- array(0, dim = c(nvar))
betan <- array(0, dim = c(nvar))
## 階層ベイズ2項ロジットモデルによる推定
for (iter in 1:R) {
rej <- 0
logl <- 0
# sV <- sbeta * oldVbeta # 里村先生のオリジナル
sV <- sbeta * V
root <- t(chol(sV))
## MH法による個人別betaのサンプリング
for (i in 1:nhh) {
betad <- oldbetas[i, ]
# betan <- betad + root %*% rnorm(nvar) # 里村先生のオリジナル
betan <- betad + root %*% rnorm(nvar) * 0.1
lognew <- loglike(hhdata[[i]]$y, hhdata[[i]]$X, betan)
logold <- loglike(hhdata[[i]]$y, hhdata[[i]]$X, betad)
logknew <- -0.5 * (t(betan) - Z[i, ] %*% oldDelta) %*% oldVbetai %*% (betan - t(Z[i, ] %*% oldDelta))
logkold <- -0.5 * (t(betad) - Z[i, ] %*% oldDelta) %*% oldVbetai %*% (betad - t(Z[i, ] %*% oldDelta))
alpha <- exp(lognew + logknew - logold - logkold)
if (alpha == "NaN")
alpha = -1
u <- runif(1)
if (u < alpha) {
oldbetas[i, ] <- betan
logl <- logl + lognew
}
else {
logl <- logl + logold
rej <- rej + 1
}
}
## 多変量回帰によるDeltaのギブスサンプリング(bayesmのrmultiregを利用)
out <- rmultireg(oldbetas, Z, Deltabar, ADelta, nu, V)
oldDelta <- out$B
oldVbeta <- out$Sigma
oldVbetai <- solve(oldVbeta)
## 現在までのサンプリング数の表示
if ((iter%%100) == 0) {
cat("繰り返し数", iter, fill = TRUE)
}
## keep回毎にサンプリング結果を保存
mkeep <- iter/keep
if (iter%%keep == 0){
Deltadraw[mkeep, ] <- as.vector(oldDelta)
Vbetadraw[mkeep, ] <- as.vector(oldVbeta)
betadraw[, , mkeep] <- oldbetas
llike[mkeep] <- logl
reject[mkeep] <- rej/nhh
}
}
1万回のシミュレーションなので、数分待てば、計算が終了します。
続けて、出てきた結果を確認。
## log-likelihoodのプロット
plot(llike,type="l",xlab="Iterations/10",ylab="",main="Log Likelihood")
## rejectのプロット
plot(reject,type="l",xlab="Iterations/10",ylab="", main="Rejection Rate of MH Algorithm")
## 平均パラメータのプロット
Deltas<-as.vector(Delta)
matplot(Deltadraw,type="l",xlab="Iterations/10",ylab="", main="Draw of Delta")
## burn-in 期間
R0<-floor(R/2)+1
R1=R
## Deltaの統計値
apply(Deltadraw[R0:R1,],2,mean)
rbind(Deltas,apply(Deltadraw[R0:R1,],2,quantile,probs=c(0.025,0.5,0.975)))
ここまでは、オリジナルコードとほとんど同じです。
続けて、階層ベイズで作ったモデルを予測で使う場合、どうなるかです。