順序制約のあるrstan

推定するパラメータに順序制約があるとき、rstan ではordered、正に限定するならpositive_ordered が使える。
　
例えばYaegashi (J Hum Genet. 1998;43(2):85-90.)らは、出生時の染色体異常が母体年齢に応じてどう推移するかをデータを取って調べている。
ここで、筆者らは というモデルを考えている。

txt <- "
   cases +21
35   736   1
36   850   4
37   824   8
38   800   2
39   713   8
40   626   6
41   417   3
42   255   6
43   146   2
44    87   2
45    22   0
46     7   0
47     1   0
"
yaegashi <- read.table(text=txt, check.names=FALSE)
y <- head(yaegashi$"+21"/yaegashi$cases, -3)
x <- 35:44
lm0 <- lm(log(y) ~ x)

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, cex.axis=1.5, las=1)
plot(x, y, log="y", ylim=c(0.5, 100)/1000, yaxt="n", pch=16, xlab="Maternal age", ylab="Rate per 1000", cex=2)
axis(2, at=c(1,5,10,50,100)/1000, labels=c(1,5,10,50,100))
lines(x, exp(cbind(1, x) %*% lm0$coefficients), lwd=3)
# exp(lm0$fitted.value)

lm0

Call:
lm(formula = log(y) ~ x)

Coefficients:
(Intercept)            x  
   -14.5127       0.2447  

なので、係数 は論文のTrisomy 21 の値 0.245 と一致している。
図も回帰線とプロットをそれなりに真似るとこんな感じ。

　
ここで、45週以降のTrisomy 21 の発生は0件なので、 の回帰では45週以降の発生率の観測値が0 になるので、線形回帰ではよろしくない。というわけで45週以降は欠測値として回帰したが、rstan による状態空間モデルにより、45週以降の発生率も推定してみる。
　
芸がなく2階差分でモデルを作る。
母体年齢 に応じた発生率は で表されるとして、ある年齢の妊婦 人から実際に観測される症例数は二項分布に従う、とすれば、以下のように書ける。
（ここではまだ順序制約はいれていない）

code <- "
data{
  int Nage;
  int age[Nage];
  int Pos[Nage];
  int Ntest[Nage];
}
parameters{
  real<lower=0, upper=1> p[Nage];
  real<lower=0> s;
}
model{
  for(i in 3:Nage){
    p[i] ~ normal(2*p[i-1] + p[i-2], s);
  }
  for(i in 1:Nage){
    Pos[i] ~ binomial(Ntest[i], p[i]);
  }
}

"

# sampling
standata <- list(Nage=nrow(yaegashi), age=as.numeric(rownames(yaegashi)), Pos=yaegashi$"+21", Ntest=yaegashi$cases)
fit_yaegashi <- sampling(m, standata, warmup=1000, iter=5000, seed=1234)

　
推定結果は下に図示するが、年齢に応じて凸凹している。というのも、38歳で突如、観測された発生率が落ち込むので、これに引きずられて推定値も下がる。
というわけで、年齢に応じて発生率が凸凹するのはなんか気持ち悪い、と思ったときに、年齢に応じた順序制約を入れる。
これはrstan ではordered もしくはpositive_ordered で作ることができる。発生率 は0以上1以下の確率なので、positive_ordered が良さそう、と思って

positive_ordered<lower=0, upper=1>

みたいなことをするのだが、rstan はこの書き方を受け付けないらしい。
とういわけでordered からtransformed parameters で変換するのが手らしい。しかしコンパイル時に警告をはく。無視。
　
実数値を[0,1] に変換したければ、ロジットの出番である。

# positive ordered に書き換えたコード
code <- "
data{
  int Nage;
  int age[Nage];
  int Pos[Nage];
  int Ntest[Nage];
}
parameters{
  //real<lower=0, upper=1> p[Nage];
  ordered[Nage] pp;
  real<lower=0> s;
}
transformed parameters{
  real<lower=0, upper=1> p[Nage];  // 変換されたパラメータ
  for(i in 1:Nage){
    p[i] = 1/(1+exp(-pp[i]));  // ロジットにすれば[0,1] に強制的に収まる
  }
}
model{
    for(i in 3:Nage){
      p[i] ~ normal(2*p[i-1] + p[i-2], s);
  }
  for(i in 1:Nage){
      Pos[i] ~ binomial(Ntest[i], p[i]);
    }
}

　
順序制約のないときはジグザグ凸凹な推定だったが、順序制約をいれると年齢に応じて発生率が上昇するわりときれいめな推定結果になった。
また、実際の観測数が0 である45週以上の推定も何とかできる。ただし、ここはということを利用しているだけなので、数が少ないため信用区間はかなり広い。

援助交際してそうなアニメキャラランキング2017をstanで考える

注意：本解析の結果と、実際に声優がそうであるかはまったく関係がありません。
　
援助交際してそうなアニメキャラランキング2016をstanで考える - 驚異のアニヲタ社会復帰への道
援助交際してそうなアニメキャラランキング2015をstanで考える - 驚異のアニヲタ社会復帰への道
これらの続き。
　
援助交際してそうなアニメキャラランキング2017 が発表されたので、投票数とアニメと円盤売上と声優から援助交際のしやすさをstanで推定する。
モデルは2016のページをそのまま丸コピしつつ、データも2013-2016 年度分を追加する。
　
キャラについて
新子憧が前人未到??の6連覇を達成したが、新子憧の得票率としてはもう均衡に達していて、7%程度の得票率にとどまるようだ。

　
声優について
新子憧のせいで東山奈央がごにょごにょ声優なのではというひどい風評被害があるが、2017年に限っていえば東山奈央より係数が高かったのは201人中45人いた。
ラブライブサンシャインの声優や、NEW GAME声優が挙がった。

  [1] "本渡楓"                 "福圓美里"               "小山茉美"              
  [4] "長縄まりあ"             "高田憂希"               "山口愛"                
  [7] "朝日奈丸佳"             "田中美海"               "與那嶺里都"            
 [10] "杉本ゆう"               "長久友紀"               "加隈亜衣"              
 [13] "斉藤朱夏"               "逢田梨香子"             "和久井優"              
 [16] "降幡愛"                 "小林愛香"               "日岡なつみ"            
 [19] "諏訪ななか"             "伊波杏樹"               "高槻かなこ"            
 [22] "田口宏子"               "和氣あず未"             "楠田亜衣奈"            
 [25] "南條愛乃"               "高部あい"               "名塚佳織"              
 [28] "能登麻美子"             "赤崎千夏"               "村川梨衣"              
 [31] "千菅春香"               "福原綾香"               "鳴海杏子"              
 [34] "木戸衣吹"               "千本木彩花"             "上坂すみれ"            
 [37] "新田恵海"               "内田真礼"               "藤村鼓乃美"            
 [40] "徳井青空"               "内田彩"                 "渡部優衣"              
 [43] "高森奈津美"             "茅原実里"               "野中藍"    

　
