2023-01-29

機械学習：Rでやってみる決定木

機械学習

機械学習をJuliaでやろうと思ったものの、新しい言語を覚えるのもそれなりに労力がいる。
最近ようやくRに慣れてきたところなので、これを定着させることを目指してRでやることにした。

機械学習の勉強に役立つサイトとして有名な「kaggle」を活用させてもらうことにして、以下のデータセットを使用して決定木を試してみる。

Titanic: Machine Learning from Disaster | Kaggle

まず、データの読み込み。

#csvデータ読み込み
train <- read.csv('../all/train.csv')
test <- read.csv('../all/test.csv')

データの中身を軽く確認してみる。

summary(train)

>  PassengerId       Survived          Pclass     
> Min.   :  1.0   Min.   :0.0000   Min.   :1.000  
> 1st Qu.:223.5   1st Qu.:0.0000   1st Qu.:2.000  
> Median :446.0   Median :0.0000   Median :3.000  
> Mean   :446.0   Mean   :0.3838   Mean   :2.309  
> 3rd Qu.:668.5   3rd Qu.:1.0000   3rd Qu.:3.000  
> Max.   :891.0   Max.   :1.0000   Max.   :3.000  
>                                                 
>                                    Name         Sex           Age       
> Abbing, Mr. Anthony                  :  1   female:314   Min.   : 0.42  
> Abbott, Mr. Rossmore Edward          :  1   male  :577   1st Qu.:20.12  
> Abbott, Mrs. Stanton (Rosa Hunt)     :  1                Median :28.00  
> Abelson, Mr. Samuel                  :  1                Mean   :29.70  
> Abelson, Mrs. Samuel (Hannah Wizosky):  1                3rd Qu.:38.00  
> Adahl, Mr. Mauritz Nils Martin       :  1                Max.   :80.00  
> (Other)                              :885                NA's   :177    
>     SibSp           Parch             Ticket         Fare       
> Min.   :0.000   Min.   :0.0000   1601    :  7   Min.   :  0.00  
> 1st Qu.:0.000   1st Qu.:0.0000   347082  :  7   1st Qu.:  7.91  
> Median :0.000   Median :0.0000   CA. 2343:  7   Median : 14.45  
> Mean   :0.523   Mean   :0.3816   3101295 :  6   Mean   : 32.20  
> 3rd Qu.:1.000   3rd Qu.:0.0000   347088  :  6   3rd Qu.: 31.00  
> Max.   :8.000   Max.   :6.0000   CA 2144 :  6   Max.   :512.33  
>                                  (Other) :852                   
>         Cabin     Embarked
>            :687    :  2   
> B96 B98    :  4   C:168   
> C23 C25 C27:  4   Q: 77   
> G6         :  4   S:644   
> C22 C26    :  3           
> D          :  3           
> (Other)    :186

年齢にNAがある模様。見にくいので欠損値の数のみ確認する。

#欠損値の数確認
is_na_train <- sapply(train, function(y) sum(is.na(y)))
is_na_train

>PassengerId    Survived      Pclass        Name         Sex         Age 
>          0           0           0           0           0         177 
>      SibSp       Parch      Ticket        Fare       Cabin    Embarked 
>          0           0           0           0           0           0

欠損値があるのはAgeだけなのかな？と思って、データを目視確認したところ、ところどころブランク（空白）がある。これはNA判定されていない模様。そこで、空白はNAにするようにデータ読み込みし直すことに。

#ブランクもNAにして読み込む
train <- read.csv('../all/train.csv', stringsAsFactors=F, na.strings = c('NA', ''))
test <- read.csv('../all/test.csv', stringsAsFactors=F, na.strings = c('NA', ''))

「stringsAsFactors=F」というのは「文字列型をFactor型にしないように」というオプション。

再び欠損値数の確認。

#欠損値の数確認（再）
is_na_train <- sapply(train, function(y) sum(is.na(y)))
is_na_train

