So-net無料ブログ作成
  • ブログをはじめる
  • ログイン

RでJAGSとやらを使うつもり [統計学もどき]

以前に検索で見つかったこちらのPDFから
http://www.rri.kyoto-u.ac.jp/NSRG/FSEB/abstracts2015-8.pdf

村瀬香・名古屋市立大学大学院システム自然科学研究科

階層ベイズと統計モデリングを通じて、原発事故が野生動物とヒトに与えた影響を考える

という発表が気になりまして、これらを勉強するために

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)   久保 拓弥 (著)

というのが定評があり、幸い県立図書館に置いてましたのでリクエストして借りました。
なかなか簡単に読み進めませんが、階層ベイズモデルを使うにはRだけではだめで、WinBUGSというのを使わないといけないと書いてありました。一応インストだけしてみるかとチャレンジしたところ、いまどきパッチをあてるなどと言うことを要求したり、Win7で使うにはいろいろややこしいことを要求したのであきらめかけたところ、村瀬香氏の論文を探すと#原子力発電_原爆の子さんブログに翻訳があり、ここを見ると分析に使ったソフトが書いていました。

http://besobernow-yuima.blogspot.jp/2015/04/httpt_3.html
#Nature 誌サイエンティフィック・リポーツ「オオタカの繁殖に対する福島第一原発事故の影響」

「われわれは、R 3.0.2上でJAGS 3.3.0とrjags 3–12を用いてMCMCサンプリングを実施した。サンプリング過程で、最初の10,000回試行はしくじりとして破棄され、MCMC連鎖ごとに100回試行をサンプリングとして、200,000回試行を実施した。5連鎖に対して、同じ手順が適用され、10,000件のサンプルが生成された。われわれは、収束兆候を求め、また値が1.0に非常に近接することを確かめるために潜在スケール減少係数R30の多変数版を計算した 30。」

JAGSというのがWinBUGSの代わりで、rjagsがRからJAGSを動かすものとわかりました。
久保 拓弥氏のHPにちゃんとインストの仕方が出ておりまして

生態学のデータ解析 - R で JAGS (rjags 編)
http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoJags.html

一応 RやEZRから JAGSが使えるようになったはずです。サンプルを試してみたけどうまくいかずあきらめかけていたところ、こちらにあったサンプルは

JAGSを使ってギブスサンプリングを試してみた (1)
http://www.singularpoint.org/blog/r/mcmc-jags-1/
JAGSサンプル.png
意味はわかっていませんが、動いたみたい。

ちょっと画像をアップしたいので、コメントでなく追記の形で書かせていただきます。

「道具としてのベイズ統計」にも出てくるのですが、尤度をポアソン分布にして、事前分布をガンマ分布にして掛け合わせると、事後分布はガンマ分布というおまじないみたいな言葉ですけど、まあなにか実例で試したいと思って探したら次の所が見つかりました。

http://rmecab.jp/wiki/index.php?Bayes_Poisson

ポアソンの共役事前分布はガンマ分布とする

Rのスクリプトはエラーなく動いたのですが、グラフが変になってしまいました。
ポアソン共役分布.png

原因はなんだったのでしょう。

ポアソン回帰
R_MCMC.png
一般線形化モデル
R_GML.png

は上の数行をコピペするだけでうまくいきましたが、意味がほとんどわかっておりません。R修行もやっぱせんといかんのかな。
nice!(1)  コメント(24)  トラックバック(0) 
共通テーマ:blog

nice! 1

コメント 24

ryu-ron

ギブスサンプリングという言葉は、道具としてのベイズ統計で一応読んではいるのですが、まだまだ難解ですね。

結局 3級の勉強は進まず、どうも秋に3級、来年6月に統計調査士となりそうな気配です。

なんとなく先にベイズが気になりだして、2級や準一級に出てくる問題の知識があってもわかりそうにないことに関心を持ったものだから、検定試験への執念がダウンしつつある。これではいかんと、3、4級問題集に再び立ち向かう日を決めていかねばなりませんが。
by ryu-ron (2016-05-01 10:21) 