東山奈央については、係数はほぼ常に援交側に正である。他には、2012年から2017年まですべての年度でリストアップされる声優については喜多村英梨、竹達彩奈などがいるが、彼女らについても、2012年はデータそのものが少ないため推定範囲が広く、2017年に近づくと推定がよくなる、という傾向だった。

　
アニメについて
219作品がリストアップされたうち、有意に援交しないアニメと援交するアニメをプロットした。
ここに挙がっているアニメはほとんどが1-2回登場しただけなので、例えば咲とかラブライブとか俺ガイルとか艦これとか、そういう登場頻度が多いやつは意外と有意にならなかった。
なお、2016年度の解析ではハイスクール・フリートとマクロスF は援交しにくそうなアニメだったようだが、今回は援交アニメになってしまった。

　
円盤売上について
2015年では売上が多いアニメは援交しなさそうだったが、今回は特にどちらでもない結果だった。
　
i7, 4core 並列、1000 burnin, 1000 iteration で30分くらい。 は全然収束しませんでした。
2017年度分を追加

vote,name,anime,cv,bd,year
1021,新子憧,咲-Saki-,東山奈央,4939,2017
501,渋谷凛,アイドルマスターシンデレラガールズ,福原綾香,14952,2017
462,烏丸千歳,ガーリッシュナンバー,千本木彩花,659,2017
356,高坂穂乃果,ラブライブ!,新田恵海,64911,2017
311,スノーホワイト/姫河小雪,魔法少女育成計画,東山奈央,1329,2017
306,一色いろは,やはり俺の青春ラブコメはまちがっている。,佐倉綾音,12968,2017
268,神野めぐみ,エロマンガ先生,木戸衣吹,8723,2017
265,由比ヶ浜結衣,やはり俺の青春ラブコメはまちがっている。,東山奈央,12968,2017
264,傘木希美,響け!ユーフォニアム,東山奈央,7854,2017
227,青羽ここな,ヤマノススメ,小倉唯,4249,2017
212,桐間紗路,ご注文はうさぎですか?,内田真礼,11308,2017
209,九条カレン,きんいろモザイク,東山奈央,6582,2017
209,ゆんゆん,この素晴らしい世界に祝福を!,豊崎愛生,10120,2017
184,西川葉子,三者三葉,和久井優,2617,2017
181,中川夏紀,響け!ユーフォニアム,藤村鼓乃美,7854,2017
177,丹生谷森夏,中二病でも恋がしたい!,赤崎千夏,15466,2017
172,アクア,この素晴らしい世界に祝福を!,雨宮天,10120,2017
171,沖田紗羽,TARITARI,早見沙織,8389,2017
171,伊万莉まりあ,この美術部には問題がある!,東山奈央,1143,2017
170,森川千夏,Lostorage incited WIXOSS,井口裕香,3000,2017
170,川神舞,無彩限のファントム・ワールド,上坂すみれ,2336,2017
169,ジータ,グランブルーファンタジージ・アニメーション,金元寿子,52553,2017
167,めぐみん,この素晴らしい世界に祝福を!,高橋李衣,10120,2017
167,渡辺曜,ラブライブ!サンシャイン!!,斉藤朱夏,54312,2017
166,涼風青葉,NEWGAME!,高田憂希,5864,2017
164,雪ノ下雪乃,やはり俺の青春ラブコメはまちがっている。,早見沙織,12968,2017
163,天真=ガヴリール=ホワイト,ガヴリールドロップアウト,富田美憂,3769,2017
162,矢澤にこ,ラブライブ!,徳井青空,64911,2017
160,南ことり,ラブライブ!,内田彩,64911,2017
159,桜内梨子,ラブライブ!サンシャイン!!,逢田梨香子,54312,2017
153,園田海未,ラブライブ!,三森すずこ,64911,2017
134,磯野フネ,サザエさん,麻生美代子,3000,2017
117,ギャル子,おしえて!ギャル子ちゃん,和氣あず未,1394,2017
117,神代小蒔,咲-Saki-,早見沙織,4565,2017
110,前川みく,アイドルマスターシンデレラガールズ,高森奈津美,14952,2017
109,鈴乃木凜,ばくおん!!,東山奈央,2826,2017
109,源内あお,フレームアームズ・ガール,日笠陽子,8398,2017
105,常木耀,セイレン,佐倉綾音,933,2017
105,一条蛍,のんのんびより,村川梨衣,8375,2017
105,中世古香織,響け!ユーフォニアム,茅原実里,7854,2017
103,坂木しずか,SHIROBAKO,千菅春香,15210,2017
96,倉田亜美,ろんぐらいだぁす！,東山奈央,1131,2017
94,呉織あぎり,キルミーベイベー,高部あい,993,2017
91,田中あすか,響け!ユーフォニアム,寿美菜子,7854,2017
91,宮永咲,咲-Saki-,植田佳奈,7758,2017
90,城ヶ崎美嘉,アイドルマスターシンデレラガールズ,佳村はるか,14952,2017
87,本田珠輝,ステラのまほう,長縄まりあ,631,2017
86,蒼井晶,selector infected WIXOSS,赤崎千夏,3873,2017
86,四ノ宮しおり,サクラクエスト,上田麗奈,1392,2017
86,高坂雪穂,ラブライブ!,東山奈央,64911,2017
80,ルリア,グランブルーファンタジージ・アニメーション,東山奈央,52553,2017
80,千反田える,氷菓,佐藤聡美,9795,2017
79,サーバル,けものフレンズ,野中藍,12650,2017
74,保登心愛,ご注文はうさぎですか?,佐倉綾音,11038,2017
74,松実玄,咲-Saki-,花澤香菜,4939,2017
73,ざんげちゃん,かんなぎ,花澤香菜,10880,2017
67,黒澤ルビィ,ラブライブ!サンシャイン!!,降幡愛,54312,2017
62,香風智乃,ご注文はうさぎですか?,水瀬いのり,11038,2017
62,津島善子,ラブライブ!サンシャイン!!,小林愛香,54312,2017
61,胡桃沢＝サタニキア＝マクドウェル,ガヴリールドロップアウト,大空直美,3769,2017
60,滝本ひふみ,NEWGAME!,山口愛,5864,2017
53,ハードゴア・アリス/鳩田亜子,魔法少女育成計画,日高里菜,1329,2017
50,リリルカ＝アーデ,ダンジョンに出会いを求めるのは間違っているだろうか,内田真礼,4941,2017
49,坂凪綾名,魔法少女育成計画,水瀬いのり,1329,2017
47,小鳥遊美羽,パパのいうことを聞きなさい!,喜多村英梨,3470,2017
45,白玉みかん,プリパラ,渡部優衣,728,2017
44,ガルーラ,ポケットモンスター,その他,2514,2017
43,西木野真姫,ラブライブ!,Pile,64911,2017
43,金剛,艦隊これくしょん-艦これ-,東山奈央,18031,2017
43,加藤恵,冴えない彼女の育てかた,安野希世乃,9602,2017
41,相楽エミ,ガールフレンド仮,東山奈央,2282,2017
40,清水谷竜華,咲-Saki-,石原夏織,4939,2017
39,ラ・ピュセル,魔法少女育成計画,佐倉綾音,1329,2017
36,皆川茜,クズの本懐,豊崎愛生,1329,2017
36,高山春香,桜Trick,戸松遥,2110,2017
35,佐久間まゆ,アイドルマスターシンデレラガールズ,牧野由依,14952,2017
35,西条ほのか,おくさまが生徒会長!+！,徳井青空,1055,2017
35,小泉花陽,ラブライブ!,飯田里穂,64911,2017
34,ナルメア,グランブルーファンタジージ・アニメーション,M・A・O,52553,2017
34,周防天音,グリザイアの果実,田口宏子,3422,2017
34,原村和,咲-Saki-,小清水亜美,7758,2017
31,綾瀬花日,12歳。,加隈亜衣,500,2017
31,山田エルフ,エロマンガ先生,高橋未奈美,8723,2017
31,結城明日奈,ソードアート・オンライン,戸松遥,35879,2017
30,藤宮香織,一週間フレンズ。,雨宮天,3718,2017
28,ミルタンク,ポケットモンスター,その他,2514,2017
28,松浦果南,ラブライブ!サンシャイン!!,諏訪ななか,54312,2017
28,高坂桐乃,俺の妹がこんなに可愛いわけがない,竹達彩奈,23888,2017
25,ベルダンディー,ああっ女神さまっ,能登麻美子,1449,2017
24,ダスティネス・フォード・ララティーナ,この素晴らしい世界に祝福を!,茅野愛衣,10120,2017
24,レンホウ,銀魂,その他,15488,2017
23,フーカ・レヴェントン,ViVid Strike!,水瀬いのり,4262,2017
22,梅沢愛子,ダンガンロンパ3TheEndof希望ヶ峰学園絶望編,本渡楓,3458,2017
22,ひまわり,リルリルフェアリル,内田彩,500,2017
22,大場カオル,逆転裁判〜その「真実」、異議あり！〜,杉本ゆう,500,2017
21,高槻やよい,アイドルマスター,仁後真耶子,8684,2017
21,松実宥,咲-Saki-,MAKO,4939,2017
20,音羽歌苗,クラシカロイド,小松未可子,235,2017
20,萩野千秋,ひなこのーと,東城日沙子,1327,2017
20,東條希,ラブライブ！,楠田亜衣奈,64911,2017
20,新垣あやせ,俺の妹がこんなに可愛いわけがない,早見沙織,23888,2017
20,成瀬優,私がモテないのはどう考えてもお前らが悪い!,花澤香菜,1028,2017
19,城ヶ崎莉嘉,アイドルマスターシンデレラガールズ,山本希望,14952,2017
19,新田美波,アイドルマスターシンデレラガールズ,洲崎綾,14952,2017
19,桜木ひな子,ひなこのーと,M・A・O,1327,2017
19,高鴨穏乃,咲-Saki-,悠木碧,4939,2017
18,高宮なすの,てーきゅう,鳴海杏子,3000,2017
17,吉川優子,響け!ユーフォニアム,山岡ゆり,7854,2017
17,大星淡,咲-Saki-,斎藤千和,4939,2017
16,荻布美夏,クロムクロ,瀬戸麻沙美,3000,2017
16,リーゼロッテ＝シャルロック,トリニティセブン,東山奈央,3151,2017
16,千斗いすず,甘城ブリリアントパーク,加隈亜衣,6225,2017
16,大連寺鈴鹿,東京レイヴンズ,佐倉綾音,2033,2017
15,鈴木エリカ,サクラクエスト,黒沢ともよ,1392,2017
15,萌咲いちご,それが声優!,長久友紀,1074,2017
15,ヘスティア,ダンジョンに出会いを求めるのは間違っているだろうか,水瀬いのり,4941,2017
15,中島ゆあ,ひなこのーと,高野麻里佳,1327,2017
15,来栖加奈子,俺の妹がこんなに可愛いわけがない,田村ゆかり,23888,2017
15,押水菜子,花咲くいろは,豊崎愛生,8530,2017
15,黄前久美子,響け!ユーフォニアム,黒沢ともよ,7854,2017
14,リンネ・ベルリネッタ,ViVid Strike!,小倉唯,4262,2017
14,武部沙織,ガールズ&パンツァー,茅野愛衣,35235,2017
14,喜屋武七緒,はいたい七葉,與那嶺里都,419,2017
14,夏川くいな,ひなこのーと,富田美憂,1327,2017
14,絢瀬絵里,ラブライブ！,南條愛乃,64911,2017
14,大和,艦隊これくしょん-艦これ-,竹達彩奈,18031,2017
14,中川かのん,神のみぞ知るセカイ,東山奈央,3040,2017
13,雨宿まち,くまみこ,日岡なつみ,2079,2017
13,黒須あろま,プリパラ,牧野由依,728,2017
13,ネリー・ヴィルサラーゼ,咲-Saki-,その他,7758,2017
12,シャーロット・リンリン,ONE PIECE,小山茉美,3938,2017
12,星井美希,アイドルマスター,長谷川明子,8684,2017
12,島村卯月,アイドルマスターシンデレラガールズ,大橋彩香,14952,2017
12,青山七海,さくら荘のペットな彼女,中津真莉子,2834,2017
12,アママイコ,ポケットモンスター,その他,2514,2017
12,榛名,艦隊これくしょん-艦これ-,東山奈央,18031,2017
12,小笠原晴香,響け!ユーフォニアム,早見沙織,7854,2017
12,園田涼子,月がきれい,東山奈央,3000,2017
12,園城寺怜,咲-Saki-,小倉唯,7758,2017
11,ぽのか,あいまいみー,茅野愛衣,715,2017
11,アリス・カータレット,きんいろモザイク,田中真奈美,6582,2017
11,翔鶴,艦隊これくしょん-艦これ-,野水伊織,18031,2017
11,柏崎星奈,僕は友達が少ない,伊藤かな恵,11111,2017
10,西住みほ,ガールズ&パンツァー,渕上舞,35235,2017
10,磯野ワカメ,サザエさん,津村まこと,3000,2017
10,町子リョウ,幸腹グラフィティ,佐藤利奈,1279,2017
10,石戸霞,咲-Saki-,大原さやか,4565,2017
9,赤沢泉美,Another,米澤円,2157,2017
9,千寿ムラマサ,エロマンガ先生,大西沙織,8723,2017
9,高砂智恵,エロマンガ先生,石川由依,8723,2017
9,弱井トト子,おそ松さん,遠藤綾,67001,2017
9,平沢唯,けいおん!,豊崎愛生,43883,2017
9,玉置亜子,ネトゲの嫁は女の子じゃないと思った?,日高里菜,2251,2017
9,佐々木千穂,はたらく魔王さま！,東山奈央,11511,2017
9,国木田花丸,ラブライブ!サンシャイン!!,高槻かなこ,54312,2017
9,高海千歌,ラブライブ!サンシャイン!!,伊波杏樹,54312,2017
9,那珂,艦隊これくしょん-艦これ-,佐倉綾音,18031,2017
9,鈴谷,艦隊これくしょん-艦これ-,ブリドカットセーラ恵美,18031,2017
9,宮永照,咲-Saki-,中原麻衣,7758,2017
9,上矢あがり,灼熱の卓球娘,田中美海,902,2017
9,小森しゅり,小森さんは断れない!,内田彩,500,2017
8,天海春香,アイドルマスター,井口裕香,8684,2017
8,桐崎千棘,ニセコイ,東山奈央,5794,2017
8,勝田聡子,ハイスクール・フリート,五月ちさと,7909,2017
8,星空凛,ラブライブ!,飯田里穂,64911,2017
8,イリーナ・イェラビッチ,暗殺教室,伊藤静,2184,2017
8,愛宕,艦隊これくしょん-艦これ-,東山奈央,18031,2017
8,無名,甲鉄城のカバネリ,千本木彩花,7512,2017
8,ミナエル,魔法少女育成計画,松田利冴,1329,2017
8,ユナエル,魔法少女育成計画,松田颯水,1329,2017
8,ベルモット,名探偵コナン,小山茉美,500,2017
7,バラライカ,BLACK LAGOON,小山茉美,5491,2017
7,森園立夏,D.C.III〜ダ・カーポIII〜,新田恵海,1937,2017
7,桜ねね,NEWGAME!,朝日奈丸佳,5864,2017
7,安城鳴子,あの日見た花の名前を僕達はまだ知らない。,戸松遥,31019,2017
7,田井中律,けいおん!,佐藤聡美,43883,2017
7,瀧エリ,けいおん!,山村響,43883,2017
7,アライグマ,けものフレンズ,大空直美,12650,2017
7,コツメカワウソ,けものフレンズ,大和田仁美,12650,2017
7,しばいぬ子,しばいぬ子さん,佐々木未来,1000,2017
7,日奈森亜夢,しゅごキャラ!,伊藤かな恵,1000,2017
7,宮藤芳佳,ストライクウィッチーズ,福圓美里,13789,2017
7,佐倉羽音,ばくおん!!,上田麗奈,2826,2017
7,袴田ひなた,ロウきゅーぶ！,小倉唯,6532,2017
7,湯音,異国迷路のクロワーゼ,東山奈央,1234,2017
7,夕立,艦隊これくしょん-艦これ-,谷邊由美,18031,2017
7,鎧塚みぞれ,響け!ユーフォニアム,種粼敦美,7854,2017
7,薄墨初美,咲-Saki-,辻あゆみ,4939,2017
7,たま,魔法少女育成計画,西明日香,1329,2017

