「マルコフ連鎖モンテカルロ法 (統計ライブラリー)」の事例を勉強がてら練習してみた。
興味があったのは5.8節の変曲点のある回帰分析。途中で回帰線が折れ曲がったらどうなるの、っていうことですね。
ご存じの通り、私の体重はある日を境に落ちなくなっているので、この分析法で「いつから落ちなくなったのか」を推定してみた。もちろんMCMCです。
結果は次の通り。
mean sd alpha <span class="synConstant">72.25000</span> <span class="synConstant">0.151700</span> beta<span class="synSpecial">[</span><span class="synConstant">1</span><span class="synSpecial">]</span> -<span class="synConstant">0.07194</span> <span class="synConstant">0.001648</span> beta<span class="synSpecial">[</span><span class="synConstant">2</span><span class="synSpecial">]</span> -<span class="synConstant">0.01412</span> <span class="synConstant">0.002922</span> sigma <span class="synConstant">0.13750</span> <span class="synConstant">0.015860</span> x.change <span class="synConstant">93.09000</span> <span class="synConstant">2.921000</span>
alphaと書いてあるのは切片みたいなもんです。正確には、体重が落ちなくなったXデーのポイント。Xデーは最後のx.changeにあるとおり、ほぼ93日目だとわかります。つまり、93日目に72.25kgだったのを境として、それまでは日に719g落ちていたのが、それ以後日に141gしか落ちなくなった、ということです。しっかり推定されると悲しいな（笑）
ところで、MCMCの精度等々から考えると、実はこのシミュレーションは十分収束してないのです。500万回もやったのにな（笑）ちなみに、500万回の計算に728秒＝12分ちょっと。まぁ耐えられるかな。
この本に載っている事例をちょくちょく実行してみるのだが、数値がピッタリ一緒になることは、まぁない。シミュレーションだしね。ということは、だいたいのところがあっていればいいや、収束指標もOKだし、という考え方になっていくのかな。それとも、全てにおいてサンプル数の少なさの問題だろうか。テキストに載っているのは、少ないものであればN=20ぐらいのデータなので、そこまでの精度を求めるのは無理なのかもしれない。ちなみに私の体重データは158サンプルあります。
興味がある人は、私の体重データも公開するので、お試しください。コードは朝倉書店からDLしてね。よい子のみんなは、本も買おうね。
データ
<span class="synType">list</span><span class="synSpecial">(</span> Y = c<span class="synSpecial">(</span><span class="synConstant">71.6</span><span class="synSpecial">,</span><span class="synConstant">71.4</span><span class="synSpecial">,</span><span class="synConstant">71</span><span class="synSpecial">,</span><span class="synConstant">71.2</span><span class="synSpecial">,</span><span class="synConstant">71.2</span><span class="synSpecial">,</span><span class="synConstant">70.9</span><span class="synSpecial">,</span><span class="synConstant">71.2</span><span class="synSpecial">,</span><span class="synConstant">71.3</span><span class="synSpecial">,</span><span class="synConstant">71.2</span><span class="synSpecial">,</span><span class="synConstant">71.6</span><span class="synSpecial">,</span><span class="synConstant">71.3</span><span class="synSpecial">,</span><span class="synConstant">71.5</span><span class="synSpecial">,</span> <span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">71.6</span><span class="synSpecial">,</span><span class="synConstant">71.5</span><span class="synSpecial">,</span><span class="synConstant">72</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span><span class="synConstant">71.1</span><span class="synSpecial">,</span><span class="synConstant">71.5</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span><span class="synConstant">72.4</span><span class="synSpecial">,</span><span class="synConstant">72.1</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">72</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">72.2</span><span class="synSpecial">,</span> <span class="synConstant">72.3</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span><span class="synConstant">71.5</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span><span class="synConstant">72</span><span class="synSpecial">,</span><span class="synConstant">71</span><span class="synSpecial">,</span><span class="synConstant">71.5</span><span class="synSpecial">,</span><span class="synConstant">71.4</span><span class="synSpecial">,</span><span class="synConstant">72.1</span><span class="synSpecial">,</span><span class="synConstant">71.3</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span> <span class="synConstant">72.3</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">71.9</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">72.1</span><span class="synSpecial">,</span><span class="synConstant">72.5</span><span class="synSpecial">,</span><span class="synConstant">73.1</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span><span class="synConstant">71.6</span><span class="synSpecial">,</span><span class="synConstant">71.9</span><span class="synSpecial">,</span><span class="synConstant">71.7</span><span class="synSpecial">,</span> <span class="synConstant">71.5</span><span class="synSpecial">,</span><span class="synConstant">72.1</span><span class="synSpecial">,</span><span class="synConstant">72.6</span><span class="synSpecial">,</span><span class="synConstant">71.9</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">72</span><span class="synSpecial">,</span><span class="synConstant">72</span><span class="synSpecial">,</span><span class="synConstant">72.3</span><span class="synSpecial">,</span><span class="synConstant">71.8</span><span class="synSpecial">,</span><span class="synConstant">72.8</span><span class="synSpecial">,</span><span class="synConstant">72.3</span><span class="synSpecial">,</span><span class="synConstant">72.6</span><span class="synSpecial">,</span><span class="synConstant">72.9</span><span class="synSpecial">,</span><span class="synConstant">72.4</span><span class="synSpecial">,</span> <span class="synConstant">72.9</span><span class="synSpecial">,</span><span class="synConstant">72.9</span><span class="synSpecial">,</span><span class="synConstant">72.1</span><span class="synSpecial">,</span><span class="synConstant">72.6</span><span class="synSpecial">,</span><span class="synConstant">72.4</span><span class="synSpecial">,</span><span class="synConstant">72.6</span><span class="synSpecial">,</span><span class="synConstant">72.2</span><span class="synSpecial">,</span><span class="synConstant">72</span><span class="synSpecial">,</span><span class="synConstant">72.8</span><span class="synSpecial">,</span><span class="synConstant">72.2</span><span class="synSpecial">,</span><span class="synConstant">72.7</span><span class="synSpecial">,</span><span class="synConstant">73.3</span><span class="synSpecial">,</span><span class="synConstant">73.7</span><span class="synSpecial">,</span><span class="synConstant">73.3</span><span class="synSpecial">,</span> <span class="synConstant">73.3</span><span class="synSpecial">,</span><span class="synConstant">73.5</span><span class="synSpecial">,</span><span class="synConstant">73.2</span><span class="synSpecial">,</span><span class="synConstant">73.6</span><span class="synSpecial">,</span><span class="synConstant">73.9</span><span class="synSpecial">,</span><span class="synConstant">74.2</span><span class="synSpecial">,</span><span class="synConstant">73.3</span><span class="synSpecial">,</span><span class="synConstant">73.7</span><span class="synSpecial">,</span><span class="synConstant">73.7</span><span class="synSpecial">,</span><span class="synConstant">73.9</span><span class="synSpecial">,</span><span class="synConstant">74</span><span class="synSpecial">,</span><span class="synConstant">74.3</span><span class="synSpecial">,</span><span class="synConstant">74.5</span><span class="synSpecial">,</span><span class="synConstant">74.5</span><span class="synSpecial">,</span><span class="synConstant">75</span><span class="synSpecial">,</span> <span class="synConstant">74</span><span class="synSpecial">,</span><span class="synConstant">74.6</span><span class="synSpecial">,</span><span class="synConstant">75.2</span><span class="synSpecial">,</span><span class="synConstant">75.3</span><span class="synSpecial">,</span><span class="synConstant">74.3</span><span class="synSpecial">,</span><span class="synConstant">75</span><span class="synSpecial">,</span><span class="synConstant">74.6</span><span class="synSpecial">,</span><span class="synConstant">74.8</span><span class="synSpecial">,</span><span class="synConstant">75.3</span><span class="synSpecial">,</span><span class="synConstant">75.6</span><span class="synSpecial">,</span><span class="synConstant">75.1</span><span class="synSpecial">,</span><span class="synConstant">75.4</span><span class="synSpecial">,</span><span class="synConstant">76</span><span class="synSpecial">,</span><span class="synConstant">76</span><span class="synSpecial">,</span><span class="synConstant">75.1</span><span class="synSpecial">,</span><span class="synConstant">75.9</span><span class="synSpecial">,</span> <span class="synConstant">75.8</span><span class="synSpecial">,</span><span class="synConstant">76.2</span><span class="synSpecial">,</span><span class="synConstant">76.1</span><span class="synSpecial">,</span><span class="synConstant">76.3</span><span class="synSpecial">,</span><span class="synConstant">76.5</span><span class="synSpecial">,</span><span class="synConstant">76.4</span><span class="synSpecial">,</span><span class="synConstant">76.4</span><span class="synSpecial">,</span><span class="synConstant">75.6</span><span class="synSpecial">,</span><span class="synConstant">76.2</span><span class="synSpecial">,</span><span class="synConstant">76.4</span><span class="synSpecial">,</span><span class="synConstant">76.3</span><span class="synSpecial">,</span><span class="synConstant">76</span><span class="synSpecial">,</span><span class="synConstant">76.4</span><span class="synSpecial">,</span><span class="synConstant">76.8</span><span class="synSpecial">,</span><span class="synConstant">76.8</span><span class="synSpecial">,</span> <span class="synConstant">76.8</span><span class="synSpecial">,</span><span class="synConstant">77.4</span><span class="synSpecial">,</span><span class="synConstant">77</span><span class="synSpecial">,</span><span class="synConstant">76.7</span><span class="synSpecial">,</span><span class="synConstant">77.3</span><span class="synSpecial">,</span><span class="synConstant">77.4</span><span class="synSpecial">,</span><span class="synConstant">77.2</span><span class="synSpecial">,</span><span class="synConstant">77.5</span><span class="synSpecial">,</span><span class="synConstant">77.5</span><span class="synSpecial">,</span><span class="synConstant">78</span><span class="synSpecial">,</span><span class="synConstant">77.8</span><span class="synSpecial">,</span><span class="synConstant">77.1</span><span class="synSpecial">,</span><span class="synConstant">77.2</span><span class="synSpecial">,</span><span class="synConstant">77.4</span><span class="synSpecial">,</span><span class="synConstant">77</span><span class="synSpecial">,</span><span class="synConstant">77.5</span><span class="synSpecial">,</span> <span class="synConstant">77.4</span><span class="synSpecial">,</span><span class="synConstant">77.6</span><span class="synSpecial">,</span><span class="synConstant">77.7</span><span class="synSpecial">,</span><span class="synConstant">77.9</span><span class="synSpecial">,</span><span class="synConstant">77.6</span><span class="synSpecial">,</span><span class="synConstant">78.1</span><span class="synSpecial">,</span><span class="synConstant">78.3</span><span class="synSpecial">,</span><span class="synConstant">78.1</span><span class="synSpecial">,</span><span class="synConstant">78.4</span><span class="synSpecial">,</span><span class="synConstant">78.3</span><span class="synSpecial">,</span><span class="synConstant">78.5</span><span class="synSpecial">,</span><span class="synConstant">78.5</span><span class="synSpecial">,</span><span class="synConstant">78.4</span><span class="synSpecial">,</span><span class="synConstant">79</span><span class="synSpecial">,</span><span class="synConstant">79.5</span><span class="synSpecial">),</span> X = c<span class="synSpecial">(</span><span class="synConstant">158</span><span class="synSpecial">,</span><span class="synConstant">157</span><span class="synSpecial">,</span><span class="synConstant">156</span><span class="synSpecial">,</span><span class="synConstant">155</span><span class="synSpecial">,</span><span class="synConstant">154</span><span class="synSpecial">,</span><span class="synConstant">153</span><span class="synSpecial">,</span><span class="synConstant">152</span><span class="synSpecial">,</span><span class="synConstant">151</span><span class="synSpecial">,</span><span class="synConstant">150</span><span class="synSpecial">,</span><span class="synConstant">149</span><span class="synSpecial">,</span><span class="synConstant">148</span><span class="synSpecial">,</span><span class="synConstant">147</span><span class="synSpecial">,</span><span class="synConstant">146</span><span class="synSpecial">,</span><span class="synConstant">145</span><span class="synSpecial">,</span><span class="synConstant">144</span><span class="synSpecial">,</span><span class="synConstant">143</span><span class="synSpecial">,</span> <span class="synConstant">142</span><span class="synSpecial">,</span><span class="synConstant">141</span><span class="synSpecial">,</span><span class="synConstant">140</span><span class="synSpecial">,</span><span class="synConstant">139</span><span class="synSpecial">,</span><span class="synConstant">138</span><span class="synSpecial">,</span><span class="synConstant">137</span><span class="synSpecial">,</span><span class="synConstant">136</span><span class="synSpecial">,</span><span class="synConstant">135</span><span class="synSpecial">,</span><span class="synConstant">134</span><span class="synSpecial">,</span><span class="synConstant">133</span><span class="synSpecial">,</span><span class="synConstant">132</span><span class="synSpecial">,</span><span class="synConstant">131</span><span class="synSpecial">,</span><span class="synConstant">130</span><span class="synSpecial">,</span><span class="synConstant">129</span><span class="synSpecial">,</span><span class="synConstant">128</span><span class="synSpecial">,</span><span class="synConstant">127</span><span class="synSpecial">,</span><span class="synConstant">126</span><span class="synSpecial">,</span><span class="synConstant">125</span><span class="synSpecial">,</span><span class="synConstant">124</span><span class="synSpecial">,</span> <span class="synConstant">123</span><span class="synSpecial">,</span><span class="synConstant">122</span><span class="synSpecial">,</span><span class="synConstant">121</span><span class="synSpecial">,</span><span class="synConstant">120</span><span class="synSpecial">,</span><span class="synConstant">119</span><span class="synSpecial">,</span><span class="synConstant">118</span><span class="synSpecial">,</span><span class="synConstant">117</span><span class="synSpecial">,</span><span class="synConstant">116</span><span class="synSpecial">,</span><span class="synConstant">115</span><span class="synSpecial">,</span><span class="synConstant">114</span><span class="synSpecial">,</span><span class="synConstant">113</span><span class="synSpecial">,</span><span class="synConstant">112</span><span class="synSpecial">,</span><span class="synConstant">111</span><span class="synSpecial">,</span><span class="synConstant">110</span><span class="synSpecial">,</span><span class="synConstant">109</span><span class="synSpecial">,</span><span class="synConstant">108</span><span class="synSpecial">,</span><span class="synConstant">107</span><span class="synSpecial">,</span><span class="synConstant">106</span><span class="synSpecial">,</span><span class="synConstant">105</span><span class="synSpecial">,</span><span class="synConstant">104</span><span class="synSpecial">,</span> <span class="synConstant">103</span><span class="synSpecial">,</span><span class="synConstant">102</span><span class="synSpecial">,</span><span class="synConstant">101</span><span class="synSpecial">,</span><span class="synConstant">100</span><span class="synSpecial">,</span><span class="synConstant">99</span><span class="synSpecial">,</span><span class="synConstant">98</span><span class="synSpecial">,</span><span class="synConstant">97</span><span class="synSpecial">,</span><span class="synConstant">96</span><span class="synSpecial">,</span><span class="synConstant">95</span><span class="synSpecial">,</span><span class="synConstant">94</span><span class="synSpecial">,</span><span class="synConstant">93</span><span class="synSpecial">,</span><span class="synConstant">92</span><span class="synSpecial">,</span><span class="synConstant">91</span><span class="synSpecial">,</span><span class="synConstant">90</span><span class="synSpecial">,</span><span class="synConstant">89</span><span class="synSpecial">,</span><span class="synConstant">88</span><span class="synSpecial">,</span><span class="synConstant">87</span><span class="synSpecial">,</span><span class="synConstant">86</span><span class="synSpecial">,</span><span class="synConstant">85</span><span class="synSpecial">,</span><span class="synConstant">84</span><span class="synSpecial">,</span><span class="synConstant">83</span><span class="synSpecial">,</span><span class="synConstant">82</span><span class="synSpecial">,</span><span class="synConstant">81</span><span class="synSpecial">,</span><span class="synConstant">80</span><span class="synSpecial">,</span><span class="synConstant">79</span><span class="synSpecial">,</span><span class="synConstant">78</span><span class="synSpecial">,</span> <span class="synConstant">77</span><span class="synSpecial">,</span><span class="synConstant">76</span><span class="synSpecial">,</span><span class="synConstant">75</span><span class="synSpecial">,</span><span class="synConstant">74</span><span class="synSpecial">,</span><span class="synConstant">73</span><span class="synSpecial">,</span><span class="synConstant">72</span><span class="synSpecial">,</span><span class="synConstant">71</span><span class="synSpecial">,</span><span class="synConstant">70</span><span class="synSpecial">,</span><span class="synConstant">69</span><span class="synSpecial">,</span><span class="synConstant">68</span><span class="synSpecial">,</span><span class="synConstant">67</span><span class="synSpecial">,</span><span class="synConstant">66</span><span class="synSpecial">,</span><span class="synConstant">65</span><span class="synSpecial">,</span><span class="synConstant">64</span><span class="synSpecial">,</span><span class="synConstant">63</span><span class="synSpecial">,</span><span class="synConstant">62</span><span class="synSpecial">,</span><span class="synConstant">61</span><span class="synSpecial">,</span><span class="synConstant">60</span><span class="synSpecial">,</span><span class="synConstant">59</span><span class="synSpecial">,</span><span class="synConstant">58</span><span class="synSpecial">,</span><span class="synConstant">57</span><span class="synSpecial">,</span><span class="synConstant">56</span><span class="synSpecial">,</span><span class="synConstant">55</span><span class="synSpecial">,</span><span class="synConstant">54</span><span class="synSpecial">,</span><span class="synConstant">53</span><span class="synSpecial">,</span><span class="synConstant">52</span><span class="synSpecial">,</span><span class="synConstant">51</span><span class="synSpecial">,</span><span class="synConstant">50</span><span class="synSpecial">,</span> <span class="synConstant">49</span><span class="synSpecial">,</span><span class="synConstant">48</span><span class="synSpecial">,</span><span class="synConstant">47</span><span class="synSpecial">,</span><span class="synConstant">46</span><span class="synSpecial">,</span><span class="synConstant">45</span><span class="synSpecial">,</span><span class="synConstant">44</span><span class="synSpecial">,</span><span class="synConstant">43</span><span class="synSpecial">,</span><span class="synConstant">42</span><span class="synSpecial">,</span><span class="synConstant">41</span><span class="synSpecial">,</span><span class="synConstant">40</span><span class="synSpecial">,</span><span class="synConstant">39</span><span class="synSpecial">,</span><span class="synConstant">38</span><span class="synSpecial">,</span><span class="synConstant">37</span><span class="synSpecial">,</span><span class="synConstant">36</span><span class="synSpecial">,</span><span class="synConstant">35</span><span class="synSpecial">,</span><span class="synConstant">34</span><span class="synSpecial">,</span><span class="synConstant">33</span><span class="synSpecial">,</span><span class="synConstant">32</span><span class="synSpecial">,</span><span class="synConstant">31</span><span class="synSpecial">,</span><span class="synConstant">30</span><span class="synSpecial">,</span><span class="synConstant">29</span><span class="synSpecial">,</span><span class="synConstant">28</span><span class="synSpecial">,</span><span class="synConstant">27</span><span class="synSpecial">,</span><span class="synConstant">26</span><span class="synSpecial">,</span><span class="synConstant">25</span><span class="synSpecial">,</span><span class="synConstant">24</span><span class="synSpecial">,</span><span class="synConstant">23</span><span class="synSpecial">,</span><span class="synConstant">22</span><span class="synSpecial">,</span> <span class="synConstant">21</span><span class="synSpecial">,</span><span class="synConstant">20</span><span class="synSpecial">,</span><span class="synConstant">19</span><span class="synSpecial">,</span><span class="synConstant">18</span><span class="synSpecial">,</span><span class="synConstant">17</span><span class="synSpecial">,</span><span class="synConstant">16</span><span class="synSpecial">,</span><span class="synConstant">15</span><span class="synSpecial">,</span><span class="synConstant">14</span><span class="synSpecial">,</span><span class="synConstant">13</span><span class="synSpecial">,</span><span class="synConstant">12</span><span class="synSpecial">,</span><span class="synConstant">11</span><span class="synSpecial">,</span><span class="synConstant">10</span><span class="synSpecial">,</span><span class="synConstant">9</span><span class="synSpecial">,</span><span class="synConstant">8</span><span class="synSpecial">,</span><span class="synConstant">7</span><span class="synSpecial">,</span><span class="synConstant">6</span><span class="synSpecial">,</span><span class="synConstant">5</span><span class="synSpecial">,</span><span class="synConstant">4</span><span class="synSpecial">,</span><span class="synConstant">3</span><span class="synSpecial">,</span><span class="synConstant">2</span><span class="synSpecial">,</span><span class="synConstant">1</span><span class="synSpecial">),</span> N = <span class="synConstant">158</span><span class="synSpecial">)</span>
ソースコード(model.txt)の一部は次のように変えましょう。
x.change ~ dunif<span class="synSpecial">(</span><span class="synConstant">50</span><span class="synSpecial">,</span><span class="synConstant">150</span><span class="synSpecial">)</span>