ryu-ron

大学のプレスは残念ながら削除されてる模様。
http://www.nagoya-cu.ac.jp/secure/149568/270324.pdf

この記事が刺激したのか?
http://www.huffingtonpost.jp/science-portal/falcon_b_6969272.html?utm_hp_ref=tw

0.1uSv/hの空間線量の上昇は、最大で 10%の繁殖成功率の低下に寄与していると試算

これが論文中から読みとれないのだが。
by ryu-ron (2016-05-01 12:02) 

ryu-ron

こっちか
http://www.nature.com/articles/srep09405
グーグル翻訳しますと、「討論」の前に

次の質問は、空間線量率が巣の成功に影響を与えている程度でした。事後平均使用してβ 2、-3.909(表5)、およびprequake期間中の巣の成功は繁殖地の中で最も高かったサイトの第10号のための巣の成功の事後分布を(図3b、補足表S1)、巣の成功の最大の低下は、空間線量率(0.1マイクロシーベルト/ hの増加と0.100と算出された補足図S5)。

ああ井上さんの訳でも「考察」の前ですね。

次の問題は、空間線量率が繁殖成功率に影響をおよぼす程度だった。β2の事後平均値[−3.909](表5)、それに繁殖地のうち、震災前期間における繁殖成功率が最高だった繁殖地No. 10の繁殖成功率・事後分布(図3b、補足表S1)を用いて、繁殖成功率低下の最大幅を、空間線量率の増加幅0.1 μSv/hあたり0.100と計算した(補足図S5)。

これでした。
by ryu-ron (2016-05-01 12:09) 

ryu-ron

今月は福島県の県民健康調査の発表もあります。まだ日程は公表されていませんが、増加がゼロでなければ、やはりそろそろ因果関係を認めるべき時期が近付いているのではないか?

原爆の初期被ばくで、短命の核種の影響が言われてましたが、今ヨウ素が少ないと言っても、事故当時のヨウ素量は半端出ない。

事故当初の内部被曝量なんて測定できていない以上、ごたくを並べず、はっきりと事故の関係を調査すべきだ。
by ryu-ron (2016-05-02 10:01) 

ryu-ron

長すぎて裁判所がどう判断したかは読みとれてないのですが
http://www.courts.go.jp/app/files/hanrei_jp/706/037706_hanrei.pdf

139ページに

9「原爆放射線の人体影響1992」350頁によれば,経過時間による線量率への寄与は,約30分後からマンガン56(2.6時間)とナトリウム24(15時間)が,約1週間後からは鉄56(44日)とスカンジウム46(83日)が,約1年後からはマンガン54(312日)とセシウム134(2.05年)がそれぞ
れ主であるとされており(かっこ内は半減期。以下同じ。),DS86報告書第6章(乙A16)では,これらのほかにコバルト60(5.3年)が挙げられ,アラカワによる「広島および長崎における残留放射能」(乙A11)では,花崗岩標本の分析に基づき,熱中性子によって広島の土壌中に発生した最も重要なガンマ線放出同位元素として硅素31(半減期2.65時間)等が挙げられている(8頁)。
よって,被告らの引用が半減期の極端に短い元素を抽出した恣意的なものであることは明らかである。また,マンガン56やナトリウム24は,減りやすい反面,その間に放出する単位時間のガンマ線の放出量が大きいため,早期に爆心地付近に入った者は,マンガン56やナトリウム24が急速にガンマ線を放出しつつある時期に被曝し,比較的遅く爆心地付近に入った者については,被告らが意図的に例示から除外した比較的長半減期の同位体に由来する誘導放射線の寄与を無視できないのである。

短いものは「単位時間のガンマ線の放出量が大きい」というところが重要ですね。
by ryu-ron (2016-05-02 14:31) 

ryu-ron

前の方のコメントで紹介した
http://www.hiroshimapeacemedia.jp/?p=11675
ですが


