R作業は統計環境『R』を用いるものとする。(MCMCpack)を用いる(その使い方についてはこちら)サイコロ今、サイコロの「できのよさ」を調べながら、そのサイコロで博打をうつことを想定するサイコロであるからには、6つの目がほぼ同確率で出るようになっている。その確率はRで言えば、次のように6つの目にそれぞれ同確率をもつベクターとして作ることができるdice<-c(1/6,1/6,1/6,1/6,1/6,1/6)理想的なサイコロ今、サイコロを振って、何の目が出るかを予想している。あるでは、このように「理想的なサイコロ」を振って次に1から6のそれぞれの目が出る確率を予想しよう。正確に一致する数値を6つ返してくる(はずである。集中するかが図示できる。してやれば、tの目が出る確率の分布が得られる現実のサイコロでは、現実のサイコロで、1から6の目の出る確率が均一でないものとする。(足して100回)となったとする。この結果を見て、このサイコロは、2の目が一番でやすい傾向があり、4の目が一番でにくい傾向があるように思われるただし、それは、1の目と5の目の出る確率が、きっかり一緒だと予想することが正しくないだろうということや、それぞれの目の出る確率がそれぞれ正確に0.18,0.20,0.16,0.12,0.18,0.16であると予想するのが危険だろうことからもわかる通り、「おおまかな傾向」を示しているだけで、正確なことはわからないしかしながら、この100回の実験をしてみる前よりもなにがしかの情報が得られていることは確かであるこのような事前情報が得られているときに、次(101回目)に1から6の目のどれがどのくらいの確率で出ると上述の「理想のサイコロ」の伝で行くと次のコマンドで、1の目が出る確率をk回出してみた場合の分布が得られる。hist(rdirichlet(k,c(18,20,16,12,18,16))[,1],xlim=range(0,1))もちろん、1の目が出る確率の分布の平均値は0.18に近くなる。そのことは次のコマンドで確認できるmean(rdirichlet(100000,c(18,20,16,12,18,16))[,1])理想的なサイコロ、再び今、実は、「理想的」であるようなサイコロがあるとする。残念ながら、それが理想的なサイコロかどうかは人智では知る由もない。実際にこのサイコロのよさを調べるとなると、上述の「現実のサイコロ」と同様な実験をするしかない。集計結果が出たとする。その集計の度に、そのとき得られている情報でそのサイコロの目の出る確率に見当をつけ、次に出る目の確率を予想するとするとその作業は次のようにrdirichlet関数を使って示すことができる1の目が出る確率を、全60回の実験後、全600回の実験後、...、全600000回の実験後に予想して分布をとるとするとhist(rdirichlet(k,dice*100)[,1],xlim=range(0,1))hist(rdirichlet(k,dice*1000)[,1],xlim=range(0,1))hist(rdirichlet(k,dice*1000000)[,1],xlim=range(0,1))となる。』に近づいていくことがわかるより現実的にはかなり精度のよいサイコロであると信じるに足るサイコロを持っているとする。みるからに立方体だし、いかさま用に作られたとは思えない、小学2年生1月号の付録のすごろくについて来るようなサイコロのような場合である「かなり精度のよいサイコロであると信じるに足る」ということはどういうことかというと『実験してみるまでもなく』dice<-c(100,100,100,100,100,100)と言っていいんじゃないかと思わせられる、ということ意訳すると、ためしに100回振ってみたら、1が出る確率がhist(rdirichlet(100,c(100,100,100,100,100,100)[,1]))で得られるくらいには正確なサイコロ、だろうと信じている、ということこのようなサイコロだけれども、やはり、本当はどれくらい「よくできたサイコロ」かを知りたいから、実験してみるとする100回投げたら、(18,20,16,12,18,16)と出た。このときは、rdirichlet(k,c(18,20,16,12,18,16))で101回目からの目の出る確率を予想してよさそうに思われるでは、10回だけ投げて(0,4,2,1,2,1)と出たとしよう。このときはどうだろうかそもそもrdirichlet(k,c(0,4,2,1,2,1))はエラーメッセージが出る。1の目が出た回数が0(1の目の出る事前情報が確率0)だからである。Rの関数的には、rdirichlet関数が呼び出しているrgamma分布の関数にこの0が渡されるのだが、rgamma関数の引数shapeが0で定義されていないからである。genotypeデータからのハプロタイプ推定にもあり、それについてのコメントは、アプリケーションsnphapの解説文にもある。GWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク複合遺伝性疾患とは----------News----------20万アクセス2008/07/2040万アクセス2010/09/12-------------------------はてなブログをはじめよう!ryamada22さんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。