Stanによる分布の近似
Stanによる乱数の近似；パラメータリカバリを通じて
乱数のタネを設定
# シード値を固定すると再現できる（擬似乱数なので）
set.seed(201808)データを作ってみます
これが今回の真値です。データを平均59，標準偏差9の正規分布から作っています。
Y <- round(rnorm(49, 53, 9))仮想データの中身
Y## [1] 56 44 51 51 53 51 41 67 62 60 40 41 63 39 54 57 65 61 44 66 43 45 55
## [24] 52 59 32 52 51 60 53 50 50 49 52 55 64 49 60 58 46 50 61 56 50 52 65
## [47] 62 60 45仮想データの平均値
mean(Y)## [1] 53.10204理論的真値は53ですが，確率変数の実現値は揺れるので53ジャストになることはありませんね。
我々はこのようなデータしか手に入らない中で，平均値はどこにあるのかを探ることになります。
stanモデルコンパイル
ここでR から Stanへ，「乱数発生器を作れ」との命令が行きます。 命令を受けたStanはコンパイルを始めます。少し時間がかかるのでお待ちください。
コンパイルが終わっても特に返事はして来ませんが，コンソール画面が「待ち」の状態になれば終わったと思ってください。
警告がたくさん出たとしても，それはコンパイル時の非常にテクニカルな指摘なので気にしなくて結構です。
# stan_model関数によるコンパイル
model1 <- stan_model("yoshinaya.stan")## data{
## int<lower=1> N;
## real X[N];
## }
##
## parameters{
## real mu;
## }
##
## model{
## X ~ normal(mu,9);
## mu ~ uniform(0,100);
## }stanに与えるデータセットを準備
# データセットはリスト型で渡す
dataset <- list(N=49, X=Y)stanをつかってサンプリング
# サンプリング関数でMCMCサンプルを作成(；´Д｀)ﾊｱﾊｱs
fit <- sampling(model1,dataset)結果の出力
print(fit)## Inference for Stan model: yoshinaya.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 53.11 0.03 1.27 50.62 52.25 53.12 53.98 55.65 1658 1
## lp__ -19.53 0.02 0.69 -21.57 -19.68 -19.27 -19.09 -19.04 1774 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:29 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).今回は真値として53の世界からデータが得られていると考えて，モデルを作り，推定値として53.11です。50% confidential inter
結果の図示
plot(fit,pars="mu",show_density=T)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
パッケージを使って美しく描きましょう
# ベイズプロットパッケージ
## stanfitオブジェクトをarray型にして渡す必要がある
bayesplot::mcmc_dens(as.array(fit),pars="mu")
bayesplot::mcmc_trace(as.array(fit),pars="mu")
bayesplot::mcmc_areas(as.array(fit,pars="mu"))
MCMCサンプルを取り出す
# rstanパッケージのextract関数
## extract関数は他のパッケージでも同名のものがある
rstan::extract(fit,pars="mu")## $mu
## [1] 52.74955 53.18230 53.39228 51.68559 53.72351 54.33388 50.68733
## [8] 53.27521 54.93270 53.35098 52.66755 53.09428 52.25171 53.27103
## [15] 53.02182 53.42924 52.71429 52.22592 53.44350 51.74899 54.33473
## [22] 53.37257 54.62530 56.66383 54.48915 52.84438 52.47471 51.40143
## [29] 53.17650 53.69437 51.39878 53.37823 52.41817 53.39685 53.43118
## [36] 52.22190 54.17689 53.56095 52.45640 53.97074 53.19910 50.96702
## [43] 54.64474 52.58582 52.54010 53.41467 54.29039 55.94868 52.28911
## [50] 53.51640 55.15224 52.47443 53.18230 54.99832 53.35342 52.98675
## [57] 52.45202 52.49212 54.59151 51.07078 52.23904 50.41275 54.59188
## [64] 55.80126 53.04215 54.10144 53.48933 53.60994 52.50286 54.62484
## [71] 53.26246 53.27172 53.69679 53.81879 53.01303 53.47297 52.83406
## [78] 52.88762 53.04763 51.37407 53.06059 51.73854 52.55881 51.42708
## [85] 52.31372 51.76312 53.11892 53.69694 50.63956 53.36775 52.76811
## [92] 52.46649 51.58681 51.53332 53.19213 50.93551 53.71960 53.82190
## [99] 52.77413 53.58318 55.79073 51.99463 54.45294 52.50873 52.00901
## [106] 53.17498 53.63792 54.17357 52.13422 51.95292 52.57878 53.02915
## [113] 51.89707 53.92403 52.59603 54.42947 49.72677 53.10339 55.16273
## [120] 52.01121 54.01287 51.39496 53.36114 52.68495 54.18482 54.44913
## [127] 51.69962 52.76277 52.19384 54.67245 54.76493 52.66803 53.10465
## [134] 53.78353 54.12675 53.04215 53.09560 52.71538 51.43857 51.47048
## [141] 51.97621 51.79664 53.48463 53.42074 55.00136 51.12993 54.24978
## [148] 53.31235 57.10370 52.76394 52.88821 54.27357 54.94349 52.65052
## [155] 52.44994 51.81375 54.01699 52.90130 54.35397 52.19955 53.82057
## [162] 53.44665 51.73854 54.10422 52.29316 52.66267 52.90982 52.05524
## [169] 54.10422 52.95148 52.38130 54.13102 52.16863 52.79317 55.21740
## [176] 53.55251 53.95083 53.15433 53.14765 52.16776 56.21343 52.84980
## [183] 51.35405 54.63458 55.22467 53.74101 52.02098 51.73137 53.60040
## [190] 53.25841 53.37409 51.09947 54.95847 50.74651 52.68871 53.03152
## [197] 55.23828 52.30672 53.86970 55.39824 54.03598 53.56168 53.76202
## [204] 53.61042 55.14148 54.25674 54.27357 53.18409 52.98261 52.01010
## [211] 53.25423 53.49479 54.51340 52.41884 51.65921 52.04041 51.70705
## [218] 53.14188 53.50736 55.16082 54.36130 54.61498 52.11979 51.30924
## [225] 52.54077 53.14188 51.95646 49.69358 52.63162 52.52367 52.16410
## [232] 51.28335 53.33392 52.94190 55.55813 54.72667 53.53091 53.00357
## [239] 53.59783 52.82984 52.33727 51.92577 53.23147 53.61015 51.80982
## [246] 55.24965 51.82454 54.64190 51.84370 53.88004 52.56147 53.62581
## [253] 54.55257 53.41234 53.57329 51.97143 54.28306 53.37110 52.35396
## [260] 55.56806 53.27521 53.49634 53.34376 51.61882 53.15907 53.06059
## [267] 54.87717 53.68984 54.75237 52.97501 53.20168 54.12129 52.28029
## [274] 53.15907 54.90216 53.22087 52.20997 53.35878 51.16187 52.73200
## [281] 53.04104 53.28835 53.54077 53.00470 55.27737 52.08452 51.02281
## [288] 52.26730 52.46349 52.73940 53.53091 53.41467 53.91806 54.44410
## [295] 52.45856 52.33013 53.42881 54.29466 53.36483 53.80258 54.87235
## [302] 56.18560 53.47902 51.65203 53.85640 54.60500 52.42302 51.69732
## [309] 50.88715 52.81118 53.16770 53.23312 52.21672 52.23904 54.80001
## [316] 54.16109 51.74389 56.20214 50.81482 52.41716 53.37881 53.98184
## [323] 52.88645 51.71119 54.70907 53.67619 55.44958 53.80635 53.58328
## [330] 53.64197 52.64310 54.80977 50.61382 54.46874 53.16738 54.52653
## [337] 54.26021 54.57574 52.40727 52.38291 54.85975 53.85531 54.61928
## [344] 51.60049 53.34827 52.78137 53.61436 53.27811 53.80595 52.73868
## [351] 53.01720 52.82650 52.22855 53.03981 52.69247 50.76202 51.53453
## [358] 54.40547 53.27432 55.32553 54.88304 53.21104 50.66158 55.68353
## [365] 53.06847 52.96438 53.62911 53.95599 53.49298 53.00627 55.62361
## [372] 53.37021 54.70170 53.35656 51.73282 53.23509 52.58802 53.10912
## [379] 52.47745 55.80165 52.07285 53.27530 53.08876 54.11950 53.20823
## [386] 54.12066 49.85265 53.76661 51.98009 51.67114 52.14595 54.49119
## [393] 53.25392 53.37881 53.16924 51.07078 54.60245 53.20846 53.01205
## [400] 54.07097 52.17270 51.17614 53.94573 53.47114 55.58770 54.39633
## [407] 53.11159 51.19517 54.55843 51.79025 52.26882 52.63351 53.93675
## [414] 53.64866 52.70275 54.51924 52.28306 51.81902 54.31129 52.47736
## [421] 53.31820 50.61932 51.42708 54.89168 54.09157 53.46285 53.36940
## [428] 53.26822 52.68470 50.85300 52.58152 52.26204 51.53911 52.83862
## [435] 52.05749 54.98498 53.08866 53.80849 52.49191 50.17923 53.71744
## [442] 53.93204 50.85378 53.26054 54.00752 52.43409 55.08140 53.20669
## [449] 53.93821 52.43357 51.41763 51.66488 50.16022 54.76820 51.62450
## [456] 50.29907 53.71017 52.61826 53.72443 51.17614 51.73854 52.74313
## [463] 53.97157 54.34427 53.66821 55.15668 53.20098 53.35624 52.95501
## [470] 53.47336 54.47851 53.71565 54.70828 54.78846 55.31856 52.16102
## [477] 51.19120 54.23724 54.38416 54.47294 52.77594 51.56686 52.48996
## [484] 51.77993 52.90115 51.69074 52.10242 53.95267 52.61359 51.88534
## [491] 53.03990 52.06871 53.00454 54.88694 53.92830 52.69606 51.95567
## [498] 52.78451 54.39428 54.98503 53.17636 53.59819 52.69120 53.22599
## [505] 55.29188 53.27464 51.59894 54.26197 54.18647 53.01233 49.81430
## [512] 53.56121 52.62478 53.95084 53.23106 54.30010 57.24455 53.41168
## [519] 55.28575 53.96878 52.84680 52.65520 52.84339 52.89037 53.19614
## [526] 53.28015 53.33013 53.09361 53.93650 54.66818 52.00913 53.42887
## [533] 53.43466 53.85293 52.26448 54.36397 52.89539 51.86729 52.43547
## [540] 53.05893 52.59494 54.04393 53.14866 52.55220 52.83732 53.38825
## [547] 50.62737 53.57220 53.79951 53.94527 52.58399 53.93903 53.16977
## [554] 54.49134 51.77867 53.91796 54.86152 50.81703 52.10242 52.68471
## [561] 53.39130 52.02627 53.19326 50.41291 52.37394 55.45317 52.56347
## [568] 52.99187 54.16290 51.51374 54.86629 54.87842 52.99326 53.67295
## [575] 54.04044 52.90206 54.93929 54.11941 54.71004 53.92626 55.96138
## [582] 53.39974 52.09167 54.69232 52.94600 53.11537 51.61503 53.06083
## [589] 53.24689 52.83838 52.53122 50.88620 55.41214 55.15387 55.28646
## [596] 53.24908 54.48681 53.50677 53.81644 52.82219 52.95376 53.70804
## [603] 51.20022 54.57568 53.75250 53.15238 53.49045 52.53017 52.58863
## [610] 52.39068 53.08279 52.06629 53.05992 51.68825 51.63025 52.35721
## [617] 53.14367 52.72884 53.92626 53.19301 49.64270 54.21124 54.06510
## [624] 51.97621 53.24728 53.68515 53.36660 50.93853 54.23360 56.19976
## [631] 51.44985 50.69352 52.24549 53.54714 53.94755 51.46017 54.55334
## [638] 54.51388 51.86963 55.98341 52.77604 53.46259 52.16047 53.49834
## [645] 52.33301 53.18409 52.73444 54.87740 53.85306 51.30595 53.26702
## [652] 52.59494 51.94569 51.34602 51.22831 53.08218 52.95148 51.95710
## [659] 54.87933 53.01720 51.89873 53.24084 53.89929 52.90759 53.54514
## [666] 52.43724 54.80846 52.42344 52.53508 51.87261 53.78838 52.17722
## [673] 52.36412 52.65224 55.15170 52.02468 51.53852 53.31045 53.24918
## [680] 54.53457 53.01163 55.95957 53.27432 55.63599 53.08218 54.94349
## [687] 55.13720 53.61517 50.48083 52.28360 53.91806 54.86914 54.54870
## [694] 52.80554 53.84714 51.43405 55.57375 54.94247 51.78303 51.17883
## [701] 53.77279 53.38641 52.10531 55.35608 54.05201 51.82361 54.61642
## [708] 54.47851 52.30292 53.41939 53.32484 51.70257 54.44470 54.03423
## [715] 51.47316 53.46825 54.16265 52.65579 54.51343 52.98462 52.55057
## [722] 51.58380 51.28616 52.93434 52.41512 52.86333 53.40900 52.20179
## [729] 51.58638 52.57189 54.17227 52.85427 51.50140 54.77457 51.31762
## [736] 54.31002 52.53887 53.72641 51.68819 51.87095 53.40539 50.74651
## [743] 52.29697 52.89341 54.36218 53.40593 52.89899 51.58689 53.83981
## [750] 51.97811 51.28616 53.88023 51.84895 54.19676 54.09463 54.63163
## [757] 53.16321 53.65222 53.04599 52.94625 51.89693 53.89929 51.91600
## [764] 53.03859 53.16485 53.54721 54.80977 57.74066 53.63109 52.76084
## [771] 51.32558 53.31294 53.81632 54.75491 52.82913 52.36381 52.39370
## [778] 55.02759 53.09809 52.02938 52.65520 55.31880 52.88876 53.00402
## [785] 53.80849 52.72980 51.21758 52.97114 54.39291 52.20525 52.29634
## [792] 53.36396 52.26511 53.81916 54.52130 54.33899 52.57746 51.10222
## [799] 53.02149 51.72246 50.53751 54.08371 50.74375 55.28033 53.19926
## [806] 53.27392 54.36218 51.69174 54.18564 53.81558 54.16109 51.99770
## [813] 50.19258 54.29910 52.51729 51.39180 51.71119 53.22694 52.78564
## [820] 51.02358 51.43064 54.21680 53.40564 53.34759 53.21291 50.43488
## [827] 54.41923 54.03742 54.62498 53.23509 53.55639 51.15658 53.11970
## [834] 50.83254 54.54469 53.20375 51.41172 54.61642 55.99040 55.65342
## [841] 52.05712 53.34715 53.65701 50.57546 52.43409 52.44956 51.21897
## [848] 49.98304 51.74190 51.92891 54.62498 53.47028 54.11788 52.42511
## [855] 53.33907 54.69705 52.72006 51.46058 52.38177 52.80554 53.39723
## [862] 52.29234 51.99801 53.08610 53.92626 53.10913 55.91401 53.43118
## [869] 55.83798 53.57354 54.29039 52.65393 54.69095 52.19886 52.57069
## [876] 53.49758 52.23221 53.11732 52.95950 52.22839 51.36814 53.26622
## [883] 52.52244 53.83745 52.26511 53.50680 52.18207 52.77190 53.11178
## [890] 52.66186 52.95327 51.91302 52.52434 54.01881 52.42252 52.27191
## [897] 56.67332 52.13756 53.46643 51.40966 52.72744 54.14101 52.57448
## [904] 53.76832 50.17824 53.46138 52.46838 54.00385 53.77054 54.80846
## [911] 53.67404 52.49212 52.05039 52.84359 51.16042 51.38709 52.20536
## [918] 54.69535 52.68370 54.18850 51.56823 55.14436 53.06033 53.93016
## [925] 54.44913 52.57876 53.55850 53.05301 54.35181 54.42947 52.37513
## [932] 53.30357 52.74021 53.55787 52.90289 51.48591 55.98497 52.74119
## [939] 53.34715 52.88876 53.30315 52.44665 53.49419 52.90906 53.29414
## [946] 53.30961 53.10580 52.80872 53.11202 51.47805 54.82978 53.96878
## [953] 53.53511 56.06447 55.12167 54.92808 52.43055 54.16591 51.21897
## [960] 52.84745 53.23804 53.45703 54.64083 49.59378 54.10675 53.84830
## [967] 52.89635 53.05098 53.47290 51.78592 50.81335 52.73149 52.63277
## [974] 50.28118 53.59678 53.96810 52.72666 53.72280 53.20098 53.48036
## [981] 52.16730 51.25171 51.66679 53.71325 52.71160 53.69898 52.62841
## [988] 50.80959 53.32022 50.17846 53.12005 52.80616 54.11409 55.50279
## [995] 53.55942 51.76756 54.80846 52.23780 52.10468 52.34907 52.78198
## [1002] 52.47261 53.55215 55.28506 54.93729 53.49822 56.16581 56.04998
## [1009] 53.52295 53.61550 51.85539 54.01705 52.42110 54.98049 52.21605
## [1016] 53.82357 54.95360 53.01332 54.61290 52.42572 54.31235 53.38052
## [1023] 50.40839 50.43446 52.52019 53.58815 53.73828 54.11443 53.82978
## [1030] 51.68638 53.75386 52.21087 52.28042 55.06905 53.29117 51.91302
## [1037] 50.68696 52.06115 50.08569 51.94209 52.98298 52.85254 53.63951
## [1044] 51.74802 52.48757 51.98472 53.46735 53.35146 54.27812 53.62699
## [1051] 51.42835 55.44864 51.30546 54.20887 54.04298 53.75200 53.21153
## [1058] 54.05926 53.53846 54.58596 52.10247 51.95030 53.05610 52.44854
## [1065] 53.97988 52.82857 53.62594 54.18084 51.66973 54.97478 52.57885
## [1072] 55.04727 52.15122 51.59098 52.17242 55.51460 53.84486 50.21094
## [1079] 52.08693 54.41785 55.18098 55.99171 51.60394 52.74099 53.56897
## [1086] 51.84186 54.35492 51.59830 52.89890 51.83675 53.51433 53.90090
## [1093] 53.25444 56.27863 55.03975 52.97026 51.05078 53.35009 54.54848
## [1100] 52.70373 54.84892 53.62650 52.34629 54.42105 52.81502 52.92440
## [1107] 50.83232 50.95673 52.09068 52.82161 51.21394 52.09317 54.12363
## [1114] 52.05323 52.56911 55.29911 51.56411 54.22379 53.92229 51.59237
## [1121] 51.65173 52.62744 55.35981 52.15317 54.23862 52.57588 50.79389
## [1128] 52.94355 53.72955 53.75080 51.30208 53.79834 52.43411 52.40300
## [1135] 53.41328 53.60408 53.06755 56.25541 53.03389 53.11917 52.73913
## [1142] 51.52755 53.82779 56.23849 54.27936 55.65624 51.67049 51.97567
## [1149] 52.64476 53.55552 54.55482 52.65058 53.25386 53.67418 52.17546
## [1156] 52.01449 54.52661 53.04603 55.15442 54.11411 52.48531 52.33604
## [1163] 53.09197 54.58492 52.58090 49.97230 52.96200 52.38348 50.41316
## [1170] 52.18689 53.14801 55.75782 50.75389 52.81393 55.60339 54.41419
## [1177] 51.78123 54.19225 54.12410 53.45020 52.18242 53.53193 55.47675
## [1184] 52.31507 51.54273 54.11832 53.02145 52.01648 53.60780 53.15743
## [1191] 52.45517 53.96546 52.76920 53.14089 53.41755 52.93552 52.30382
## [1198] 52.55446 52.93318 52.40674 53.10560 54.50530 55.69111 53.04077
## [1205] 55.35910 52.77451 52.29830 51.84974 53.93988 53.33977 53.62439
## [1212] 51.52632 50.84669 52.91299 50.65908 53.02593 52.59162 52.53252
## [1219] 52.56216 53.94661 56.02761 52.44373 53.64897 53.82903 54.22395
## [1226] 54.92628 55.35834 53.71593 55.21865 55.76944 53.02700 51.30677
## [1233] 52.86598 51.17353 54.30204 53.69315 53.30471 53.87817 52.95634
## [1240] 52.22568 54.39971 51.98202 55.84404 52.16199 52.18899 52.74871
## [1247] 51.55994 53.75498 52.77689 51.96822 53.45009 52.53066 50.93058
## [1254] 51.92991 55.73476 51.64879 54.08642 53.69399 52.21962 52.88573
## [1261] 52.87645 52.40486 50.52432 52.65886 52.93135 53.03521 53.21079
## [1268] 51.36484 53.61239 52.89324 51.28499 53.54322 54.60967 52.74235
## [1275] 53.19523 52.73769 52.53192 51.96119 53.75464 54.50307 52.95383
## [1282] 51.77725 50.92126 52.26330 53.70955 54.29025 52.14219 54.11382
## [1289] 51.79405 54.32590 50.46998 50.93294 51.70120 54.55449 54.13976
## [1296] 52.46577 51.57204 52.85376 55.04388 55.04457 53.43494 54.57821
## [1303] 55.38706 52.67248 52.46738 54.33753 52.92672 51.02288 54.38828
## [1310] 53.53941 54.19493 53.44977 52.75544 54.06429 53.51596 53.65509
## [1317] 53.15539 52.92884 52.64244 55.20139 54.70170 54.59051 54.20977
## [1324] 52.73797 52.88046 53.58477 54.70399 52.64581 51.28456 51.16070
## [1331] 51.71001 53.16280 53.21769 53.95631 53.13339 52.98958 54.64298
## [1338] 52.29799 52.04730 52.96487 54.17663 53.44977 51.40485 52.92433
## [1345] 53.67434 50.86395 51.62983 53.91073 54.84277 54.88755 54.62538
## [1352] 54.05217 53.04007 51.40907 54.32966 54.08759 51.94824 54.30204
## [1359] 54.03797 55.35834 53.35045 52.67594 53.19351 53.74500 54.38955
## [1366] 53.56817 52.20974 54.97807 52.14139 53.13579 52.49023 54.05501
## [1373] 53.03298 54.10287 54.48680 52.05053 54.87241 53.54466 51.67815
## [1380] 53.13988 51.83659 53.79041 54.48237 56.00198 53.56865 53.07562
## [1387] 52.99414 54.09628 54.08714 51.86560 51.05880 52.67912 52.06841
## [1394] 53.89765 54.14865 52.57314 54.35019 54.52468 54.59419 53.65011
## [1401] 49.89146 53.63672 52.82462 52.54470 53.26534 52.92196 53.99899
## [1408] 53.75031 52.41823 54.56036 52.71675 54.05189 54.13045 53.74769
## [1415] 51.96450 52.06100 51.76920 51.17807 52.78737 51.78535 54.73116
## [1422] 53.79302 53.43710 53.75765 51.87486 52.49162 52.90829 52.18460
## [1429] 53.53413 53.83694 52.55446 52.12966 54.95612 52.19534 52.55822
## [1436] 53.96758 53.79311 52.10546 52.54101 52.05055 52.51117 53.84944
## [1443] 49.68348 52.25097 54.57151 52.51113 49.78599 53.21063 53.01509
## [1450] 52.33731 50.23861 51.49585 55.00499 53.51493 52.62111 54.12355
## [1457] 54.99530 50.96296 52.95634 52.37145 54.25925 53.89397 53.17841
## [1464] 53.23370 51.36968 54.35735 54.45306 54.61064 51.03372 53.08008
## [1471] 53.39692 53.69166 54.21068 53.66154 54.73753 51.89039 54.80352
## [1478] 53.62956 52.26494 53.38935 53.85136 55.71867 51.97567 52.53052
## [1485] 52.62174 53.58877 53.54112 52.72464 50.96221 52.31973 52.90481
## [1492] 52.20221 54.93715 50.97603 51.80451 53.11971 53.14852 51.36000
## [1499] 52.99007 54.33580 51.01812 52.04723 52.68754 52.79391 53.79047
## [1506] 51.81205 53.31681 50.23375 53.30863 52.57145 54.35921 52.28694
## [1513] 51.94819 52.91399 54.43549 54.23256 51.02574 51.09479 51.60313
## [1520] 54.11412 54.04116 53.74393 49.93053 52.46405 55.13245 51.96644
## [1527] 50.67578 54.38195 52.09773 53.13001 55.96678 53.55288 53.66843
## [1534] 53.85069 53.03733 55.16070 54.21615 53.44281 51.89446 52.74099
## [1541] 51.18790 54.79774 54.10457 55.34846 54.61331 54.02669 52.06714
## [1548] 51.70013 52.04111 52.98855 53.75351 53.04913 53.28145 53.98182
## [1555] 53.63961 54.09552 55.64385 51.47029 53.84600 52.26785 52.77883
## [1562] 54.77905 52.86235 51.44616 53.82357 54.05501 52.51483 54.21196
## [1569] 53.78569 50.79448 53.40052 54.02494 53.33914 52.34846 54.60864
## [1576] 52.60886 51.17421 54.60967 52.72396 55.02061 53.84841 53.06987
## [1583] 51.36000 51.35432 52.82173 53.22932 54.17389 54.33544 52.79114
## [1590] 53.19570 51.59139 54.84029 52.52514 53.06191 53.78768 52.56994
## [1597] 54.04074 54.21476 52.51780 52.26917 53.33230 52.64065 52.64449
## [1604] 53.77675 53.35864 52.59201 54.58098 53.67173 53.67473 54.24453
## [1611] 54.08925 54.56848 52.39969 54.76358 52.68042 54.58169 53.04354
## [1618] 53.60933 52.92822 52.68931 53.14937 54.91413 55.49448 54.48680
## [1625] 52.28317 55.61816 53.12801 50.94963 53.68486 54.22489 52.34401
## [1632] 51.96348 51.33603 51.96629 53.94170 55.53443 52.08243 54.23586
## [1639] 53.88691 52.87596 52.84857 53.91081 54.36885 53.60555 53.74719
## [1646] 56.48951 53.08268 51.06436 51.96270 51.93290 52.66508 53.28877
## [1653] 55.54406 54.08735 53.32972 53.29832 54.14518 52.91379 51.97090
## [1660] 53.48573 52.99885 54.43077 51.42724 53.71440 53.08797 55.11940
## [1667] 53.33573 53.09577 51.87710 53.58484 54.49535 53.27617 55.63703
## [1674] 51.40313 53.16182 52.02697 52.27749 53.39289 53.42206 54.24508
## [1681] 55.37410 52.79152 50.05570 54.06321 55.08081 51.97453 53.08462
## [1688] 53.47959 52.50829 52.53726 52.64256 52.85617 53.72633 50.99828
## [1695] 55.18719 54.72002 53.88233 51.10522 53.52331 53.03389 50.77945
## [1702] 53.18870 51.84753 51.00524 52.49427 51.57777 53.61976 55.59324
## [1709] 52.52653 53.38539 53.65401 53.97220 52.36410 53.48037 51.84991
## [1716] 52.21883 53.36647 52.99250 54.31684 53.55644 52.37611 52.00069
## [1723] 51.07841 53.40588 51.57419 52.92481 51.39935 51.78148 52.74059
## [1730] 54.13672 54.28399 54.24866 54.22395 54.42311 54.78330 53.96034
## [1737] 55.53964 51.78943 52.09819 51.53523 51.62016 53.49822 54.04710
## [1744] 51.03708 55.31039 51.08832 53.00640 50.99222 53.13079 53.66561
## [1751] 49.48934 54.36506 53.34765 54.02791 54.77204 55.10100 55.34205
## [1758] 53.16619 52.31666 52.27076 50.29198 52.67690 53.57359 54.14470
## [1765] 52.80259 53.89110 52.51848 52.68109 53.14563 52.77847 52.09002
## [1772] 50.79754 52.27211 54.33485 53.62332 51.18210 53.50158 55.19051
## [1779] 54.67931 53.91129 51.58856 54.55771 54.48865 55.47230 52.63320
## [1786] 52.49211 52.96484 52.23241 53.37105 52.68428 53.82713 53.23024
## [1793] 53.63619 53.91081 54.15771 56.00796 53.12037 50.56672 53.20768
## [1800] 54.60801 54.48813 52.00611 54.22987 53.13914 53.88277 52.09010
## [1807] 52.31982 54.22424 56.92416 54.37997 52.76963 51.34726 52.49786
## [1814] 49.67024 52.39441 53.73184 54.60102 52.45448 53.79148 52.64853
## [1821] 51.69156 54.40326 53.78855 52.42703 52.77830 51.99601 54.05341
## [1828] 53.79317 55.10047 52.05494 52.42050 54.28046 52.40697 51.77618
## [1835] 54.17223 53.85069 52.42304 54.72914 52.36983 55.15165 51.30414
## [1842] 51.91192 55.28693 54.42493 52.60288 53.41856 52.98872 53.48787
## [1849] 53.58687 51.85410 54.13964 53.97349 52.80240 51.93439 53.20851
## [1856] 52.67451 52.60203 54.35735 50.82289 52.53699 53.54670 52.47180
## [1863] 51.45362 52.39489 54.41816 52.86841 51.76575 52.23136 55.83422
## [1870] 52.86087 52.58602 54.55148 52.01096 53.91692 51.59702 53.67298
## [1877] 54.80674 54.58188 52.37234 53.96909 51.88441 53.17526 53.49689
## [1884] 55.31994 54.65333 51.84390 51.48330 53.91135 56.12172 53.51188
## [1891] 52.96484 52.69975 52.19960 53.13579 54.09549 53.36627 53.17661
## [1898] 52.51390 54.41807 52.30375 53.46943 52.45517 53.32059 53.74963
## [1905] 53.78428 50.92748 55.35381 51.84810 53.45367 52.96919 53.94117
## [1912] 55.00695 54.65333 53.20846 50.91014 51.97422 55.81444 54.00554
## [1919] 51.69809 52.51874 52.72464 53.70736 52.57718 52.19764 53.56524
## [1926] 52.74938 54.55012 56.26931 52.84120 52.36381 51.92132 52.62991
## [1933] 54.00846 51.83952 53.53193 53.28660 51.54368 56.24048 53.51126
## [1940] 51.19569 52.65808 51.47340 51.15598 51.80451 53.31867 51.25650
## [1947] 54.56036 54.04842 54.70201 52.24162 52.39116 53.57854 53.00408
## [1954] 53.52951 54.13028 54.60585 53.69586 52.86770 52.90567 52.79370
## [1961] 53.24896 54.35878 52.02348 53.35624 52.04111 52.16733 53.53696
## [1968] 54.07692 50.14589 52.96570 54.22555 56.55602 53.19370 54.88726
## [1975] 52.35198 53.12801 52.61488 52.89509 51.76134 51.24250 53.00029
## [1982] 52.66437 53.22757 53.73109 53.36187 55.74907 52.80846 52.29899
## [1989] 51.57351 53.90357 53.41724 52.45737 53.79880 51.93496 53.60834
## [1996] 52.58255 51.98822 49.18271 53.01332 52.38258 53.35066 56.03703
## [2003] 51.89594 50.07661 55.91000 52.34063 54.96417 53.68483 56.17580
## [2010] 52.57822 51.98152 52.52997 52.54299 51.97950 54.84494 54.75885
## [2017] 53.30774 53.43826 54.90535 51.85552 52.65580 51.67532 52.64808
## [2024] 51.73736 53.91014 54.01888 55.03599 53.16112 52.75976 54.30922
## [2031] 52.47987 53.20461 53.40946 54.61341 54.31780 53.14795 52.59744
## [2038] 53.04233 51.77285 53.65014 51.60371 54.08190 53.08158 53.25269
## [2045] 54.03239 51.90653 52.13324 55.06979 52.84516 51.05287 54.82270
## [2052] 53.62106 53.36898 51.71471 54.02820 52.37449 54.82761 53.91610
## [2059] 52.60349 53.69939 54.65641 55.07323 54.04755 54.63749 53.57427
## [2066] 53.56512 51.43358 51.10781 52.40406 52.22299 51.84302 52.13989
## [2073] 52.56520 55.59270 51.89909 50.94182 53.59693 53.89963 53.41913
## [2080] 53.19123 54.04972 52.87112 52.80279 50.48329 52.15437 55.61896
## [2087] 52.01681 52.39624 53.14621 54.59669 51.73893 51.52979 51.15254
## [2094] 52.88633 55.13039 52.01403 53.81034 55.80458 52.96381 53.54002
## [2101] 51.74483 50.87567 54.02745 54.11440 51.33999 52.67641 53.81034
## [2108] 51.95239 52.51430 54.06499 53.99562 52.23735 54.03239 52.31174
## [2115] 52.95177 54.07886 51.60696 53.28170 51.95328 53.40827 52.54514
## [2122] 52.21539 54.26349 53.99525 54.21973 55.24220 53.09319 53.38936
## [2129] 51.87685 51.26777 53.18162 56.48578 53.19606 52.13537 55.68884
## [2136] 52.31631 52.83317 50.01154 55.31809 53.90780 54.37449 53.75656
## [2143] 52.26542 53.46180 53.09739 54.61643 50.76851 54.12138 54.92437
## [2150] 52.66287 51.57511 52.56389 52.40712 52.31621 53.50443 53.53147
## [2157] 54.02595 55.14052 53.17617 50.83289 53.32970 52.46676 53.81882
## [2164] 54.65931 54.11018 54.94767 53.42966 53.41171 53.04119 52.56426
## [2171] 51.25143 53.49925 54.70194 53.87554 52.55662 52.84478 53.84912
## [2178] 53.95204 54.01394 54.61341 53.43826 54.95543 52.57134 50.67816
## [2185] 53.48502 52.38378 53.24568 54.16587 53.43856 52.04632 51.56821
## [2192] 52.95539 50.90670 52.67992 53.25998 54.94229 52.16540 52.24306
## [2199] 53.19600 52.41435 51.41179 54.45723 51.45874 54.75061 53.24461
## [2206] 51.97431 54.44132 53.96114 53.99947 53.88039 54.16862 54.28620
## [2213] 52.79287 53.93610 54.13958 51.65504 53.40855 54.72356 51.42306
## [2220] 52.29049 53.92052 52.04000 51.85446 51.71159 53.24686 53.59693
## [2227] 54.15328 52.83135 52.76312 52.35598 55.96275 52.29514 54.13614
## [2234] 53.12537 52.62904 54.00300 52.24739 51.59190 52.82260 52.03139
## [2241] 53.24502 54.44772 52.86010 54.65829 52.61489 51.96298 54.80865
## [2248] 51.86088 54.23529 51.88473 51.87778 50.83892 52.42913 54.24295
## [2255] 52.30305 54.42931 55.49080 52.60034 50.96830 52.62878 51.87163
## [2262] 55.16394 54.28620 52.84516 54.65425 54.40228 52.42808 54.46290
## [2269] 49.99627 51.91120 53.03024 51.40919 53.02340 51.58713 50.78978
## [2276] 51.94075 52.36529 52.91125 54.92100 52.33442 53.28221 52.10879
## [2283] 55.06061 52.68681 53.83771 52.27699 51.73736 51.23507 53.01765
## [2290] 52.29247 54.87017 53.30763 54.94767 54.35393 54.12011 52.57919
## [2297] 52.82910 53.55783 53.65225 54.48875 52.01988 53.17259 54.05443
## [2304] 52.62831 51.88687 51.64445 54.29632 53.60898 53.47173 53.07655
## [2311] 53.25516 52.41557 53.22122 51.15056 55.13039 52.88815 53.73596
## [2318] 51.77808 54.45272 51.23798 52.24457 53.85487 52.74103 53.22228
## [2325] 55.07583 52.75175 51.47317 50.67647 56.95595 53.40267 55.33434
## [2332] 53.44134 52.03760 52.27875 53.04111 53.72925 51.52503 53.41086
## [2339] 52.70887 52.95715 52.64930 51.93480 53.47059 54.52484 54.37730
## [2346] 54.16587 54.27025 51.78275 55.37618 52.27779 54.46693 52.75899
## [2353] 53.07090 53.07113 52.78248 52.21127 53.87927 52.16236 51.25182
## [2360] 51.60873 54.19281 50.98508 52.58081 53.62139 50.43489 52.86221
## [2367] 53.00714 54.01394 51.33944 54.58269 53.30440 54.25617 52.95872
## [2374] 52.18743 52.64085 50.60620 53.29839 55.81456 53.54808 54.46836
## [2381] 52.41259 53.48182 53.99220 53.34586 53.17219 53.26744 54.27270
## [2388] 53.62734 54.41536 52.22299 53.68100 52.04247 52.79612 50.55787
## [2395] 50.69241 51.96176 54.28438 52.12579 49.99635 53.03464 52.35580
## [2402] 52.00248 51.69347 53.28315 53.27142 52.32616 54.01501 53.51316
## [2409] 54.30057 55.35166 55.42115 55.13854 55.37618 52.20625 52.72357
## [2416] 51.93795 55.02994 53.20314 54.07642 52.46243 54.77850 54.14502
## [2423] 51.78554 51.20381 54.76295 53.15096 54.01888 52.21200 51.59198
## [2430] 53.22861 52.99625 53.06828 53.82582 52.56194 51.52449 49.99924
## [2437] 52.56562 52.22509 54.48108 52.31738 54.77300 54.97288 51.59023
## [2444] 52.23076 53.19241 52.68553 50.67816 55.32880 53.33353 55.10992
## [2451] 54.87824 53.32694 50.60163 54.75963 52.57154 53.84820 52.25058
## [2458] 52.02164 54.59880 51.68546 54.33402 53.61138 52.17502 52.67100
## [2465] 53.37192 51.06793 55.33600 53.02846 51.61894 54.23224 52.09469
## [2472] 51.28658 52.77456 52.20535 56.30123 51.77498 53.81728 53.14980
## [2479] 52.80279 52.84478 51.16886 54.60283 50.95440 51.59272 52.21720
## [2486] 50.96143 54.43052 54.23020 55.42254 53.84025 52.39181 51.15503
## [2493] 52.36378 53.55729 54.50876 54.29256 54.17924 55.38152 52.51380
## [2500] 52.06821 55.11637 52.35409 54.60190 52.83279 52.94554 54.58916
## [2507] 54.35144 53.56454 52.51592 53.52320 53.80271 53.19598 55.43817
## [2514] 52.81216 53.87031 50.89063 53.23006 54.14718 53.71592 53.44121
## [2521] 53.46406 50.98508 52.94252 52.96884 51.97716 51.45484 52.83804
## [2528] 52.63294 51.91067 51.86167 52.28199 53.32635 53.02828 54.30315
## [2535] 52.35770 50.76118 51.02760 52.31312 52.37167 50.79570 52.14923
## [2542] 53.18605 55.18838 52.62573 52.72727 53.39193 53.69359 55.11102
## [2549] 51.87179 52.39696 53.29599 52.29049 53.73888 50.85890 53.28477
## [2556] 55.09325 54.64060 55.17680 51.23642 52.48578 52.29720 55.72871
## [2563] 54.13452 54.96328 50.68151 53.53457 52.26816 53.22076 52.20928
## [2570] 52.21720 52.60096 53.13811 50.80961 51.07078 51.68546 54.63840
## [2577] 51.15172 53.60898 54.74836 52.22883 52.72275 54.78787 53.42060
## [2584] 52.15399 52.24666 54.28028 52.39508 52.85756 52.11344 53.11963
## [2591] 54.50838 56.06825 51.48170 53.58320 52.10879 54.20288 53.58844
## [2598] 53.34845 54.20309 54.94984 56.44390 55.68393 54.07662 53.53806
## [2605] 52.10745 51.61884 53.53184 53.62543 52.36499 53.28138 49.67242
## [2612] 51.81201 52.93663 53.04387 51.68668 54.79489 54.77130 54.48691
## [2619] 52.95715 53.35152 53.45871 52.03979 51.90061 52.53194 53.60025
## [2626] 53.33148 51.77948 53.51134 51.42306 53.61076 53.92784 53.51338
## [2633] 53.33489 50.58473 52.39270 54.25170 52.82499 52.05734 51.55465
## [2640] 54.80702 53.01048 50.92846 54.44772 51.45638 53.72769 52.17502
## [2647] 52.21539 53.90446 53.98137 53.70675 54.19058 53.11124 54.42851
## [2654] 53.11124 53.44041 52.98086 50.90202 52.80718 53.91610 52.68868
## [2661] 53.28245 51.79964 53.74803 54.49071 55.19435 51.66916 54.40787
## [2668] 51.78275 52.16468 53.48504 53.35533 53.71200 55.56832 53.96114
## [2675] 54.20178 53.51307 52.36188 51.88725 56.41163 54.20831 54.31740
## [2682] 52.66325 53.98646 54.01459 54.55757 52.72880 54.14870 53.52204
## [2689] 53.00110 53.79350 52.62047 54.30300 54.07886 55.83988 50.53976
## [2696] 54.61981 54.15328 54.34819 52.60568 53.44858 51.71159 54.10481
## [2703] 53.51127 55.51389 52.80854 52.72757 52.30086 53.24167 53.89706
## [2710] 54.63346 53.28180 53.81544 52.80416 53.03935 53.25101 53.25045
## [2717] 55.16140 52.59933 53.43472 52.91705 50.92221 52.05850 53.84114
## [2724] 51.65550 52.25300 54.04638 53.89744 51.87928 50.60620 53.48658
## [2731] 53.10954 52.99621 52.77047 52.14869 53.56794 55.59433 50.70472
## [2738] 51.82100 53.29308 54.43393 53.18918 53.84228 50.37464 49.46259
## [2745] 51.14065 55.27989 53.66265 51.92863 54.42806 53.51999 50.24497
## [2752] 53.69683 50.86694 53.12133 53.50816 51.09746 51.64242 55.31036
## [2759] 52.88754 52.90519 53.89744 51.78489 52.91322 52.83212 52.60349
## [2766] 54.31774 53.33974 53.10954 52.65198 53.64142 52.73198 52.34294
## [2773] 51.55988 51.15503 54.04183 52.93862 53.87815 53.73888 51.56479
## [2780] 54.95126 53.53457 50.41951 53.19384 53.20485 55.10296 51.73604
## [2787] 53.66906 54.71668 52.16073 51.96345 51.92712 52.57871 52.88097
## [2794] 53.44700 53.05976 53.07743 54.03772 54.86417 54.04789 53.59013
## [2801] 51.23943 51.12048 54.70109 53.57232 52.23735 53.39357 51.74492
## [2808] 54.66534 51.14680 54.22618 52.98602 52.79287 51.67551 53.12007
## [2815] 51.91679 51.85750 51.96208 55.15902 54.27603 51.30618 53.56150
## [2822] 52.04842 54.01266 52.33860 54.28156 54.99640 52.97773 50.86138
## [2829] 54.12909 54.12341 53.23271 53.11214 52.53033 53.26882 54.01970
## [2836] 54.43103 54.42435 52.41674 54.29052 53.97850 55.87619 55.09584
## [2843] 52.91322 52.30086 53.27735 53.15283 52.26422 50.79631 51.92165
## [2850] 54.98314 50.33657 53.95155 53.97114 52.80294 54.45272 52.43603
## [2857] 55.32364 50.65974 52.27239 53.47303 54.91607 52.16024 52.58767
## [2864] 54.67091 54.74088 51.32377 53.34880 52.07993 52.90797 51.33670
## [2871] 54.55324 52.72604 54.23844 52.90797 56.41084 55.71936 50.73515
## [2878] 52.33879 53.19022 54.54458 53.74360 54.96934 52.19559 51.36235
## [2885] 54.56503 51.59301 53.29474 52.64285 54.38058 52.07869 54.71995
## [2892] 52.38159 50.45050 52.14252 53.75624 52.43748 53.81254 52.14129
## [2899] 51.58840 52.24306 54.79489 51.28881 52.16875 51.36706 50.94878
## [2906] 52.68887 53.28170 53.95488 53.03464 54.14986 53.50746 51.87837
## [2913] 52.40748 51.45197 54.77300 52.40029 54.11191 53.17208 51.25080
## [2920] 52.18606 53.74519 51.29426 52.38931 52.59878 54.29256 54.05750
## [2927] 55.04247 53.02362 56.14233 54.79782 52.66476 51.78064 51.44396
## [2934] 53.46961 52.93324 51.80877 52.11228 51.30618 53.34363 52.27259
## [2941] 55.88322 51.70078 52.35770 53.70695 52.89481 52.47987 53.92236
## [2948] 55.06480 53.25965 53.20847 53.89105 52.48933 50.98508 53.98863
## [2955] 53.66635 51.13468 52.90437 53.26311 53.67544 52.61742 52.07825
## [2962] 53.04343 52.91302 51.42257 52.90022 54.22688 52.35770 52.66287
## [2969] 52.75972 54.46476 53.17664 53.32943 51.96345 55.27939 52.72289
## [2976] 52.36440 56.78305 54.04789 51.97227 54.56213 53.39667 53.33148
## [2983] 51.75552 53.92236 54.16505 55.26999 51.01321 52.15445 53.29599
## [2990] 52.01542 50.12030 53.59986 53.34353 54.38882 55.07323 51.90816
## [2997] 51.57838 52.66247 53.22183 52.40014 53.91549 55.55227 52.35196
## [3004] 53.22187 52.92116 52.88503 52.65358 54.57560 50.91460 53.13501
## [3011] 51.53284 48.65735 52.22969 51.27617 51.05641 53.53072 53.81943
## [3018] 51.64008 52.08850 52.27983 51.13436 52.23750 52.92225 53.03148
## [3025] 53.94283 53.47234 54.28570 53.39045 54.71434 51.17453 55.91418
## [3032] 53.29726 48.65735 52.58494 50.85896 51.27302 52.57666 51.20971
## [3039] 52.34453 51.21886 50.80629 53.07753 54.17980 53.19578 51.27967
## [3046] 56.16024 52.61916 53.42111 53.59713 51.33771 53.41072 53.51863
## [3053] 53.30063 52.06084 53.85693 51.94604 54.20955 52.88045 52.60241
## [3060] 52.81001 52.87675 53.12359 51.73512 53.01921 52.88740 52.14750
## [3067] 51.93573 53.19765 51.54021 53.98002 53.83975 54.07744 53.33220
## [3074] 54.96740 55.24450 51.32148 53.39447 51.54988 53.50389 52.68384
## [3081] 51.48209 53.45470 52.60717 52.92184 53.48924 51.52942 51.01907
## [3088] 52.10406 51.89136 53.25233 52.89481 50.70468 51.28059 54.36598
## [3095] 52.54350 52.79122 53.88858 52.64365 53.75307 54.34900 52.21555
## [3102] 55.78659 52.24366 53.04389 52.94896 52.07585 51.84061 54.01548
## [3109] 52.34434 51.87918 53.32226 54.56258 53.29688 55.42118 52.94158
## [3116] 52.44694 50.79092 53.00636 53.92123 52.70508 53.05759 50.18797
## [3123] 54.41859 53.56413 52.17702 51.79769 51.27552 54.27650 52.39240
## [3130] 52.40548 54.58544 50.93214 54.19266 54.07313 53.11712 54.71317
## [3137] 55.29692 52.92763 50.64770 52.25066 51.86056 53.39395 53.75019
## [3144] 52.62074 52.75989 53.25114 53.62340 53.07310 54.10513 52.76641
## [3151] 53.29246 53.05709 53.05709 55.64668 54.45473 54.36328 53.42844
## [3158] 51.77000 53.86755 52.70413 54.49894 53.05025 51.85387 52.77389
## [3165] 51.74155 50.98168 54.57925 53.51373 54.81360 53.33765 51.96354
## [3172] 54.25850 52.62061 53.30912 51.98267 54.28536 54.26559 50.95730
## [3179] 54.52735 53.42233 53.44716 52.65353 52.63327 54.36207 53.83382
## [3186] 54.38935 53.37305 52.67442 53.22687 51.63024 54.33361 54.47382
## [3193] 52.87377 53.56246 51.16529 53.53872 51.03196 52.75963 52.83052
## [3200] 54.94066 51.80145 53.89587 53.89709 52.37144 53.54214 53.30912
## [3207] 51.31900 50.17533 53.01202 54.52509 52.93831 52.34710 53.14384
## [3214] 51.79008 52.42814 53.68169 51.84512 51.78984 54.50707 53.72604
## [3221] 54.65278 53.42558 53.29130 55.09984 52.83327 52.76869 53.21397
## [3228] 52.66660 52.72832 52.84193 53.92066 51.85068 52.37668 50.40183
## [3235] 53.74697 53.44406 52.19145 52.05139 51.33771 52.87913 53.92454
## [3242] 50.10329 52.87560 53.39203 53.21656 52.06776 53.39503 52.03123
## [3249] 54.15900 52.69543 52.36694 54.03233 53.22164 53.24510 53.85824
## [3256] 52.08413 53.22030 55.12477 54.01069 53.66902 54.80824 53.95742
## [3263] 52.62477 52.85111 54.26507 51.35183 51.83816 52.93528 52.47485
## [3270] 52.53061 52.05261 53.29547 51.79758 55.44792 53.12943 54.63845
## [3277] 54.51902 53.21419 52.43194 52.01369 50.80629 54.44973 53.09927
## [3284] 53.71540 52.61888 53.27670 52.36877 52.04983 53.45892 53.25302
## [3291] 53.55837 53.59030 53.99167 54.73105 52.12422 54.00769 50.27449
## [3298] 54.39412 54.61577 53.17876 54.01387 56.31274 51.51384 53.29071
## [3305] 52.11377 54.92821 53.45105 51.67437 51.76034 52.22466 52.22698
## [3312] 52.85489 53.58513 54.70794 53.43778 48.63836 53.40829 49.65669
## [3319] 53.36543 51.03747 52.29852 52.87141 52.67203 51.50939 53.69462
## [3326] 53.56707 52.48562 53.53997 54.32698 52.28156 52.05720 52.42322
## [3333] 54.24059 50.08710 52.22276 51.30133 52.72449 52.79912 52.28923
## [3340] 54.65587 56.02279 51.43960 55.85760 52.64299 52.39883 53.92854
## [3347] 55.90972 53.11003 53.85279 53.29974 53.21419 53.08274 56.19395
## [3354] 52.34097 51.18203 53.10816 49.77236 53.63213 55.22950 53.93548
## [3361] 54.27939 53.14015 51.34287 54.01085 55.88996 53.11291 54.06044
## [3368] 54.19592 53.03717 52.07255 52.21465 52.40122 55.23675 53.18102
## [3375] 54.45099 55.73717 52.13327 54.24229 53.09744 55.39970 52.14482
## [3382] 54.84435 53.31026 53.59148 50.03319 54.76068 53.64407 54.04258
## [3389] 54.27303 52.04436 53.10649 53.08631 56.18355 52.19887 53.50765
## [3396] 51.76599 54.35041 54.01672 54.12807 54.62279 53.08891 53.80632
## [3403] 53.56419 52.20077 54.33180 52.57616 52.91615 55.38728 54.62808
## [3410] 53.46254 54.23427 55.35151 50.91208 52.92928 52.50677 53.00588
## [3417] 52.73885 53.40848 53.20737 53.44675 53.59184 51.62753 52.36125
## [3424] 53.15326 52.58901 52.69303 54.52994 53.43494 54.68514 53.67193
## [3431] 53.77349 50.44983 55.04866 52.09309 55.74502 52.89075 54.05197
## [3438] 53.62701 53.50086 53.53556 53.42420 52.94559 54.48243 51.52012
## [3445] 53.27864 53.75065 52.78464 54.55623 53.33444 53.76296 52.67401
## [3452] 53.66503 53.83709 54.68747 53.46807 50.98539 53.38569 53.78807
## [3459] 54.21824 54.61007 54.80134 54.59621 54.52412 54.36238 49.68598
## [3466] 55.02170 53.37593 55.60880 53.75629 54.91510 51.62358 55.18966
## [3473] 52.71234 50.19354 53.92817 51.91144 52.11787 53.24891 54.31896
## [3480] 52.65056 52.94413 54.05316 53.94639 53.05394 53.94704 54.47382
## [3487] 56.06564 52.82526 54.40122 53.36137 57.10981 53.43334 54.43061
## [3494] 54.26844 52.61665 51.63988 51.88662 54.53050 51.29474 50.24918
## [3501] 53.33444 55.85030 53.67231 51.57836 52.46433 53.38787 51.94423
## [3508] 53.86799 51.10792 52.17702 52.95268 51.75385 51.85730 54.07602
## [3515] 52.95071 55.04770 54.17189 51.93752 50.59101 53.98684 51.58699
## [3522] 50.31823 52.35455 53.33873 51.85235 53.01760 51.81781 53.28322
## [3529] 54.51038 51.60379 54.74983 55.04599 50.39367 52.93941 49.97712
## [3536] 53.89597 54.14698 53.73811 54.51769 53.25959 52.79634 54.03351
## [3543] 52.81001 53.07453 52.50500 53.72667 52.82547 53.12815 52.05451
## [3550] 54.01809 52.34003 54.53614 53.08713 52.51230 50.74248 51.72373
## [3557] 53.98564 51.23506 51.89932 52.92969 55.82550 53.35862 54.44033
## [3564] 52.23975 52.92570 52.20137 55.80478 53.75684 50.49897 52.35151
## [3571] 54.17105 53.28852 53.80814 52.09457 52.92570 53.41745 50.93538
## [3578] 53.87550 50.71711 53.69628 52.88680 52.43119 51.54199 53.20608
## [3585] 53.86978 52.68649 54.16055 52.53699 54.85409 52.86052 51.14524
## [3592] 54.14228 53.11585 52.02359 52.13422 52.42354 52.81940 53.76214
## [3599] 53.62201 53.62144 51.44125 54.15311 51.40404 53.83667 54.67038
## [3606] 52.38191 52.89713 50.76374 51.92258 51.30927 52.50120 52.07210
## [3613] 51.75876 55.04670 54.17598 52.57563 52.08684 51.98618 53.05025
## [3620] 53.55298 53.18712 53.85846 53.63317 52.51034 53.97512 54.31096
## [3627] 52.56588 52.01027 51.41582 53.56403 52.66972 52.54603 50.94635
## [3634] 53.53412 54.20392 50.99531 54.45751 52.47180 54.03233 53.25126
## [3641] 54.20507 50.83720 54.15209 53.80884 53.30024 53.22778 54.04239
## [3648] 56.12591 52.19095 53.36710 54.07240 52.04544 55.12958 51.56726
## [3655] 52.09827 51.87852 54.44167 52.24056 53.52754 52.71697 49.36133
## [3662] 52.88299 51.23978 52.24425 50.54480 52.77909 53.30599 51.13318
## [3669] 53.07057 54.08109 56.06314 52.76221 52.48616 52.09982 54.90350
## [3676] 52.70082 52.70396 51.79007 51.94153 53.88296 53.50491 53.02246
## [3683] 50.64968 53.21271 54.09910 53.92028 53.36453 54.28193 55.06747
## [3690] 52.79161 55.09662 52.90721 52.67053 53.37821 51.57189 50.49117
## [3697] 51.86919 52.95281 51.83774 53.57636 53.52906 52.86432 55.43200
## [3704] 53.95742 53.09960 49.70417 53.33098 54.43175 54.52965 53.05176
## [3711] 50.78679 54.34699 53.28485 53.72620 51.85496 55.54116 53.36078
## [3718] 52.05620 52.59619 53.56361 52.94031 55.32792 51.82928 53.67050
## [3725] 54.21619 52.91277 54.15555 54.39330 53.23191 53.33082 54.51225
## [3732] 52.67403 54.60561 54.72699 53.56377 54.69483 50.74937 53.92568
## [3739] 52.12031 52.28004 52.28841 51.26517 52.12583 54.16892 53.14066
## [3746] 51.82299 52.58728 53.45905 52.35488 52.99172 52.86177 53.29755
## [3753] 55.23959 50.00576 55.10246 53.75570 53.05207 55.11586 51.99753
## [3760] 52.71897 51.17923 54.08961 51.22711 54.24761 54.72431 54.09092
## [3767] 53.15888 52.65242 50.28598 51.70687 54.10007 52.64009 51.57255
## [3774] 53.54166 51.69346 53.89174 53.99285 55.80391 54.34699 52.09187
## [3781] 54.49546 53.00557 54.93053 51.98062 51.05512 53.11952 53.90339
## [3788] 52.39283 54.15719 54.97767 55.80221 53.17317 53.03148 53.31969
## [3795] 51.54682 52.60687 54.48145 52.91652 52.36029 53.42420 52.81700
## [3802] 52.21323 52.80308 50.07555 53.12252 51.90623 53.95878 50.34277
## [3809] 53.08639 52.35187 53.36488 52.17024 53.77691 51.26137 55.12086
## [3816] 52.07240 53.15005 49.77702 49.24528 52.75399 52.59507 51.40765
## [3823] 52.04211 53.54206 53.16578 52.65466 53.36584 51.78189 52.25698
## [3830] 54.41741 50.89196 53.44380 54.11960 55.20961 52.66028 52.88764
## [3837] 52.78923 55.69497 53.16916 52.56614 52.06487 50.65899 52.40593
## [3844] 52.46987 53.84286 51.72630 54.13264 53.52545 53.79706 52.87720
## [3851] 54.66398 53.82636 53.18685 52.77513 52.62717 53.34092 52.15245
## [3858] 51.79208 52.88473 52.09155 51.48812 52.85489 51.42091 54.78899
## [3865] 54.97610 55.83222 54.11592 53.64186 55.17870 50.73447 52.17320
## [3872] 52.15572 53.37302 56.22898 53.09938 54.25348 52.96587 53.02175
## [3879] 54.36426 51.64611 55.98876 54.42713 52.50003 54.42491 55.21181
## [3886] 50.32545 53.37970 56.16624 53.99755 53.60896 54.54431 54.03453
## [3893] 54.48942 54.14401 53.17035 53.23191 52.29097 52.76556 50.46944
## [3900] 50.32593 52.75154 55.38256 52.24200 53.18712 51.83972 53.86866
## [3907] 53.72278 53.15721 53.18841 55.85353 52.62277 52.63127 53.65437
## [3914] 53.66824 52.03097 53.71976 53.51739 55.58599 50.95666 55.68059
## [3921] 50.75984 53.52427 54.21813 51.17733 52.10766 54.49546 55.17106
## [3928] 54.47462 53.83760 51.99530 53.26331 53.79543 53.25308 53.08958
## [3935] 49.36387 52.78957 54.80711 54.02620 52.84390 55.06929 51.61113
## [3942] 56.25507 52.84878 55.02285 49.77702 52.70099 52.08031 50.41786
## [3949] 52.72270 53.93489 51.89169 54.71136 52.85095 51.23147 53.49875
## [3956] 55.27246 54.17551 53.22030 52.11732 51.79248 51.24039 54.11283
## [3963] 53.22162 51.97138 51.63024 52.66874 54.61392 53.93271 52.43656
## [3970] 53.26114 51.84630 51.89748 55.23236 52.81032 52.31922 52.16474
## [3977] 53.85242 52.77851 53.43283 52.57088 52.87158 53.11572 52.94413
## [3984] 54.38524 54.03654 55.33059 52.54291 52.27850 53.08274 52.42118
## [3991] 52.69173 51.69638 52.48306 53.00787 51.48423 53.25951 53.81372
## [3998] 54.16703 51.77401 52.73958取り出したデータセットを見てみる
# パイプ演算子( %>% )は次の関数の第一引数に橋渡しする
rstan::extract(fit,pars="mu") %>% data.frame() %>% summary()## mu
## Min. :48.64
## 1st Qu.:52.25
## Median :53.12
## Mean :53.11
## 3rd Qu.:53.98
## Max. :57.74取り出したデータセットを使って
rstan::extract(fit,pars="mu") %>% data.frame() ->result.df
# データ数
NROW(result.df)## [1] 4000# 任意のパーセンタイル点
quantile(result.df$mu,probs = c(0.025,0.56,0.98))## 2.5% 56% 98%
## 50.61918 53.29570 55.80127# 任意の区間の密度
NROW(result.df[result.df$mu>55 & result.df$mu < 60,])/NROW(result.df)## [1] 0.0685# 別解
result.df %>% dplyr::filter(.$mu>55 & .$mu<60) %>% NROW(.)/4000## [1] 0.0685MCMCの様々な特徴
MCMCサンプルの数が増えると?
## 5000サンプル
fit5000 <- sampling(model1,dataset,warmup=2500,iter=5000,chains=4,thin=2)
## 10000サンプル
fit10000 <- sampling(model1,dataset,warmup=2500,iter=7500,chains=4,thin=2)
##
print(fit5000)## Inference for Stan model: yoshinaya.
## 4 chains, each with iter=5000; warmup=2500; thin=2;
## post-warmup draws per chain=1250, total post-warmup draws=5000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 53.14 0.02 1.29 50.68 52.24 53.15 54.03 55.64 3131 1
## lp__ -19.55 0.01 0.71 -21.49 -19.72 -19.28 -19.09 -19.04 3575 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:34 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).print(fit10000)## Inference for Stan model: yoshinaya.
## 4 chains, each with iter=7500; warmup=2500; thin=2;
## post-warmup draws per chain=2500, total post-warmup draws=10000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 53.07 0.02 1.29 50.54 52.20 53.08 53.93 55.58 6440 1
## lp__ -19.54 0.01 0.70 -21.52 -19.71 -19.27 -19.09 -19.04 7194 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:37 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).データ数が増えると?
# データを二つ用意。サンプルサイズ10のY1と100のY2
Y1 <- round(rnorm(10, 53, 9))
Y2 <- round(rnorm(100,53, 9))
# リスト型で渡す
dataset1 <- list(N=10,X=Y1)
dataset2 <- list(N=100,X=Y2)
# 同じモデルでサンプリング（データが異なる）
mcmc1 <- sampling(model1,dataset1)
mcmc2 <- sampling(model1,dataset2)
# 結果の表示
print(mcmc1)## Inference for Stan model: yoshinaya.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 52.02 0.07 2.82 46.55 50.14 52.00 53.87 57.43 1540 1
## lp__ -5.58 0.02 0.71 -7.49 -5.75 -5.31 -5.14 -5.09 1928 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:39 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).print(mcmc2)## Inference for Stan model: yoshinaya.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 53.63 0.02 0.90 51.91 53.00 53.63 54.24 55.39 1444 1
## lp__ -45.19 0.02 0.69 -47.17 -45.34 -44.92 -44.75 -44.69 1862 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:42 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).事前分布を活かす
# データセット全体
Y_all <- round(rnorm(30,53,9)) %>% print## [1] 51 51 67 52 65 47 48 48 65 48 46 52 55 52 61 45 36 62 51 46 63 61 56
## [24] 46 41 54 51 51 57 50# データセットを半分こ
Y_pre <- Y_all[1:15] %>% print## [1] 51 51 67 52 65 47 48 48 65 48 46 52 55 52 61Y_post <- Y_all[16:30] %>% print## [1] 45 36 62 51 46 63 61 56 46 41 54 51 51 57 50model2 <- stan_model("yoshinaya2.stan")## data{
## int<lower=1> N;
## real X[N];
## real pre_mu;
## real pre_sig;
## }
##
## parameters{
## real mu;
## }
##
## model{
## X ~ normal(mu,9);
## mu ~ normal(pre_mu,pre_sig);
## }# 全データを使って
dataAll <- list(N=30,X=Y_all,pre_mu=0,pre_sig=100)
fit_all <- sampling(model2,dataAll)
fit_all## Inference for Stan model: yoshinaya2.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 52.57 0.04 1.63 49.41 51.48 52.55 53.62 55.68 1585 1
## lp__ -10.54 0.02 0.73 -12.58 -10.71 -10.26 -10.10 -10.05 1883 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:44 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).# 前半分
data_pre <- list(N=15,X=Y_pre,pre_mu=0,pre_sig=100)
fit_pre <- sampling(model2,data_pre)
# 無情報後ろ半分
data_post <- list(N=15,X=Y_post,pre_mu=0,pre_sig=100)
fit_withNoPrior <- sampling(model2,data_post)
# 事前分布の情報を前半分データから取り出す
rstan::extract(fit_pre,pars="mu") %>% data.frame %>%
gather(key,val) %>% summarize(pre_mu=mean(val),pre_sig=sd(val)) -> Prior_val
# 事前分布の情報を追記
data_post2 <- list(N=15,X=Y_post,pre_mu=Prior_val$pre_mu,pre_sig=Prior_val$pre_sig)
fit_withPrior <- sampling(model2,data_post2)
# 各条件の出力を比較
print(fit_pre)## Inference for Stan model: yoshinaya2.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 53.72 0.06 2.34 49.16 52.15 53.69 55.33 58.29 1497 1
## lp__ -5.05 0.02 0.72 -7.05 -5.21 -4.77 -4.59 -4.54 1682 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:47 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).print(fit_withNoPrior)## Inference for Stan model: yoshinaya2.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 51.41 0.06 2.30 46.70 49.90 51.47 52.96 55.80 1571 1
## lp__ -5.84 0.02 0.71 -7.69 -5.99 -5.57 -5.40 -5.35 1878 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:49 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).print(fit_withPrior)## Inference for Stan model: yoshinaya2.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 52.56 0.04 1.68 49.36 51.41 52.52 53.69 55.86 1592 1
## lp__ -6.00 0.02 0.73 -8.10 -6.16 -5.71 -5.54 -5.48 1871 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:51 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).print(fit_all)## Inference for Stan model: yoshinaya2.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 52.57 0.04 1.63 49.41 51.48 52.55 53.62 55.68 1585 1
## lp__ -10.54 0.02 0.73 -12.58 -10.71 -10.26 -10.10 -10.05 1883 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:44 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).事後予測分布を描こう
モデルのチェックにも使える事後予測分布を描く練習です。
# 標準正規分布のデータを作成
N <- 50
Z <- rnorm(N,0,1)
model_norm <- stan_model("normal_PPD.stan")## data{
## int<lower=0> N;
## real X[N];
## }
##
## parameters{
## real mu;
## real<lower=0> sig;
## }
##
## model{
## // likelihood
## X ~ normal(mu,sig);
## //prior
## mu ~ normal(0,100);
## sig ~ cauchy(0,5);
## }
##
## generated quantities{
## real pred;
## pred = normal_rng(mu,sig);
## }# MCMC！
fit.norm <- sampling(model_norm,data=list(N=N,X=Z))
print(fit.norm)## Inference for Stan model: normal_PPD.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 0.01 0.00 0.15 -0.29 -0.10 0.01 0.11 0.32 2408 1
## sig 1.08 0.00 0.11 0.89 1.01 1.08 1.15 1.33 2789 1
## pred -0.01 0.02 1.10 -2.18 -0.74 0.01 0.72 2.13 3495 1
## lp__ -28.20 0.03 1.01 -30.96 -28.61 -27.88 -27.47 -27.21 1602 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 3 19:28:54 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).同時確率空間
今回のモデルは一つの正規分布ですが，パラメータが二つあります。\(\mu\)と\(\sigma\)の二つを同時に推定しました。それぞれの分布を描くと次のようになります。
fit.norm.df <- rstan::extract(fit.norm,pars=c("mu","sig")) %>% data.frame
ggplot(fit.norm.df,aes(x=mu))+geom_histogram(binwidth = 0.01)
ggplot(fit.norm.df,aes(x=sig))+geom_histogram(binwidth=0.01)
しかしこの二つは，同時に起こり得た乱数の組み合わせセットなのです。同時確率空間からのサンプリングであり，それぞれのパラメータのヒストグラムは「周辺化した」ものになっています。
# ポイントを有効にするためにplotlyパッケージを使う
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout# ggplotのオブジェクトをplotlyの関数に渡すだけ
gg <- rstan::extract(fit.norm,pars=c("mu","sig")) %>% data.frame %>%
ggplot(aes(x = mu, y = sig)) + geom_point()
ggplotly(gg)## 三次元表示
x_c <- rstan::extract(fit.norm,pars=c("mu")) %>% unlist() %>% cut(., 10)
y_c <- rstan::extract(fit.norm,pars=c("sig")) %>% unlist() %>% cut(., 10)
z <- table(x_c, y_c)
plot_ly(z=z,type='surface')これはパラメータの数が増えても同じことです。十数個のパラメータであっても，それらが同時に成立しうる組み合わせのセットである，という感覚を持つことは重要です（それぞれのEAP推定値が同時確率空間の中でのEAPになっていないことがある）。
また，周辺化を解析的に行うには積分が必要ですが，MCMCサンプルの場合は注目したい変数についての集計をするだけでよく，これがベイズ統計学を一気に実用的なものにした理由でもあります。
データの分布
ggplot(data.frame(val=Z),aes(x=val))+geom_histogram(binwidth = 0.1)
事後予測分布
rstan::extract(fit.norm,pars="pred") %>% data.frame %>%
ggplot(aes(x=pred))+geom_histogram(binwidth = 0.1)