マンガンは土中に豊富に含まれる物質。原爆で放射された中性子をマンガンの原子核が吸収すると、放射性マンガンに変わる。

 これまで候補に挙がっていたナトリウムやスカンジウムは半減期がそれぞれ15時間、83・8日と長く、6日入市者だけリスクが高い理由を説明できないため、候補から外れた。

こういう絞り込みができれば、裁判の証拠として活きてくるのではないでしょうか。
by ryu-ron (2016-05-02 15:29) 

ryu-ron

アラカワによる「広島および長崎における残留放射能」(乙A11)では,花崗岩標本の分析に基づき,熱中性子によって広島の土壌中に発生した最も重要なガンマ線放出同位元素として硅素31(半減期2.65時間)等が挙げられている(8頁)。

ここに出てきた硅素31がひょっとして、福島のセシウムボールにも含まれていなかったが気になったのですが、ケイ素31Si31は崩壊するとリン31 P31に変わって、そのあとはこれで安定するとありますので、いろいろ探しましたがp31が検出された様子はないので、まあ大丈夫なのでしょう。

それとこちらの判決書では
http://www.courts.go.jp/app/files/hanrei_jp/125/082125_hanrei.pdf
53ページに
「現在の線量推定に関連があると考えられる放射性核種は,アルミニ
ウム28(半減期2分),マンガン56(半減期2.6時間),ケ
イ素31(半減期2.6時間),カリウム42(半減期12.4時
間),ナトリウム24(半減期15時間),鉄59(半減期44.
5時間),スカンジウム46(半減期83.8時間),セシウム1
34(半減期2.1年),コバルト60(半減期5.3年)である。
このうち,ケイ素31,カリウム42及び鉄59は重要でない。」

とありますので、裁判所としては短命核種ケイ素31を問題視してないようです。

by ryu-ron (2016-05-06 08:53) 

ryu-ron

アラカワによる「広島および長崎における残留放射能」(乙A11)とは
http://www.rerf.or.jp/library/scidata/tr_all/TR1962-02.pdf
これのようですね。

14ページ目が、8ページで表にSi30,Si31と出てきますね。

by ryu-ron (2016-05-06 09:03) 

ryu-ron

なに規制委員会のサイトにある
http://www.nsr.go.jp/archive/jnes/content/000017127.pdf

ベイズ用語がずらずら出てくる。福島のデータ計算したら因果関係わかるでしょ。使い方知らんのか?
by ryu-ron (2016-05-06 17:01) 

ryu-ron

オオタカ論文に対するツイッターのまとめ
http://togetter.com/li/804315?page=3

あ またこの人か。見るのやめた。
by ryu-ron (2016-05-06 17:19) 

ryu-ron

原子力安全基盤機構の生涯がんリスク解析の報告書は、RやWinBUGUなんかを使っているので勉強にはなりますが、高線量ひばくのパターンなので、これを低線量の場合で考え直すのはなかなか大変なようです。

最も読むのもかなり大変なページ数ですが。

さて、データ解析のための統計モデリング入門もだいたいコピーをとったので、そろそろ返して、来週は仕事に励みます。月曜からハードになりそう。
by ryu-ron (2016-05-07 09:54) 

ryu-ron

ここに階層ベイズが出てきますね。
http://rerf.or.jp/intro/report/planj(2016).pdf
2016 年(平成 28 年)度事業計画

放影研ですか。

5ページ

被爆位置:被爆位置、線量推定における円形対称性等の検定の前提条件、両市の広い地域におけるバックグラウンド率の空間的均一性などに関する既存の二次元データを最大限に活用するため、地理空間的手法を放影研データに適用する方法の開発を目的とする調査に着手する。このため、階層ベイズ法や経験ベイズ法などの最新の方法を用いる。
原子爆弾から直接受けた放射線量や他の共変量の影響が強いので放影研データは空間解析に関し難があり、その影響についてはシミュレーション解析の一環として同一のデータから推定する必要がある。最終的には、特定のプロジェクトに急性脱毛、染色体異常、およびがん罹患率に関する空間解析を含めるかもしれない。
by ryu-ron (2016-05-07 11:12) 