100年間の気温変化

こんな話を見かけた。
「昔の夏はもっと涼しかっただろ！」と思い、50年前の日本の気温を調べてみたところ、まさかの結果が…！ | netgeek
これについて言及していた人がいて、データソースがおかしい、ということのようだ。
気象庁 | 日本の年平均気温
変化がない、と言ってるほうも気象庁を参照している、といっている。どこからデータをとってきたのだか…
これをみて、実際に気象庁からデータを取ってきてプロットしたという人がいた。

　
ここでよくあるのが、気温というデータは典型的な自己回帰がある時系列データなので、単純に線形回帰するのはどうなのか、という話。
　
自分は統計素人なので時系列データの解析は詳しくない。状態空間モデルでrstanを使って適当にやってみる。
あと、上のツイートのデータ、どこから取ってきたのかが分からない。最高気温が39度の年があるが、気象庁のどのリンクをたどって取得できたのかがわからない。
気象庁の　ホーム > 各種データ・資料 > 過去の気象データ検索　から、東京・東京を選択して、8月の100年間の平均・最高・最低気温データを使うことにした。
データリンク
　
モデルとしては、平均・最高・最低気温の状態は2階差分で正規分布で決まり、実際の観測気温は適当にコーシー分布からサンプリングした。
3つの気温状態のサンプリングに使う分散パラメータは、平均・最高・最低気温で共通としていたが の収束が悪く、3つ別々にしたら改善したので結局別々にしている。
　
結果としては、1910年までは最高・平均気温の変化はほとんどないが、最低気温は常に上昇している傾向である。
1920年から1940年頃は、最高・平均気温が上昇している。最高気温については、1940年代移行、横ばいであるが、2000年移行また上昇している。
平均気温は緩やかに上昇傾向。
最低気温は最高・平均気温に比べてかなり上昇幅が大きいが、2010年くらいから頭打ち、なのかもしれない。

　