>PassengerId    Survived      Pclass        Name         Sex         Age 
>          0           0           0           0           0         177 
>      SibSp       Parch      Ticket        Fare       Cabin    Embarked 
>          0           0           0           0         687           2

Age以外にも欠損値がありました。Cabinにかなりの数の欠損が。でも、モデルに使わないので無視。Embarkedの欠損が2件。
テストデータの方も確認しておく。

is_na_test <-  sapply(test, function(y) sum(is.na(y)))
is_na_test

>PassengerId      Pclass        Name         Sex         Age       SibSp 
>          0           0           0           0          86           0 
>      Parch      Ticket        Fare       Cabin    Embarked 
>          0           0           1         327           0

モデルを作る前にデータのクリーニングをしなくてはならない。まず欠損値を埋める処理を行う。

#Ageの欠損値を埋める。とりあえずmedianを入れる
train$Age[is.na(train$Age)] <- median(train$Age, na.rm=T)
test$Age[is.na(test$Age)] <- median(test$Age, na.rm=T)

#Embarkedの欠損値を埋める。とりあえず数が最も多いSを入れる
#testには欠損値なし
train$Embarked[is.na(train$Embarked)] <- 'S'

#Fareの欠損値を埋める。medianを入れる
#trainには欠損値なし
test$Fare[is.na(test$Fare)] <- median(test$Fare, na.rm=T)

欠損値がなくなったかどうか確認。

#欠損値の数確認（再々）
is_na_train <- sapply(train, function(y) sum(is.na(y)))
is_na_train

>PassengerId    Survived      Pclass        Name         Sex         Age 
>          0           0           0           0           0           0 
>      SibSp       Parch      Ticket        Fare       Cabin    Embarked 
>          0           0           0           0         687           0

Age、Embarkedの欠損値は0になった。テストデータの確認は省略。

次に、文字列型のデータをカテゴリ値に修正する。

#文字列をカテゴリ値に
train$Sex[train$Sex=='male'] <- 0
train$Sex[train$Sex=='female'] <- 1
train$Embarked[train$Embarked=='S'] <- 0
train$Embarked[train$Embarked=='C'] <- 1
train$Embarked[train$Embarked=='Q'] <- 2

test$Sex[test$Sex=='male'] <- 0
test$Sex[test$Sex=='female'] <- 1
test$Embarked[test$Embarked=='S'] <- 0
test$Embarked[test$Embarked=='C'] <- 1
test$Embarked[test$Embarked=='Q'] <- 2

これでデータクリーニングが終わったので、いよいよモデルを作る。

#決定木でモデルを作る
library(rpart)

model1 <- rpart(Survived ~ Pclass + Age + Sex + Fare, data=train)
model1

>n= 891 
>
>node), split, n, deviance, yval
>      * denotes terminal node
>
> 1) root 891 210.727300 0.3838384  
>   2) Sex=0 577  88.409010 0.1889081  
>     4) Age>=6.5 553  77.359860 0.1681736  
>       8) Pclass>=1.5 433  44.226330 0.1154734 *
>       9) Pclass< 1.5 120  27.591670 0.3583333 *
>     5) Age< 6.5 24   5.333333 0.6666667 *
>   3) Sex=1 314  60.105100 0.7420382  
>     6) Pclass>=2.5 144  36.000000 0.5000000  
>      12) Fare>=23.35 27   2.666667 0.1111111 *
>      13) Fare< 23.35 117  28.307690 0.5897436 *
>     7) Pclass< 2.5 170   8.523529 0.9470588 *

うーん、よくわからない。なので、図示してみる。

library(rpart.plot)

#図示してみる
rpart.plot(model1, extra = 1)

出力結果は以下の図。

f:id:prism0081:20180919125511j:plain

これの評価はさておき、作ったモデルをテストデータにあてはめて予測してみる。

#モデルで予測する
pred1 <- predict(model1, test)
pred1