ryu-ron

いろいろ探しているとこういう論文が
https://www.jstage.jst.go.jp/article/mammalianscience/54/1/54_53/_pdf

状態空間モデルを用いた階層ベイズ推定法によるキョン(Muntiacus reevesi)の個体数推定

科研DBで見た記憶はあったのですが、オオタカ論文に似た手法で分析してますね。

キョンといえば、八丈島ではなかったか?これはこまわりくんの世界だったか。今の若い人は知らないでしょう。

摘     要
千葉県房総半島に生息するキョン(Muntiacus reevesi)
の 2006年度~ 2011 年度末時点における 59 管理ユニット
の個体数について,糞粒法,区画法の調査結果およびユ
ニット別,年別捕獲数を用いて状態空間モデルを構築し,
階層ベイズ推定法で推定した.推定の結果,2011 年度末
の合計個体数は中央値 19,826 頭(95%信用区間:14,542
~ 26,422 頭)となった.従来の糞粒区画法と出生数捕獲
数法による総個体数推定値は,やや過少評価していた.
ベイズ法を利用することで,より推定幅が狭く精度の高
い推定が可能となった.個体群増加率は平均 1.294であっ
た.野外で捕獲された個体の性齢構成と妊娠率から,年
1 回の繁殖を仮定すると,個体群増加率は 1.356 ないし
1.407 となることから,無視できない程度に大きい死亡率
が個体群動態に寄与していることが推測された.

両面印刷をしてしっかり読みます。WinBUGSを使えるようにせんといかんのか、読むだけにしておきましょう。
by ryu-ron (2016-05-07 11:15) 

ryu-ron

きょん論文は、原発事故との関係はいっさい触れていませんが、11年度末(12年3月)時点のデータで死亡率が上がったことを示しています。
その後、4,5年のデータを追い掛けた場合にどんな結果が出ているか。
気になるところです。
by ryu-ron (2016-05-07 11:26) 

ryu-ron

あれま事故後に異常繁殖してる。
http://futonsend.net/entry/chiba-kiong-breeding-over-40000

こりゃ参考にならない。まあベイズの勉強用ですね。
by ryu-ron (2016-05-07 11:31) 

ryu-ron

道具としてのベイズ統計は、連休中も持ち歩いていたのですが、今一つ読めておりませんで、あいかわらずふとん専用です。
MCMCや階層ベイズがエクセルでも試せるというので、出版社のサイトからエクセルサンプルファイルをダウンロードはしております。

一応これを最後まで読みあげてから、サンプルを動かしてみますかね。

きょんの論文も読んでは見ましたが、準備段階の大変さがとんでもなくて、まあ、同じようなことができるようになるには何年かかるかわかりません。
by ryu-ron (2016-05-07 16:53) 

ryu-ron

MCMCのわかりやすいところはないかと探すと
http://tjo.hatenablog.com/entry/2014/02/08/173324

みどりぼんの解説になるのですが、昨日の夜iPadで読んだときはわかったような気になったが、今朝、パソコンで読みなおすと、やっぱ難しい。

まあ繰り返しですね。
by ryu-ron (2016-05-08 09:43) 

ryu-ron

ちとどうでもいいことながら統計ネタ。
おととい紀伊国屋に行ったときに「身につく ベイズ統計学 (ファーストブックSTEP)」しかと眺めてみたのですが、やはり1冊の本で、階層ベイズやMCMCまでを解説したものとしてはなかなか優れていますね。
税込みで2000円は越えてしまうけど、コストパフォーマンスはなかなかいいです。
道具としてのベイズ統計のわかりにくいところを、図解・ベイズ統計「超」入門でカバーしながらやってます。
これ涌井兄弟、別々に書いていたのですね。”身につく”のほうは、二人の共著です。
by ryu-ron (2016-05-10 09:17) 

ryu-ron