# モデル
data{
  int N_year;
  vector<lower=0, upper=40>[3] Temp[N_year];
}
parameters{
  vector<lower=0, upper=40>[3] S[N_year]; // state 
  real<lower=0> s_s[3]; // state variance
  real<lower=0> s_o[3]; // observation variance
}
model{
  for(i in 3:N_year){
    for(j in 1:3){
      S[i, j] ~ normal(2*S[i-1, j] - S[i-2, j], s_s[j]);
    }
  }
  for(i in 1:3){
    Temp[, i] ~ cauchy(S[, i], s_o[i]);
  }
}

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

dat <- read.table("temperature.txt", header=F, row=1)
m02 <- stan_model("tokyo02.stan")
standata <- list(N_year=nrow(dat), Temp=dat)
fit <- sampling(m02, standata, iter=1500)
ex <- extract(fit)

s_s <- apply(ex$s_s, 2, median)
s_o <- apply(ex$s_o, 2, median)
S <- mapply(function(z) apply(ex$S[,,z], 2, median), 1:3)
cols <- c("black", "red", "blue")
year <- as.numeric(rownames(dat))
yl <- range(dat)
par(mar=c(4, 4.5, 2, 2), cex.lab=1.5)
matplot(year, S, type="l", col=cols, lwd=3, las=1, lty=1, ylim=yl, xlab="year", ylab="temperature")
for(i in 1:3) points(year, dat[,i], col=cols[i], pch=16, cex=0.5)
abline(h=0:40, lty=3, col=grey(0.7))
abline(v=seq(1800, 2030, by=10), lty=3, col=grey(0.7))
# legend("bottomright", legend=sprintf("%s気温", c("平均", "最高", "最低")), col=cols, pch=15, cex=1.4)