>        1         2         3         4         5         6         7 
>0.1154734 0.5897436 0.1154734 0.1154734 0.5897436 0.1154734 0.5897436 
>        8         9        10        11        12        13        14 
>0.1154734 0.5897436 0.1154734 0.1154734 0.3583333 0.9470588 0.1154734 
（以下略）

0,1での出力じゃなくて、確率で出してくれるのね。kaggeleへの提出は0,1での予測が必要なので変換する。

#0,1で表現する
pred1 <- round(pred1)
pred1

>  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
>  0   1   0   0   1   0   1   0   1   0   0   0   1   0   1   1   0   0   1 
> 20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38 
>  1   0   0   1   0   1   0   1   0   0   0   0   0   1   0   0   0   1   1 
（以下略）

予測結果をkaggleで評価してもらうために、csvファイルに書き出す。

# 予測結果とPassengerIdをデータフレームへ
my_solution1 = data.frame(
  PassengerId = test$PassengerId,
  Survived = pred1
)

#結果のCSV書き出し
write.csv(my_solution1, file='./model1_solution.csv', row.names = F)

出力したcsvファイルをkaggleでアップロードして判定してもらったところ・・・

f:id:prism0081:20180919130530j:plain

まあまあの数字。
まったく同じ事をPythonでやってみたけど、Rの方がスコアは高かった。ライブラリの違いでしょうかね。

少しモデルの説明変数を増やしてみる。

#モデルの説明変数を増やしてみる
model2 <- rpart(Survived ~ Pclass + Age + Sex + Fare + SibSp + Parch + Embarked, data=train)

#図示してみる
rpart.plot(model2, extra = 1)

f:id:prism0081:20180919132255j:plain

大きくは変わってない印象・・・モデル１と同様に予測してみる。

#モデルで予測する
pred2 <- predict(model2, test)

#0,1で表現する
pred2 <- round(pred2)

# 予測結果とPassengerIdをデータフレームへ
my_solution2 = data.frame(
  PassengerId = test$PassengerId,
  Survived = pred2
)

#結果のCSV書き出し
write.csv(my_solution2, file='./model2_solution.csv', row.names = F)

kaggleにアップロードして評価してもらうと

f:id:prism0081:20180919132239j:plain

少しだけ改善しましたが、もっとチューニングが必要です。

2018-09-07

Julia入門：環境構築

Julia

ベイズ学習を勉強するにあたって、何の言語でやろうかなと考えていたところ、Juliaが良いという話を聞いたので勉強してみることにした。

ひとまず環境構築しないといけないのだけれど、あまり流行ってない言語なのか参考情報が少ない＆やや古め・・・
なので、自分用のメモとしてまとめておく。

①Juliaのインストール
公式サイトで自分の環境にあったインストーラーを選んでダウンロード＆実行。
The Julia Language

②実行
インストールが終わると「Julia」というショートカットができるはずなので、それを実行。
すると、以下のようなコンソールが起動する。

f:id:prism0081:20180907171044j:plain

これにごりごりとコードを書いてもいいのでしょうが勉強にはあまりにも不便なので、IDEを導入する。

③Junoのインストール
「Juno」というIDEがあるそうです。以下を参考にさせてもらった。

st-hakky.hatenablog.com

手順としては、まず「Atom」というオープンソースのテキストエディターをインストール。
公式サイトは以下。

atom.io

インストールが完了したら起動。Juliaを扱えるようにするため「uber-juno」というパッケージをインストールする。
Atomのメニューの「File」から「Settings」を選択。
その中に「Install」があるはずなのでそれを選択して、表示される検索ボックスに上記のパッケージ名「uber-juno」を入力する。
すると、当該パッケージがヒットするはずなのでインストールを選択する。

④自分好みに設定する
あとは適宜、自分好みにテーマを変えるとかしてカスタマイズしていけばよい。

こうして箇条書きすると難しくなく見えるが、ちょこちょこ躓いた。
そんなときは、アプリケーションをいったん終了して再起動すると先に進めた。
調べ物の時間含めて、なんだかんだで３時間くらい使ってしまった。