ちょっと追記をしました。
津田先生もポアソン分布は使われるので、ベイズでひねりを入れるとどうなるか試したくて、自然な共役分布なる手法を使おうとしたけど、さっぱりわからん。
Rのサンプルはあったものの、グラフがうまく表示されず、なにかパッケージを読みこまないとだめだったのかな。
by ryu-ron (2016-05-13 08:37) 

ryu-ron

ここの最初の方のコメントを見たら、統計調査士を6月とか書いてしまいましたが、これは11月にしか試験がないので、来年になりますね。

今度の6月の試験日から3級の過去問解きを始めましょう。
by ryu-ron (2016-05-15 13:14) 

ryu-ron

昨日 駅ビルの紀伊国屋にいけば
https://www.kinokuniya.co.jp/f/dsg-01-9784534047694
Excelでスッキリわかるベイズ統計入門
があるかと思ったのですが、在庫はありませんでした。

道具としてのベイズ統計の自然な共役分布の説明はそれなりに詳しいのですが、残念ながらエクセルのサンプルがない。

まあ買うつもりはないのですが、サンプルファイルのあるところのURLを暗記して帰ろうと思っていたのに残念です。
内容的に道具としてのと同じ部分が多いので、古本を買うこともないのですが、県立図書館にも置いてなくて、自分で考えるにはちと難しい。
by ryu-ron (2016-05-15 14:55) 

ryu-ron

さて 統計ネタをあとひとつだけ。道具ベイズは二回目の通読に入りました。今度は少し頭に入れて行かねばと思ってるのですが、みどりぼんのコピーを並行して読んでると、また違った角度から考えられるのでいい感じです。

ちょうど去年の今頃だったのか?探偵の探偵で北川景子の部屋に分厚いベイズ統計解析ハンドブックが置いてあったのを。

修行して探偵をめざすか。体力が追いつかない。あんな命がけの事ほんとにやってるのかな?
by ryu-ron (2016-05-17 08:51) 

ryu-ron

おっとあとひとつとしながらもうひとつ。
これは統計ネタではなく甲状腺がんですが、今月ある予定の県民健康調査が遅れています。

http://www.pref.fukushima.lg.jp/site/portal/kenkocyosa-kentoiinkai.html

去年は5月19日にはやっているので、年度が変わってあわただしいは通用しない。
まあ、まだ数日あるが、だいたい事前アナウンスがあってたはず。
by ryu-ron (2016-05-17 09:55) 

ryu-ron

ここにマンガンのことを書いていたのですが、以前の記事でJCO事故のときのマンガンを書いたのを思い出しました。
http://ryukei-rondan.blog.so-net.ne.jp/2016-01-18

この時ですね。kokikokiyaさんに中性子のことを指摘されると、ノーコメントと言わざるをえなくなるのですが、
短命のマンガン54は、
http://www.cnic.jp/knowledge/2586

軽水炉から放出されるマンガン-54

1970年代後半に、敦賀・福島第一・浜岡の各原発の周辺で採取された松葉からマンガン‐54とコバルト-60(60Co、5.27年)が検出され、話題になった。粉塵として放出され、松葉の上に沈積したのである。初期に建設された沸騰水型軽水炉の周辺では日常的に起こっていた。
原発作業員の体内にある放射能の測定によってマンガン‐54とコバルト-60が検出されている(森江信「原子炉被曝日記」、技術と人間社、1979年刊、pp.79-81)。定期検査の時などに体内に入ったのであろう。
一次冷却水の配管の近くで作業をすると、ガンマ線による被曝を受ける恐れがある。手や足での被曝線量がポケット線量計やフィルムバッジが装着されている胸部の線量の100倍に達する場合もあろう。原子炉の周辺での作業は、非常にきびしい条件下にあると想定できる。

まあ即発臨界でなくても出るのかな。

軌道電子を捕獲して、クロム-54(54Cr)となる。ガンマ線が放出される。

ということでクロムであればセシウムボールから出ているので、新しい記事に引っ越ししてなんかコメントするとしましょう。
by ryu-ron (2016-05-17 12:25) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

※ブログオーナーが承認したコメントのみ表示されます。

トラックバック 0