1876 26.6 31.6 21
1877 25.9 30.5 21.1
1878 24.6 28.5 21.2
1879 26.6 31.1 22.5
1880 25.5 29.5 21.4
1881 26.7 31.1 22.7
1882 24.7 28.5 21.2
1883 25.1 29.3 21
1884 24.1 28.1 20.4
1885 25.4 29 22
1886 26.5 31.2 22.5
1887 25.3 29.3 22.2
1888 25.6 29.7 22.2
1889 25.8 30 22.3
1890 25.4 29 22.8
1891 25.5 29.6 22
1892 26.4 31.1 22.5
1893 26.2 30.6 22.7
1894 27 31.5 23.4
1895 25.5 29.7 22.8
1896 25.9 30.4 22.5
1897 25 29.2 21.7
1898 26.1 30.6 22.7
1899 26.1 30.2 22.8
1900 26.1 30.5 22.6
1901 25.1 29.9 21.6
1902 22.9 26.2 20.2
1903 25.7 31 21.6
1904 25.1 29.9 21.2
1905 22.2 25.7 19.4
1906 24.5 28.5 21.1
1907 25.8 29.6 23
1908 25.4 29.5 22.3
1909 25.2 29.7 21.7
1910 24.1 27.8 21.2
1911 25.6 29.7 22.7
1912 25.2 30.1 21.6
1913 23.8 28.8 20
1914 26.4 31.3 22.7
1915 25.7 29.7 22.9
1916 25 29.2 22.1
1917 25 29.4 21.5
1918 26.1 30.6 22.5
1919 25 29.3 21.5
1920 25.7 29.3 23
1921 25.3 29.2 22.2
1922 27.3 32.2 23.6
1923 27.2 31.7 23.6
1924 26.2 30.7 22.4
1925 25.7 29.7 22.5
1926 26.3 30.8 22.5
1927 26.6 31.4 22.8
1928 24.1 27.7 21
1929 27.1 31.4 23.2
1930 26.8 30.6 23.4
1931 26.4 30.4 22.8
1932 26.7 30.6 23.6
1933 27.5 31.6 24.1
1934 26.1 30.4 22.5
1935 24.8 28.2 21.9
1936 26.5 30.8 22.8
1937 28.2 32.8 24.3
1938 26 29.7 23.2
1939 25.8 29.8 22.4
1940 24.9 29 21.4
1941 25.4 29.3 22.5
1942 27 31.3 23.5
1943 27.4 31.9 23.7
1944 27.5 32.1 23.9
1945 26.7 31.3 23.1
1946 26.7 31.4 23.3
1947 28 33.3 23.9
1948 25.4 29.7 22.4
1949 26.6 31.4 22.9
1950 26.2 30.4 23.2
1951 26.7 31.6 23.5
1952 26.8 31.4 23.7
1953 25 28.9 22.4
1954 27 31 24.2
1955 26.3 30.2 23.3
1956 25.4 29.7 22.2
1957 27.3 31.2 24.6
1958 25.8 29.7 23
1959 26.7 30.7 23.7
1960 26.4 30.4 23.4
1961 26.8 31.2 23.7
1962 28.1 32.9 24.4
1963 26.6 30.9 23.6
1964 27.8 32.2 24.3
1965 26.7 30.8 23.1
1966 26.9 30.5 24
1967 28 32 24.8
1968 26.6 30.2 23.4
1969 27.2 31.4 23.8
1970 27.4 31.7 24
1971 26.7 30.6 23.7
1972 26.6 30.7 23.3
1973 28.5 32.7 25.5
1974 27.1 31.1 24.1
1975 27.3 31.5 24
1976 25.1 29.1 21.9
1977 25 28.4 22.3
1978 28.9 33 25.6
1979 27.4 31 24.6
1980 23.4 26.6 20.7
1981 26.2 30 23.4
1982 27.1 30.2 24.6
1983 27.5 31 24.6
1984 28.6 32.4 25.7
1985 27.9 31.6 25.1
1986 26.8 30.4 24.1
1987 27.3 30.8 24.3
1988 27 30.2 24.4
1989 27.1 30.5 24.5
1990 28.6 32.4 25.7
1991 25.5 29.1 22.6
1992 27 30.4 24.3
1993 24.8 28 22.3
1994 28.9 32.9 25.9
1995 29.4 33.7 26.1
1996 26 30 23
1997 27 30.7 23.9
1998 27.2 31 24.4
1999 28.5 32.3 25.8
2000 28.3 32.4 25.4
2001 26.4 30 23.6
2002 28 32.1 24.8
2003 26 29.5 23.1
2004 27.2 31 24.1
2005 28.1 31.8 25.1
2006 27.5 31.1 24.8
2007 29 33 25.9
2008 26.8 30.7 24.1
2009 26.6 30.1 23.8
2010 29.6 33.5 27
2011 27.5 31.2 24.6
2012 29.1 33.1 26.3
2013 29.2 33.2 26
2014 27.7 31.2 24.8
2015 26.7 30.5 23.9
2016 27.1 31.6 23.9

賭ケグルイの投票じゃんけんで蛇喰さんはどうやって芽亜里の奴隷であるクラスメイトの割合を推定したか

賭ケグルイを見ている。
黒髪ロングのプリーツスカート黒タイツのはやみんなので視聴意欲がやばい。
　
ここで、賭ケグルイの1話で、投票じゃんけんという変則じゃんけんで勝負する。
クラス 人がグー 、チョキ 、パー のいずれかの手を投票する。
蛇喰さんと芽亜里はその投票箱から重複なしで3枚カードを引く。
各々の持ち手はその3枚となり、あとはじゃんけんのルールに従って勝敗が決まる。出した手は捨てられる。
　
ということである。ここで、両者120万円(120枚のチップ)から始めて、出した手の順に(他はアニメ演出上、蛇喰さんが把握できた手について)
1回戦：2枚賭け　蛇喰勝(G??)　芽亜里(C??)　アニメ絵では残りの2枚はG であることが分かっているが蛇喰さんには見えない
2回戦：50枚賭け　蛇喰負(GCC)　芽亜里(GG?)
3回戦：2枚賭け　蛇喰勝(P??)　芽亜里(G??) アニメ絵では残りの2枚はC と(たぶん) G であることが分かっているが蛇喰さんには見えない。
4回戦：50枚賭け　蛇喰負(PG?)　芽亜里(PP?)
その後、じゃんけんぽん、が3回は行われて、蛇喰さんがすべてのチップを奪われるが、即金で
N回戦：1000万賭け　蛇喰勝(C??)　芽亜里(PPG)
というふうになって最終的に蛇喰さんが勝利する。
　
ここでこの投票じゃんけんイカサマのネタバレだが、芽亜里はクラスメイト30人中21人の投票を操り、主人公の鈴井くんの提示した手(N回戦ではチョキ) 以外の手(つまりグー とパー)のどちらかに投票する。こうすると、クラスの の割合が投票を操作されているとき、 の投票確率は
となる。 ですべて、 のときグーとパーのみ出るようになって各々 である。
このとき、 が出る確率はかなり低いから、 を出しておけば少なくとも負けることはないだろう、というのが芽亜里の作戦である、たぶん。
　
ここで蛇喰さんは、2枚賭けのときは芽亜里が適当に手を出しているのに50枚賭けという大勝負のときには策を講じた手を出していること、自分の真後ろにたっていた主人公の手を手鏡で観察していたことから、N回戦のときに芽亜里にイカサマの存在を伝え、なおかつそれが、「10人から20人が芽亜里の指示通りに投票する」ということを見抜いた。
実際にこれを見抜けるかやってみる。
　
いま、 がある確率で出るとき(等確率 とする)、自分の手札の揃え方は10通りある。この10通りが出る確率は多項分布であり、以下の通りである。