2018-09-06

クラウドソーシング

雑記

最近、研究者の間でクラウドソーシングを使った調査が流行っているという話を聞いた。
気になって調べて見ると、海外で問題があるケースがあったそうで。

Amazonで労働力を購入するMechanical Turkでサーベイを行ったら回答が全くランダムな入力であることが発覚。同じGPS情報でリピートされている。少なくともこれまでの90の研究に影響。

# つまり自動回答で荒稼ぎしている詐欺師がいるということか https://t.co/7rAX8l1NRd
— 高梨陣平 (@jingbay) 2018年8月18日

元記事はこちら。
www.maxhuibai.com

日本でも同様のケースは十分起こりえるので、使用する際には十分注意しないと。

2018-07-17

MRPをめぐるツイート群

雑記

Love seeing more MRP in the wild, but count me skeptical after the fact. Would love to see more pollsters (besides @doug_rivers and @benlauderdale of course), do this in real time! https://t.co/y7QeRFPdzy
— Jonathan Robinson (@jon_m_rob) 2018年7月15日

This is just a function of their low sample size resulting in their biased midwestern estimates being regularized down and their biased southern estimates being regularized up. If they applied this to a 200k dataset then it wouldn’t have worked (as someone who has tried).
— *1 (@davidshor) 2018年7月15日

I RTed this but I have no earthly idea what you mean, I'll just say that up front.
— Xenocrypt (@xenocryptsite) 2018年7月15日

What you mean is that they were right by accident because stuff cancelled out?
— Xenocrypt (@xenocryptsite) 2018年7月15日

I *think* what he means is that the national data was more predictive than the state level data, so using an MR P-d by national data is more accurate post facto. But if they'd had more state data to put in their model, it would have been just as wrong as state poll-heavy models.
— 𝔸lfred 𝕁ones (@rAlfredJones) 2018年7月15日

Because practically they're not really adding anything, just weighting their model more towards the (in 2016) more accurate national polls/census data and away from the (in 2016) less accurate state polls. Which, duh, because the national polls were more accurate.
— 𝔸lfred 𝕁ones (@rAlfredJones) 2018年7月15日

But that regularization only happened because they didn't have much data. If they had 10X the amount of polling data, it wouldn't have done that and the result would have been wrong.
— *2 (@davidshor) 2018年7月15日

You're doing some high level speculation here !
— Jonathan Robinson (@jon_m_rob) 2018年7月15日

but this has always been true. you can go back and find gelman having clinton winning MS in 96 or whatever for this reason
— Nate Cohn (@Nate_Cohn) 2018年7月16日

I'm aware of that as a feature of these kinds of models, I'm just willing to at least speculate about the differences in the DGP's across different survey datasets.
— Jonathan Robinson (@jon_m_rob) 2018年7月16日

Fun thought experiment/real thing that Gallup sort of did at AAPOR '18. Only build an MRP model off a sample/subset of your data so you are more conservative/leaning on the prior to hedge your data having all sorts of biases
— Jonathan Robinson (@jon_m_rob) 2018年7月16日

I think in most cycles the optimal amount of data to throw away is going to be zero and that 2016 was a weird fluke.
— *3 (@davidshor) 2018年7月16日

Well Gallup was analyzing their 2012 data so...
— Jonathan Robinson (@jon_m_rob) 2018年7月16日

2012 seems like a perfect example of how focusing on correcting non-response bias and collecting as much data as you can is going to do better than messing around with MRP pic.twitter.com/FrLppsZFn1
— *4 (@davidshor) 2018年7月16日

I mean you can do both
— Jonathan Robinson (@jon_m_rob) 2018年7月16日

that MRP isn't a magic fix for poor sampling seems rather obvious to me?
— (((Will Cubbison))) וויל קביסון (@wccubbison) 2018年7月16日

