{"id":565,"date":"2018-06-24T11:37:27","date_gmt":"2018-06-24T02:37:27","guid":{"rendered":"https:\/\/best-biostatistics.com\/toukei-er\/entry\/gibbs-sampler-by-r\/"},"modified":"2024-10-14T11:55:33","modified_gmt":"2024-10-14T02:55:33","slug":"gibbs-sampler-by-r","status":"publish","type":"post","link":"https:\/\/best-biostatistics.com\/toukei-er\/entry\/gibbs-sampler-by-r\/","title":{"rendered":"R \u3067\u30ae\u30d6\u30b9\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3092\u5b9f\u65bd\u3059\u308b\u65b9\u6cd5"},"content":{"rendered":"\n

\u30ae\u30d6\u30b9\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\uff08\u30ae\u30d6\u30b9\u6cd5\uff09\u3092 R \u3067\u5b9f\u65bd\u3059\u308b\u65b9\u6cd5\u3002<\/p>\n\n\n\n

\u521d\u671f\u5024\u3092\u4e0e\u3048\u3066\u3001\u4f55\u5ea6\u3082\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3057\u3066\u3044\u308b\u3068\u3001\u3060\u3093\u3060\u3093\u306b\u5f62\u306b\u306a\u3063\u3066\u304f\u308b\u3002<\/p>\n\n\n\n\n\n\n\n

\u30ae\u30d6\u30b9\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u7528\u306e\u30c7\u30fc\u30bf\u306e\u7528\u610f<\/h2>\n\n\n\n

heikin\u3068\u3044\u3046\u540d\u524d\u306e\u30c7\u30fc\u30bf\u3092\u5165\u529b\u3059\u308b\u3002<\/p>\n\n\n\n

heikin <-<\/span> c<\/span>(<\/span>6.0<\/span>,<\/span>10.0<\/span>,<\/span>7.6<\/span>,<\/span>3.5<\/span>,<\/span>1.4<\/span>,<\/span>2.5<\/span>,<\/span>5.6<\/span>,<\/span>3.0<\/span>,<\/span>2.2<\/span>,<\/span>5.0<\/span>,<\/span>\n3.3<\/span>,<\/span>7.6<\/span>,<\/span>5.8<\/span>,<\/span>6.7<\/span>,<\/span>2.8<\/span>,<\/span>4.8<\/span>,<\/span>6.3<\/span>,<\/span>5.3<\/span>,<\/span>5.4<\/span>,<\/span>3.3<\/span>,<\/span>\n3.4<\/span>,<\/span>3.8<\/span>,<\/span>3.3<\/span>,<\/span>5.7<\/span>,<\/span>6.3<\/span>,<\/span>8.4<\/span>,<\/span>4.6<\/span>,<\/span>2.8<\/span>,<\/span>7.9<\/span>,<\/span>8.9<\/span>)<\/span>\nmean1 <-<\/span> mean<\/span>(<\/span>heikin)<\/span>\nn <-<\/span> length<\/span>(<\/span>heikin)<\/span>\nmean<\/span>(<\/span>heikin)<\/span>\nvar<\/span>(<\/span>heikin)<\/span>\n<\/code><\/pre>\n\n\n\n
heikin\u30c7\u30fc\u30bf\u306e\u5e73\u5747\u3001\u5206\u6563\u3001\u30b5\u30f3\u30d7\u30eb\u30b5\u30a4\u30ba\u306f\u3001\u4ee5\u4e0b\u306e\u901a\u308a\u3002<\/p>\n\n\n\n
><\/span> mean<\/span>(<\/span>heikin)<\/span>\n[<\/span>1<\/span>]<\/span> 5.106667<\/span>\n><\/span> var<\/span>(<\/span>heikin)<\/span>\n[<\/span>1<\/span>]<\/span> 4.766161<\/span>\n><\/span> length<\/span>(<\/span>heikin)<\/span>\n[<\/span>1<\/span>]<\/span> 30<\/span>\n<\/code><\/pre>\n\n\n\n
\u30ae\u30d6\u30b9\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3092\u884c\u3046\u305f\u3081\u306e\u30d1\u30c3\u30b1\u30fc\u30b8\u3092\u30a4\u30f3\u30b9\u30c8\u30fc\u30eb\u3059\u308b<\/h2>\n\n\n\n
pscl\u30d1\u30c3\u30b1\u30fc\u30b8\u3092\u30a4\u30f3\u30b9\u30c8\u30fc\u30eb\u3059\u308b\u3002\u4e00\u56de\u3060\u3051\u3002<\/p>\n\n\n\n
install.packages<\/span>(<\/span>\"pscl\"<\/span>)<\/span>\n<\/code><\/pre>\n\n\n\n
pscl\u30d1\u30c3\u30b1\u30fc\u30b8\u3092\u3001\u547c\u3073\u51fa\u3059\u3002<\/p>\n\n\n\n
library<\/span>(<\/span>pscl)<\/span>\n<\/code><\/pre>\n\n\n\n
\uff1e\uff1e\u3082\u3046\u7d71\u8a08\u3067\u60a9\u3080\u306e\u306f\u7d42\u308f\u308a\u306b\u3057\u307e\u305b\u3093\u304b\uff1f\u00a0<\/a><\/strong><\/span><\/p>\r\n\"\"<\/a>\r\n
\u21911\u4e07\u4eba\u4ee5\u4e0a\u306e\u533b\u7642\u5f93\u4e8b\u8005\u304c\u8cfc\u8aad\u4e2d<\/span><\/strong><\/span><\/p><\/div>
\u30ae\u30d6\u30b9\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u306e\u4e8b\u524d\u5206\u5e03\u306e\u56f3\u793a<\/h2>\n\n\n\n
0.05\u304b\u308910\u307e\u3067\u30010.05\u523b\u307f\u306e\u5909\u6570x\u3092\u4f5c\u308b\u3002<\/p>\n\n\n\n
alpha=0.01, beta=0.01\u306e\u9006\u30ac\u30f3\u30de\u5206\u5e03\u3092\u63cf\u304f\u3002<\/p>\n\n\n\n
\u3053\u308c\u304c\u4e8b\u524d\u5206\u5e03\u3002<\/p>\n\n\n\n
x <-<\/span> seq<\/span>(<\/span>0.05<\/span>,<\/span>10<\/span>,<\/span>0.05<\/span>)<\/span>\ncurve<\/span>(<\/span>densigamma<\/span>(<\/span>x,<\/span>alpha=<\/span>0.01<\/span>,<\/span> beta=<\/span>0.01<\/span>),<\/span>from=<\/span>0.05<\/span>,<\/span>to=<\/span>10<\/span>,<\/span>lwd=<\/span>3<\/span>)<\/span>\n<\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/span><\/p>\n\n\n\n
\u30ae\u30d6\u30b9\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u306e\u305f\u3081\u306e\u6761\u4ef6\u8a2d\u5b9a<\/h2>\n\n\n\n
\u3053\u3053\u3067\u6761\u4ef6\u8a2d\u5b9a\u3002<\/p>\n\n\n\n
n0 <-<\/span> 2<\/span>*<\/span>0.01<\/span>\nn0S0 <-<\/span> 2<\/span>*<\/span>0.01<\/span>\nS0 <-<\/span> n0S0\/<\/span>n0\nmu0 <-<\/span> 5<\/span>\nm0 <-<\/span> 1<\/span>\/<\/span>4<\/span>\nm1 <-<\/span> m0+<\/span>n\nn1 <-<\/span> n0+<\/span>n\nmu1 <-<\/span> (<\/span>n*<\/span>mean1+<\/span>m0*<\/span>mu0)<\/span>\/<\/span>(<\/span>m0+<\/span>n)<\/span>\nQ <-<\/span> sum<\/span>((<\/span>heikin-<\/span>mean1)<\/span>^<\/span>2<\/span>)<\/span>\nn1S1 <-<\/span> n0S0+<\/span>Q+<\/span>m0*<\/span>n\/<\/span>(<\/span>m0+<\/span>n)<\/span>*<\/span>(<\/span>mean1-<\/span>mu0)<\/span>^<\/span>2<\/span>\n<\/code><\/pre>\n\n\n\n
\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u306e\u305f\u3081\u306e\u95a2\u6570\u3092\u6e96\u5099\u3059\u308b<\/h2>\n\n\n\n
\u5206\u6563\u306e\u4e8b\u5f8c\u5206\u5e03\u304b\u3089\u6bcd\u5206\u6563 $ \\sigma^2 $ \u3092\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3059\u308b\u95a2\u6570\u3002<\/p>\n\n\n\n
my.sigma2 <-<\/span> function<\/span>(<\/span>mu){<\/span>\nsigma2 <-<\/span> rigamma<\/span>(<\/span>1<\/span>,<\/span> alpha=<\/span>(<\/span>n1+<\/span>1<\/span>)<\/span>\/<\/span>2<\/span>,<\/span> beta=<\/span>(<\/span>n1S1+<\/span>m1*<\/span>(<\/span>mu-<\/span>mu1)<\/span>^<\/span>2<\/span>)<\/span>\/<\/span>2<\/span>)<\/span>\nreturn<\/span>(<\/span>sigma2)<\/span>\n}<\/span>\n<\/code><\/pre>\n\n\n\n
\u5e73\u5747\u306e\u4e8b\u5f8c\u5206\u5e03\u304b\u3089\u6bcd\u5e73\u5747\u3092\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3059\u308b\u95a2\u6570\u3002<\/p>\n\n\n\n
my.mu <-<\/span> function<\/span>(<\/span>sigma2){<\/span>\nmu <-<\/span> rnorm<\/span>(<\/span>1<\/span>,<\/span> mean=<\/span>mu1,<\/span> sd=<\/span>sqrt<\/span>(<\/span>sigma2\/<\/span>m1))<\/span>\nreturn<\/span>(<\/span>mu)<\/span>\n}<\/span>\n<\/code><\/pre>\n\n\n\n
\u30e2\u30f3\u30c6\u30ab\u30eb\u30ed\u6cd5\u306b\u3088\u308b\u30b7\u30df\u30e5\u30ec\u30fc\u30b7\u30e7\u30f3<\/h2>\n\n\n\n
\u30e2\u30f3\u30c6\u30ab\u30eb\u30ed \u30b7\u30df\u30e5\u30ec\u30fc\u30b7\u30e7\u30f3\u3092\u5b9f\u884c\u3059\u308b\u3002<\/p>\n\n\n\n
my.montecarlo <-<\/span> function<\/span>(<\/span>n){<\/span>\nset.seed<\/span>(<\/span>20110608<\/span>)<\/span>\nmu <-<\/span> 5<\/span>\nnew.mu <-<\/span> 5<\/span>\nsigma2 <-<\/span> 4<\/span>\nfor<\/span> (<\/span>i in<\/span> 1<\/span>:<\/span>n){<\/span>\nnew.sigma2 <-<\/span> my.sigma2<\/span>(<\/span>new.mu)<\/span>\nsigma2 <-<\/span> sigma2 +<\/span> new.sigma2\nnew.mu <-<\/span> my.mu<\/span>(<\/span>new.sigma2)<\/span>\nmu <-<\/span> mu +<\/span> new.mu\n}<\/span>\nlist<\/span>(<\/span>\"mu\"<\/span>=<\/span>mu\/<\/span>n,<\/span> \"sigma2\"<\/span>=<\/span>sigma2\/<\/span>n)<\/span>\n}<\/span>\n<\/code><\/pre>\n\n\n\n
10\u4e07\u56de\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3057\u305f\u7d50\u679c\u306f\u4ee5\u4e0b\u306e\u901a\u308a\u3002<\/p>\n\n\n\n
><\/span> my.montecarlo<\/span>(<\/span>100000<\/span>)<\/span>\n$<\/span>`mu`\n[<\/span>1<\/span>]<\/span> 5.102863<\/span>\n$<\/span>sigma2\n[<\/span>1<\/span>]<\/span> 4.936363<\/span>\n<\/code><\/pre>\n\n\n\n
\u30ae\u30d6\u30b9\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3067\u5f97\u3089\u308c\u305f\u4e8b\u5f8c\u5206\u5e03\u306e\u56f3\u793a<\/h2>\n\n\n\n
\u5206\u6563\u306e\u4e8b\u5f8c\u5206\u5e03\u3092\u63cf\u3044\u3066\u307f\u308b\u3002<\/p>\n\n\n\n
curve<\/span>(<\/span>densigamma<\/span>(<\/span>x,<\/span> alpha=<\/span>(<\/span>n1+<\/span>1<\/span>)<\/span>\/<\/span>2<\/span>,<\/span> beta=<\/span>(<\/span>n1S1+<\/span>m1*<\/span>(<\/span>5.102863<\/span>-<\/span>mu1)<\/span>^<\/span>2<\/span>)<\/span>\/<\/span>2<\/span>),<\/span>from=<\/span>0.01<\/span>,<\/span>to=<\/span>8.0<\/span>,<\/span> ylim=<\/span>c<\/span>(<\/span>0<\/span>,<\/span>0.4<\/span>),<\/span>lwd=<\/span>3<\/span>)<\/span>\n<\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/span><\/p>\n\n\n\n
\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u7d50\u679c\u3092\u4f7f\u3063\u3066\u5e73\u5747\u306e\u5206\u5e03\u3092\u63cf\u3044\u3066\u307f\u308b\u3002<\/p>\n\n\n\n
curve<\/span>(<\/span>dnorm<\/span>(<\/span>x,<\/span> mean=<\/span>mu1,<\/span> sd=<\/span>sqrt<\/span>(<\/span>4.936363<\/span>\/<\/span>m1)),<\/span>from=<\/span>0<\/span>,<\/span>to=<\/span>10<\/span>,<\/span>lwd=<\/span>3<\/span>)<\/span>\n<\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/span><\/p>\n\n\n\n
\u307e\u3068\u3081<\/h2>\n\n\n\n
\u30ae\u30d6\u30b9\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3092\u53c2\u8003\u66f8\u7c4d\u901a\u308a\u306b R \u3067\u884c\u3063\u3066\u307f\u305f\u3002<\/p>\n\n\n\n
\u53c2\u8003\u66f8\u7c4d<\/h2>\n\n\n\n
\"\u9053\u5177\u3068\u3057\u3066\u306e\u30d9\u30a4\u30ba\u7d71\u8a08\"<\/a><\/figure>\n\n\n\n
<\/div>\n\n\n\n
\n
\u9053\u5177\u3068\u3057\u3066\u306e\u30d9\u30a4\u30ba\u7d71\u8a08<\/a><\/p>\n