x <- expand.grid(0:3, 0:3, 0:3)
x <- x[rowSums(x) == 3,]
colnames(x) <- c("G", "C", "P")
prob <- c(1, 1, 1)/3
dm <- apply(x, 1, dmultinom, prob=prob)
cbind(x, dm)

   G C P         dm
4  3 0 0 0.03703704
7  2 1 0 0.11111111
10 1 2 0 0.11111111
13 0 3 0 0.03703704
19 2 0 1 0.11111111
22 1 1 1 0.22222222
25 0 2 1 0.11111111
34 1 0 2 0.11111111
37 0 1 2 0.11111111
49 0 0 3 0.03703704

ここで、蛇喰さんがイカサマを見抜けたのは、「の手札が揃う勝負が極端に少ない」からという仮説にしておく。というのも、投票が偏らなければ投票箱は各々 になるはずだから、とする。たぶん。
ここで、上の各手の出方のうち、以外に絞れば、余事象的に「手札がすべて同じ(3枚)」と「2枚と1枚(と0枚)」は計算できる。
そしてこれはか、それ以外にすれば、二項分布で計算できる。
二項分布から95%信頼区間を計算してもよいが、蛇喰さんは「10人から、万全を期すなら20人」といっているので、適当にrstan を使って事後分布で評価してみる。
　
例えば、1回の勝負で が来たときの、投票箱の手の分布はこのような感じになる。
ここで、 の出る確率は足して1 なので、 の3次元平面になるが、平面なので2次元平面に写像した図にしている。三角形頂点のGCP の点に近いほど、その手の出る確率が1 に近くなる。

1回の勝負で2枚と1枚、例えば のときはこうなる。

1回の勝負で2枚と1枚、例えば のときはこうなる。

library(rstan)
code <- "
data{
  int trial;
  int gcp[3, trial];
}
parameters{
  simplex[3] a;
}
model{
  for(i in 1:trial)
    gcp[, i] ~ multinomial(a);
}

"
m <- stan_model(model_code=code)
trial <- 3
standata <- list(trial=trial, gcp=replicate(trial, c(1, 1, 1)))
fit <- sampling(m, standata, iter=30000)


M <- cbind(c(0,1,0.5), c(0,0,sqrt(3)/2)) # 2次元に写像する行列
triangle <- diag(1, 3) %*% M
xy <- ex$a %*% M
kd <- kde2d(xy[,1], xy[,2], c(bandwidth.nrd(xy[,1]), bandwidth.nrd(xy[,2])), n=500)

cols <- matlab::jet.colors(100)
i <- which(kd$z == max(kd$z), arr.ind=TRUE)
out <- c(kd$x[ i[1] ], kd$y[ i[2] ], 1) %*% solve(cbind(M, 1))
par(mar=c(3, 3, 2, 2), cex.axis=1.5)
image(kd$z, col=cols)
polygon(rbind(triangle[c(1:3,1),],c(-5,0),c(-5,5),c(5,5),c(5,-5),c(-5,-5),c(-5,0),triangle[1,]), col="white", border=NA)
polygon(triangle)
points(triangle, pch=16, cex=3, xpd=TRUE, col="hotpink")
text(triangle, c("G", "C", "P"), col="yellow", font=2, xpd=TRUE)
points(out %*% M, pch="★", col="hotpink", cex=1.5, xpd=TRUE)
legend("topright", legend=mapply(sprintf, "%s: %.4f", c("G", "C", "P"), out), bty="n", cex=2)

　
ここで大事なのが、確率的な問題なので、たとえ が出る確率が低くても(青の部分)、引く手が になることはありえる、ということである。ましてや、等確率 はそこそこありうる確率点である。
しかし、この分布も、勝負の回数が増えることで、等確率 はおかしい、という状況が出てくる。
　
さてここで、買収されている(奴隷)の割合を とすると、 が出る確率は、それ以外(2枚同じか、3枚とも同じ) 手が出る確率は となる。勝負が進んで蛇喰さんが手の内訳を確認できた、すなわち、蛇喰さん自身の手(これは勝負の回数分)と、芽亜里の手(初手で勝負が決まれば不明だが、2回目以降じゃんけんが進めばわかる)を足した回数、ということにする。つまり、アニメの4回戦が終わったときには、手を6回確認できていることとする。
芽亜里の奴隷の割合 が一様な事前分布を仮定したときの、N回手を確認したあとの事後分布は以下のようになる。


回が増えると分布は極端に、つまりみんな奴隷に傾く。
　
% 信用区間をプロットした。12回くらい手のうちがわかると、10/30 人が95% CI、20/30 人が50% CI、つまり中央値になるので、ここらへんのことを言っているのだろうか。

　
とはいっても、クラスメイト全員が芽亜里の奴隷であることはあるのはあるだろうが、ちょっと現実的ではない、かもしれない。このとき、奴隷割合 の事前分布をちょっといじってみる。
いま、一様分布を事前分布として仮定すると、も も等確率で生じてしまう。30人のクラスのうち、さすがに芽亜里が奴隷にできるのは半分もいかないだろう、という仮定をおけば、違う事前分布を設定できる。ここで、 は確率なので、適当にベータ分布からサンプリングしてみよう。ベータ分布のパラメータはふたつで決まる。ここで、ふたつのパラメータをどちらも1 にすれば、一様分布になる。
今回は適当に3と10を選んだ。おそらく あたりが最頻値になる。

　
これで同じようにすると、50回手のうちを確認してもたかだか 程度までしか事後確率は変動しない。

　
% CI をプロットしてみると、10/30 人の奴隷、20/30 人の奴隷を推定値として採用するのは難しいような気がする。というのも、50枚賭けを2回負けた時点で、持ちチップが20枚くらいしかないので、1枚ずつ賭けて負けても50回手のうちを確認することができないからである。
とすれば、 より多い割合で奴隷がいる、と事前分布を仮定するのが無難っぽい。というのも、転校してきて最初の勝負を自信満々にふっかけられたら何か裏があると思っておいたほうがいいっぽい。

　
統計ガチ勢の人がいたらぜひとも数式的にごにょごにょ解いてください。

code <- "
data{
  int n;
}
parameters{
  real<lower=0, upper=1> p[n];
}
model{
  p ~ beta(3, 10);
  for(i in 1:n)
    i ~ binomial(i, (2+p[i])/3);
}
"