There's also just that MRP itself is like, a 2009-era methodology. The core insight of "data-pooling and regularization is good" was revolutionary at the time, but we have better tools and a very different data situation now that do that better.
— *5 (@davidshor) 2018年7月16日

MRP strikes me as an already antiquated methodology applicable primarily to pollsters who still refuse to use voter files
— Patrick Ruffini (@PatrickRuffini) 2018年7月16日

What is stopping you from using voter files on MRP?
— doug rivers (@doug_rivers) 2018年7月16日

Or put differently, having a voter file and using MRP are solutions to different parts of the problem.
— Benjamin Lauderdale (@benlauderdale) 2018年7月16日

Having a voter file enables me to generate model scores for every voter, and then take geographic aggregates from individual scores
— Patrick Ruffini (@PatrickRuffini) 2018年7月16日

i think half of the people referring to 'mrp' are doing that, even if it's not formally post-stratification
— Nate Cohn (@Nate_Cohn) 2018年7月16日

I think of MRP as a convoluted workaround where you don't have a census of all voters you can apply estimates to, if you're doing RDD for example
— Patrick Ruffini (@PatrickRuffini) 2018年7月16日

sure, though you can do the same thing on ACS or CPS microdata without formally stratifying it
— Nate Cohn (@Nate_Cohn) 2018年7月16日

bringing back in a paper author (@chadkdj) in case they want to respond to the comments in this threat; in particular, the method applied to other past elections
— Matt Grossmann (@MattGrossmann) 2018年7月16日

It *is* post-stratification with cells of size one. The whole point is that the sample cells can be small or empty, but you can still use a model to make cell level estimates. https://t.co/zs7sflPYuT
— doug rivers (@doug_rivers) 2018年7月16日

Agreed. In particular, one should be wary of sparse solutions, which don't make any sense for some representations. https://t.co/eJr5EcKHg0
— doug rivers (@doug_rivers) 2018年7月16日

*1:David Shor

*2:David Shor

*3:David Shor

*4:David Shor

*5:David Shor

2017-12-13

勉強用の本

書籍・資料

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル原著第2版

ベイズ統計で実践モデリング: 認知モデルのトレーニング

プラレールトミカと遊ぼう！くるぞわたるぞ！カンカン踏切セット

Apple MD826AM/A Lightning - Digital AVアダプタ

【ふるさと納税】高評価★4.59 発送時期が選べるエンペラーサーモン【1kg】サーモン鮭刺身さけサケスモークサーモン魚人気ふるさと北海道海鮮送料無料ランキング多数入賞

【ふるさと納税】【累計2,000万個突破!!】鉄板焼ハンバーグデミソース 20個温めるだけ福岡飯塚牛冷凍小分け大容量ハンバーグ飯塚市はんばーぐ肉簡単調理デミグラスソース特製湯煎人気子供【A2-074】【1110レビューCP】（クラウドファンディング対象）

【ふるさと納税】北海道十勝スイーツ「三方六の小割」10本入り4箱本別町観光協会送料無料《60日以内に順次出荷(土日祝除く)》

【ふるさと納税】ROYCE'人気スイーツ詰め合わせ（A)

【ふるさと納税】【発送時期が選べる】中山牧場佐賀牛赤身スライス 1kg

【ふるさと納税】高レビュー4.7以上!!＜選べるお届け月＞国産うなぎ蒲焼4尾(計800g以上) 数量限定鰻蒲焼ウナギ蒲焼用たれさんしょうのセット(うなぎ1尾190g以上の鰻4尾からなるウナギの詰め合わせ)

【ふるさと納税】【2023年出荷予約】北海道産赤肉メロン大玉 2玉

【ふるさと納税】羽毛布団掛け布団ゴールドラベルシングル日本製合い掛け洗えるコロナ支援訳あり選べる配送時期ホワイトダックダウン90％洗濯可寝具布団国産甲州羽毛ふとんダンボールハウスつながるふじよしだ（クラウドファンディング対象）