m <- stan_model(model_code=code)
standata <- list(n=50)
fit <- sampling(m, standata, iter=6000, warmup=1000)
ex <- extract(fit)
med <- colMeans(ex$a)
cint <- c(99.9, 99, 95, 90, 80, 75, 50)/100 # 信用区間
ci <- apply(ex$p, 2, quantile, (1-cint)/2)

# 累積確率分布
cm <- apply(ex$p, 2, sort)
par(mar=c(5, 5, 2, 2), cex.axis=1.5, cex.lab=1.5)
matplot(seq(0, 1, length=nrow(cm)), cm, type="l", col=matlab::jet.colors(ncol(cm)), lwd=3, lty=1, xlab="index", ylab="cumulative probability")
abline(0, 1, lty=3)

# 信用区間幅のプロット
meari <- 10:20
cols <- matlab::jet.colors(length(cint))
par(mar=c(5, 5, 2, 2), cex.axis=1.5, cex.lab=1.5)
matplot(t(ci), type="l", lwd=3, col=cols, lty=1, ylim=c(0, 1), xlab="カード内訳を確認した回数", ylab="クラス内の寝返り割合")
abline(h=meari/N, lty=3)
abline(h=n0/N)
legend("topleft", legend=sprintf("%.1f%s", cint*100, "%"), title="CI", lty=1, lwd=3, col=cols, bg="white")

# 分布のベタ書き
dxy <- mapply(function(z){
  dx <- density(unlist(z))
  i0 <- 0 < dx$x & dx$x < 1
  return(list(x=dx$x[i0], y=dx$y[i0]))
}, apply(ex$p, 2, list), SIMPLIFY=FALSE)

par(mar=c(5, 5, 2, 2), cex.axis=1.5, cex.lab=1.5)
yl <- range(0, dxy)
cols <- jet.colors(length(dxy))
plot(0, type="n", ylim=yl, xlim=c(0, 1), xlab="probability", ylab="density")
for(i in seq(dxy)){
  lines(dxy[[i]]$x, dxy[[i]]$y, col=cols[i], lwd=2)
}

糖質制限ダイエットを始めたらたった1日で体重が2kg減った話をしたらもっと詳しい体重推移データをもらった

ダイエットのデータを昔使って遊んだけど、元データを持っている人が更に長期のデータを公開してくれた。

さっそく遊んでみる。
　
2013年は測っていない時期が多すぎて、推定がクソ。
以前遊んだときは、ダイエットを意識して計画的に体重を測定していた2016年のものだったらしい。2015年とかは体重が増えている。
2016年のダイエット意識高い系()以降は油断して体重がまた増えている。
　
ダイエットさぼったとかそういう異常検知系のことができませんかね? というディスカッションになったけど頑張ればできるかも?

single cell RNA-seq のdropout

読んだ。
MAGIC: A diffusion-based imputation method reveals gene-gene interactions in single-cell RNA-sequencing data
コードはPython で書かれている。

computational flowcytometry のDana Peer ラボ。RNA-seq のデータ行列が取るであろう高次元空間の形からデータ点(細胞)を復元しようという雰囲気。
single cell RNA-seq ではmRNA の読みの精度の悪さのせいで、全体の80-90% くらいは0 という結果が返ってくるdropout という現象に悩まされる。有効なデータが10% 程度しかなく、他が0 のときに元のデータがどのようなものだったかを推定するのが必要な作業だが、MAGIC(Markov Affinity-based Graph Imputation of Cells) という、マルコフ連鎖と近傍のデータ点からもともと読まれるはずだったmRNA のデータと細胞を再構成してデータを復元する方法があるらしい。
　

図はNat Methods. 2014 Jul; 11(7): 740–742.のFig.1a
　
n 遺伝子 m 細胞の行列 がある。このなかのほとんどはdropout のせいで0 であり、データを復元したがほしい。細胞で眺めれば遺伝子発現のベクトルになるので、このベクトルで表現型が分類でき、ベクトルが近いものはその近傍に細胞が固まるはずである(diffusion)。グラフ理論的にある表現型(ベクトルパターン)をもつ細胞の近傍k 個を考える。
からは簡単に距離行列が計算できる(適当にユークリッド)。そして、を適切なadaptive ガウスカーネル関数で変換するのだが、

というように変換してaffinity matrix を作る。このとき、 はkernel width と呼ばれ、このアルゴリズムでは非常に重要である。 以下のは速やかに0 affinity になるため、相当近い細胞同士でないとneighborhood が構成できない。
は対称行列でないので として強制的に対称行列化して、マルコフ過程を満たすために行和を1 になるようにする。つまり確率推移行列になるようにする。
ここで、 というべき乗パラメータを使って、 とすると、目的だったが得られる。ここでは、dropout して80-90% が0 だったRNA-seq データが、もともと得られるはずだったRNA-seq リードカウントデータになっている、ということらしい。
( は行列積ではなく、ただ単に要素のべき乗である)
　
scImpute: Accurate And Robust Imputation For Single Cell RNA-Seq Data
MAGIC とは異なり、発現の多寡とdropout するかしないかを確率分布からモデル化して、それを説明するパラメータを推定する。
データの取得のされ方について以下のようなモデルを考える。

：dropout rate
：ガンマ分布のパラメータ。ガンマ分布はdropout をモデル化しており、遺伝子発現としては高くてばらつきの少ないものだが、これが0になるのはdropout によるもの、という感じ。
：正規分布のパラメータ。正規分布は通常の遺伝子発現をモデル化しており、低めから中程度で、ばらつきの大きな遺伝子発現をしており、これが0 になるということはまあ普通に考えて遺伝子発現がないから0 である、という感じ。
これをEM アルゴリズムでゴリゴリ解く。ベイズでやってもいいと思う。
実際に 遺伝子が 細胞でdropout する確率は

となる(各パラメータはEM アルゴリズムで推定されたあとのハット付きである)
　
dropout 確率と、細胞遺伝子発現の隠れモデルが推定されたので、実際にdropout しているデータをimputation する。ある適当なdropout 確率閾値 により、細胞 についてimputation されないといけない遺伝子発現セット と、dropout 確率閾値を下回っているため確実にデータとして信頼できる遺伝子発現セット に分けて、 を用いてimputation する。これは 正則化で行い、

となる発現量を復元する。
dropout するかしないかというのは、例えば行動定量学で鳥の生態を観察するような、鳥は生きているが観察されないのか、鳥は死んでいるので観察されないのか、というのをモデル化するzero inflated model 的な感じがした。
　
Nat Methods. 2014 Jul; 11(7): 740–742.
もともとはこれがいっぱい引用されている。
ベイズ的にやるのだが、モデル式としては、dropout はポアソン分布、amplification は負の二項分布からサンプリングされるとしている。負の二項分布を使うと裾が伸びて分散が大きい値を表現できるから、これが高発現というかamplification を意図しているのだと思う。
あるsubpopulation で平均して の発現量がある遺伝子の事後確率を