【ふるさと納税】ついに復活！※数量限定※九州産黒毛和牛切り落とし 2000g（500g×4パック）小分け合計2kg 期間限定切落し牛肉黒毛和牛冷凍国産九州産送料無料応援緊急支援おすすめランキング

【ふるさと納税】【背ワタなし】大型むきえび（高級ブラックタイガー）約1kg/約50〜80尾【甲羅組蟹カニ海老エビ】【発送時期をお選びください】

【ふるさと納税】シュークリーム【最大6か月待ち】かみのやまシュー 6個お菓子おやつデザートスイーツスウィーツ洋菓子小分け一人暮らし少人数お取り寄せグルメ 3000円冷凍配送山形県 0048-2209

【ふるさと納税】大地の黄金干し芋計200g(100g×2袋)鹿児島県産のさつまいも「紅はるか」を使用した、しっとり、こく甘、やわらかな黄金色のほしいも!無添加・無着色な健康的なおいものスイーツ! 【末永商店】p3-003

【ふるさと納税】【ルタオ】ドゥーブルフロマージュお菓子チーズケーキルタオドゥーブルフロマージュスイーツ TV メディア北海道ふるさと納税千歳市ふるさと納税【北海道千歳市】お中元北海道セットギフトプレゼント

【ふるさと納税】柳月「三方六」詰合せ秋冬限定セット（珈琲＆ミルク）本別町観光協会《60日以内に順次出荷(土日祝除く)》三方六菓子バームクーヘン洋菓子コーヒーミルク送料無料

早稲田大学（社会科学部）（2023年版大学入試シリーズ） [ 教学社編集部 ]

上智大学（神学部・文学部・総合人間科学部）（2023年版大学入試シリーズ） [ 教学社編集部 ]

ポイント最大26倍！！ Panasonic ななめドラム式洗濯乾燥機乾燥フィルター用おそうじブラシ AXW22R-9DA0 掃除ブラシパナソニック AXW22R9DA0 純正品送料無料【SK01559】

ふと浮かぶ記憶と思考の心理学: 無意図的な心的活動の基礎と臨床

NS67 バスティンコンピレーションアルバムレベル2

赤青鉛筆

ビット・バイ・ビット　デジタル社会調査入門

スーパーマリオブラザーズワンダー

iPad Air5用ペーパーライクフィルム

傷モノの花嫁 ~虐げられた私が、皇國の鬼神に見初められた理由~(1)

DABADA(ダバダ) エスボード光るタイヤ組み立て済みプロテクター 3点セット付きブレイブボードブレイドボードキャスターボードブレボースケボースケートボード ESSボード ESSBOARD 子供キッズ

世界の果ては深愛(2) (KCデラックス)

ベイブレードX　バトルエントリーセット

オペレーションズリサーチ

書道用下敷き

Z会中学受験シリーズ入試に出る動物図鑑改訂版

Z会中学受験シリーズ入試に出る植物図鑑改訂版

ウケるゴロ合わせ《日本史編》　イヤでも覚える基本重要事項98

WP70J ファーストピアノレパートリーアルバム

山善(YAMAZEN) 突っ張りポールハンガーラックブラック幅97.5-166cm 高さ140-280cm ダブル伸縮 WJ-775(BK)

【BLKP】パール金属ハンガーラック限定マットブラックシングル耐荷重 20kg 幅 80cm 極太パイプキャスター付き高さ調節可能 BLKP 黒洗濯室内干し N-7573

Q.E.D.iff -証明終了-(27) (講談社コミックス月刊マガジン)

ないない堂 ~タヌキ和尚の禍事帖~(3) (講談社コミックス月刊マガジン)

ないない堂 ~タヌキ和尚の禍事帖~(2) (講談社コミックス月刊マガジン)

BBC micro:bit V2.21（マイクロビット) バルク品正規品

iPhone SE 第3世代第2世代ケースショルダーシリコン iphone 8/7 スマホケースアイフォン7/8/se2/se3 カバー耐衝撃肩掛け斜めかけ縄掛けケース紐付き斜めがけ織りストラップ付き長さ調整可能取り外し可能傷付き防止マット質感指紋防止(グレーブルー)

こどものチェルニー (標準版)

ソードアート・オンライン28 ユナイタル・リングVII (電撃文庫)

ミステリと言う勿れ (14) (フラワーコミックスα)

中学受験国語文章読解の鉄則増補改訂版

キヤノン BCI-301+300/5MP 小

Excelで今すぐはじめる心理統計簡単ツールHADで基本を身につける (KS心理学専門書)

ふくしま式「国語の読解問題」に強くなる問題集〔小学生版〕

コクヨ鉛筆シャープ芯径1.3mm 黒 PS-P101D-1P

コクヨ(KOKUYO) シャープペンシャープ替芯 1.3mm 2B 3個セット PSR-2B13-1PX3SET

コクヨ鉛筆シャープ芯径1.3mm ダークグリーン PS-P101DG-1P

市川中学校　2025年度用 5年間（＋３年間ＨＰ掲載）スーパー過去問（声教の中学過去問シリーズ 352）

東邦大学付属東邦中学校　2025年度用 4年間（＋３年間ＨＰ掲載）スーパー過去問（声教の中学過去問シリーズ 351）

早稲田大学高等学院中学部　2025年度用 9年間（＋３年間ＨＰ掲載）スーパー過去問（声教の中学過去問シリーズ 81）

慶應義塾中等部　2025年度用 10年間（＋３年間ＨＰ掲載）スーパー過去問（声教の中学過去問シリーズ 5）

高輪中学校　2025年度用 3年間（＋３年間ＨＰ掲載）スーパー過去問（声教の中学過去問シリーズ 105）

東京学芸大学附属竹早中学校　2025年度用 7年間（＋３年間ＨＰ掲載）スーパー過去問（声教の中学過去問シリーズ 13）

金の国水の国 (フラワーコミックスアルファスペシャル)

オーエ組み合わせ風呂ふたアイボリー 2枚組:幅73×長さ158cm用防カビ抗菌日本製 L-16 2枚入

2025年中学入試用サピックス重大ニュース

昭和学院秀英中学校　2025年度用 4年間（＋３年間ＨＰ掲載）スーパー過去問（声教の中学過去問シリーズ 362）

小学 100%丸暗記白地図&用語:マンガと地図で楽しく覚える! (受験研究社)

スーパーマリオパーティジャンボリー - Switch

シャーペン

ぺんてるシャープ芯ぺんてるアイン 0.9mm HB 5個 XC289-HB

[サクラクレパス] ミーティングバッグノータムオフィストートバッグ A4 ブラック UNT-A4#49

大きなシリコンマフィン型ケーキ型プレートノンスティック 12個取りマフィンパンベーキングカップケーキモールドキッチンツール掃除が簡単グレー

コクヨペンケース筆箱大容量ツールペンケースカバコ

[ミズノ] ランドリーバッグスポーツ部活小物 33JMB002 ブラック/ゴールド

コカ・コーラやかんの麦茶 from 爽健美茶ラベルレス 650mlPET×24本

コカ・コーラ綾鷹 525mlPET×24本

[グンゼ] 半袖シャツ男の子半袖シャツセット組子供肌着オールシーズン年間綿100 コットン部屋干し抗菌防臭キッズ

BOTANIST ボタニスト | シャンプースカルプクレンズ 460ml

心理学検定一問一答問題集[A領域編]

コカ・コーラやかんの麦茶 from 爽健美茶ラベルレス 650mlPET×24本

任意のテキスト

まんがで身につくめざせ! あしたの算数王 (1) 可能性の数、確率

カラフルピーチと一緒に攻略!わくわく中学英単語記憶に残る!英単語&熟語1800

素数ゼミの謎