とする。ただし、 は のブートストラップサンプル、 はRPM (RNA-seq で得られるリードカウントを適当に補正したようなもので、なんらかの定量値と思えばいい)、 はエラー項であり、 の細胞について、 は事後確率であり、dropout を観測する確率 を用いて

と書く。
というようなことをやっているので、結局はdropout で0 になる、その0 になるのが通常発現パターンか、高発現パターンなのかを別のそれっぽい分布から出てきたようにモデル化するっぽい。

糖質制限ダイエットを始めたらたった1日で体重が2kg減った

親戚が集まったときの会で、ふいにそんなことを言われた。その人としては減量をがんばったということを言いたかったのだろうが、自分自身、糖質制限ダイエットに対してエビデンスを持っていないことと、1日で2kg 減るっていうのは、偶然「2日間連続で体重を測ったら2kg 減っていた」というポジティブな観測を強烈に覚えているだけの観測バイアスなので「水分だけでも1〜2kg 程度は体重の日内変動はある。1ヶ月間毎日体重測ってプロットしてから出直せ」と返してしまったのだが、本当に糖質制限ダイエットをするだけで1日で2kg 減ったらどうしよう、と思ったので、データで考えてみる。
　
ある程度の期間、それなりに管理された測定データが欲しいのだが、実は某氏の体重データがある。某氏はデータ取得と解析の知識があるので、データのとり方は「毎朝、起床時に自宅の体重計で測定する」というものらしい。出張などで測れないときは欠測である。170cm 半ばの某氏は、さすがにラーメンを食べ過ぎたので減量しようと思ったようである。ただし、食生活を変えたり、運動習慣を強化した、というようなダイエット行動がいまいちわからなかったのでよくわからない。
　
ここで、体重を測っていなかったとしても、決して体重が存在しなかった、というわけではない。観測していなかっただけである。そういう場合に、状態空間モデルを使って、「真の体重」というのもは常にあるが、それを観測するときに誤差が生じるし、場合によっては欠測にもなる、というのをstan でやる。
　
真の体重 は、ふたつ前までの状態の差分で決まるとする。このとき、差分に対して適当なパラメータが決まり、その分散は固定された という実数とすると

つまり

である。
実際に観測した体重 (大文字) は、測定時に着ている服とか、体内の水分バランスとか、体重計のキャリブレーション具合とかで、なんらかの観測誤差を含むと考えられるので、適当な分散 を用いて

でサンプリングされる、とする。毎回の観測の条件を揃えると、 は小さくて済む。
欠測は、rstan ではNA を受け付けないので、欠測の日を先にindex で作成しておいて、rstan のmodel ブロック内でfor loop で回避する。
　
, をcauchy 分布からサンプリングしているが、なるべくばらついて大げさな値が出て欲しかったのでやってみた。適当である。
　
でうまく収束したようだ。計算時間は1000 iter で1 chain 3分程度。

1日での体重減少の事後分布を見てみると、中央値で0.03 kg 程度減少する。95% 信用区間では、実は有意に体重が減少するという結果ではないようだ。しかし、200日程度のデータがあるので、分布の正の領域と負の領域が半々でない以上、負の方向(痩せる方向)へ傾いていくのは妥当だろう。

quantile(c(apply(ex$w, 1, diff)), c(a/2, 0.5, 1-a/2))

       2.5%         50%       97.5% 
-0.08474003 -0.03705018  0.01160505 

何も考えず、データ観測期間内での体重減少から日割りすると、-0.04 kg 程度だった。

diff(range(dat$weight, na.rm=TRUE))/nrow(dat)

[1] 0.04056604

　
さて、ここで気になるのが、「1日で2kg 痩せる」というのがどれくらいの出来事なのだろうか、ということである。体重の減少は、「真の体重減少」に加えて、そのときの偶然の誤差を加味した「体重計の指し示す値」で決まる。ここでいう偶然の誤差は、先に言った通り、体重計のキャリブレーションや、そのときの服装、直前に食事をしたり排泄をしたかどうか、が含まれる。
この「誤差」の作用は、上のモデルの によって表現される。ここから、2kg 減ったというのが95% 信頼区間(ここでは頻度主義的に) でどうか、というのを調べるには

qnorm(1-a/2, sd=s)

[1] 0.8670893

となり、0.87kg 以上の観測誤差があるとき、2kg の1日減量やばい、ということになる。よって、糖質制限ダイエットの真の1日体重減少量は1.2 kg 程度となる。
(体重の減少、が興味の対象なので、片側では? という話もあるが、このときは0.72 kg である)

単純に、糖質制限ダイエットの真の減量効果が1kg/day としよう。1ヶ月で30kg 痩せる。1ヶ月で30kg 痩せたら自分としては悪性腫瘍を疑って医療機関の受診を勧める。
　
というのはまあ笑い話なのだが、糖質制限はするがたんぱく質、脂質はどんだけ摂取してもよいらしい。糖質がないと脳内でエネルギーにはならないけど? と聞くと、午前中はツライらしく、脂質のエネルギーって9kcal/g だけど? と聞くと、脂質は肌の恒常性の維持に必要だからむしろたくさん摂取しろという。謎()
　

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
dat <- read.csv("weight.csv")


code <- "
data{
  int Day;
  real<lower=0, upper=200> W[Day];
  int Na;
  int na[Na];
}
parameters{
  real<lower=0, upper=200> w[Day]; # weight
  real<lower=0> b; # sigma for hidden model
  real<lower=0> s; # sigma for observation
}
model{
  for(i in 3:Day){
    w[i] ~ normal(2*w[i-1] - w[i-2], b);
  }
  for(i in 1:Na){
    W[ na[i] ] ~ normal(w[ na[i] ], s);
  }
  b ~ cauchy(0, 0.5);
  s ~ cauchy(0, 0.5);
}
"

m <- stan_model(model_code=code)

W = replace(dat$weight, is.na(dat$weight), 1)
standata <- list(Day=nrow(dat), W=W, Na=sum(!is.na(dat$weight)), na=which(!is.na(dat$weight)))

fit <- sampling(m, standata)

ex <- extract(fit)
a <- 0.05
w <- apply(ex$w, 2, median)
b <- median(ex$b)
s <- median(ex$s)
ci <- apply(ex$w, 2, quantile, c(a/2, 1-a/2))

par(mar=c(5, 4.5, 2, 2), cex.lab=1.5, cex.axis=1.5)
xi <- seq(nrow(dat))
yl <- range(dat$weight, ci, na.rm=TRUE)
plot(w, type="l", ylim=yl, xlab="Day", ylab="Weight [kg]")
polygon(c(xi, rev(xi)), c(ci[1,], rev(ci[2,])), col=grey(0.8), border=NA)
lines(w)
points(dat$weight, pch=16, col=4, cex=0.8)
points(which(is.na(dat$weight)), w[which(is.na(dat$weight))], pch=16, col=2)
legend("topright", legend=c("実測値", "欠測時"), col=c(4, 2), pch=16, cex=3)