【Amazon.co.jp限定】DEAR MYSTERIES(CD+Blu-ray) - TOMOO (メガジャケ日本武道館LIVEver. 19cm×19cm＋早期予約メーカー特典『DEAR MYSTERIES』オリジナルZINE）

フリーリーキッチンマットクリア透明マット拭ける抗菌防臭厚1.5mm 床保護マット傷防止マット撥水クリア水拭き PVCマットビニール台所 (80×270cm)

2017-09-06

Chapter5 練習問題 Part2

StanとRでベイズ統計モデリング

(3)5.3.2項の続きで、アルバイトが好きかどうかごとにYを集計せよ

元データはこんなもの。

data-attendance-3.txt

PersonID,A,Score,Weather,Y
1,0,69,B,1
1,0,69,A,1
1,0,69,C,1
1,0,69,A,1
1,0,69,B,1
1,0,69,B,1
1,0,69,C,0
1,0,69,B,1
1,0,69,A,1

コードは以下のとおり。

#データ読み込み
d2 <- read.csv(file="data-attendance-3.txt")
#aggregate(x=list(<表示列名>=<変数名>$<集計列>),by=list(<表示列名>=<変数名>$<キー列>,<表示列名>=<変数名>$<キー列>,……),FUN=sum)

aggregate(x=list(Y=d2$Y),by=list(A=d2$A),FUN=table) #自分の回答

aggregate(Y~A,data=d2,FUN=table) #回答はこっち

結果はどちらも一致。

  A Y.0 Y.1
1 0 288 994
2 1 386 728

(4)5.3.3項では曇りと雨の影響を固定した。晴れを０とし、曇りと雨の影響をパラメーター化したモデルも考えられる。Stanで推定せよ。

まずモデル式を本から書き写すと

$q[i]=inv\_logit(b_{1}+b_{2}A[i]+b_{3}Score[i]+b_{4}Weather[i])~~~i=1,\cdots,I$
$Y[i]=Bernoulli(q[i])~~~i=1,\cdots,I$

これに基づくStanプログラムは

model5-5.stan

data{
  int I;
  int<lower=0, upper=1> A[I];
  real<lower=0, upper=1> Score[I];
  real<lower=0, upper=1> W[I];
  int<lower=0, upper=1> Y[I];
}

parameters{
  real b[4];

}

transformed parameters{
  real q[I];
  for(i in 1:I)
    q[i]=inv_logit(b[1]+b[2]*A[i]+b[3]*Score[i]+b[4]*W[i]);
}

model{
  for(i in 1:I)
    Y[i]~bernoulli(q[i]);
}

これを修正して天気をパラメーターに入れる。

model5-5-2.stan

data{
  int I;
  int<lower=0, upper=1> A[I];
  real<lower=0, upper=1> Score[I];
  int<lower=1, upper=3> WID[I];
  int<lower=0, upper=1> Y[I];
}

parameters{
  real b[3];
  real bw2;
  real bw3;
}

transformed parameters{
  real q[I];
  real bw[3];
  bw[1]=0; //晴れは0
  bw[2]=bw2;
  bw[3]=bw3;
  for(i in 1:I)
    q[i]=inv_logit(b[1]+b[2]*A[i]+b[3]*Score[i]+bw[WID[i]]);
}

model{
  for(i in 1:I)
    Y[i]~bernoulli(q[i]);
}

実行コードは以下のとおり

conv <- c(1,2,3)
names(conv)<- c('A','B','C')
data <- list(I=nrow(d2),A=d2$A, Score=d2$Score/200, WID=conv[d2$Weather], Y=d2$Y)
fit4 <- stan(file="model5-5-2.stan", data=data, seed=1234)

トレースプロットは・・・

f:id:prism0081:20170906183700j:plain

Stanコード間違って何度かやり直したけど、ようやくちゃんと収束しました。

続きはまた後日。

勉強などのメモ

勉強用のメモ

おすすめ漫画