全ての確率はコイン投げに通ず

શͯͷ֬཰͸ ίΠϯ౤͛ʹ௨ͣ +BQBO3ൃදࢿྉ ,FO`JDIJ.BUTVJ !LFONBUTV

ࣗݾ঺հ!LFONBUTV ɾ'BDFCPPLϖʔδ ɹIUUQTXXXGBDFCPPLDPNNBUTVLFOCPPL ɾ5XJUUFSΞΧ΢ϯτ ɹ!LFONBUTV ɾ2JJUBͰϒϩάΛॻ͍͍ͯ·͢ʢ౷ܭɺػցֶशɺ1ZUIPO౳ʣ ɹɹɹIUUQRJJUBDPNLFONBUTV ɹɹɹʢDPOUSJCVUJPOΛ௒͑·ͨ͠ʂʣ ɾझຯ ɹɹɹόϯυͰϕʔεΛ஄͍ͨΓ͍ͯ͠·͢ɻ
ɹɹɹओʹ౦ೆΞδΞ΁όοΫύοΫཱྀߦʹߦͬͨΓ͠·͢ ɹɹʢΧϯϘδΞɺϛϟϯϚʔɺόϯάϥσγϡɺ৽ᙜ΢Πάϧ࣏ࣗ۠FUDʣ ཱྀߦͷࣸਅIUUQNBUTVLFOKJNEPDPN TwitterΞΠίϯ

1ZUIPOλάͰ̍Ґʹ ͳΓ·ͨ͠ʂ 㱼ਓ

http://www.slideshare.net/matsukenbook

ຊ୊Ͱ͢ɻ

http://nbviewer.ipython.org/github/bakfoo/pyconjp/blob/master/probprog3.ipynb ݁ہɼ֬཰తϓϩάϥϛϯάʹ͓͍ͯɼϓϩάϥϚ͕΍Δ͜ͱ͸ɼ ɹw σʔλΛ༻ҙ͢Δ ɹw ύϥϝʔλͷछྨΛܾΊɼ֬཰෼෍ΛܾΊΔ ɹw
ύϥϝʔλࠐΈͰࣄલ෼෍Λ਺ࣜͰද͢ ɹw ύϥϝʔλࠐΈͰ໬౓Λ਺ࣜͰද͢ ɹw .$.$αϯϓϦϯά๏ΛબͿ ͋ͱ͸ܭࢉػ͕.$.$αϯϓϦϯάΛͯ͠ɼύϥϝʔλͷਪఆ஋Λฦ͢ɽ ݁ہɼࣄલ෼෍΍໬౓ͷ֬཰෼෍Λ্खʹϞσϦϯά͢Δʹ͸ɼ ֬཰෼෍ͱͦͷ֬཰෼෍ͷؔ܎ΛΑ͘஌͍ͬͯͳ͚Ε͹ͳΒͳ͍ɽ ദ໺͞Μ!1ZDPO+1 ͳΔ΄Ͳʂ

ͦ͜Ͱ঺հ͞Ε͍ͯͨͷ͕ɺ

͜Εʂ http://www.math.wm.edu/~leemis/2008amstat.pdf

͜Εʂ http://www.math.wm.edu/~leemis/2008amstat.pdf ଟ͍ͳʂ ʀƅЧƅ

͜Εʂ http://www.math.wm.edu/~leemis/2008amstat.pdf ଟ͍ͳʂ ʀƅЧƅ ݸ͋Γ·͢

෼ؒͰग़དྷΔݶΓ આ໌ͯ͠Έ·͢ʂ ʆɾТɾ

Let’s Start!!!

ϕϧψʔΠ෼෍ #FSOPVMMJ%JTUSJCVUJPO

ϕϧψʔΠ෼෍ ද͸ ཪ͸ ࢼߦΛ̍౓ߦ͍ɺ੒ޭ͔ࣦഊͳͲ̎஋ ྫɿίΠϯ౤͛ ֬཰ؔ਺ P ( X =
x ) = px(1 p )1 x ( x = 0 , 1) ύϥϝʔλʔ Qɿ ද ͕ग़Δ֬཰

ϕϧψʔΠ෼෍ # ϕϧψʔΠ෼෍͔ΒͷαϯϓϦϯάΛ࣮ߦ # ύϥϝʔλʔ p = 0.7
trial_size = 10000 set.seed(71) # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒ data <-‐ rbern(trial_size, p) # ϕϧψʔΠ෼෍ͷ֬཰෼෍ dens <-‐ data.frame(y=c((1-‐p),p)*trial_size, x=c(0, 1)) # άϥϑඳը ggplot() + layer(data=data.frame(x=data), mapping=aes(x=x), geom="bar", stat="bin", bandwidth=0.1 ) + layer(data=dens, mapping=aes(x=x, y=y), geom="bar", stat="identity", width=0.05, fill="#777799", alpha=0.7) 3ίʔυ ࢵɿཧ࿦తͳ֬཰ ࠇɿཚ਺͔Βੜ੒ͨ͠ώετάϥϜ

͜͜Ͱɺ ΋͏Ұ౓λΠτϧΛ

શͯͷ֬཰͸ ίΠϯ౤͛ʹ௨ͣ +BQBO3ൃදࢿྉ ,FO`JDIJ.BUTVJ !LFONBUTV

ೋ߲෼෍ #JOPNJBM%JTUSJCVUJPO

ೋ߲෼෍ ϕϧψʔΠࢼߦΛOճ࣮ࢪ͠ɺԿճ੒ޭ͔ͨ͠ͷճ਺ Λ֬཰ม਺ͱͨ͠෼෍ɻ ྫɿϑϦʔεϩʔ̍ηοτճ࣮ࢪɻ ɹɹ̍ճ͋ͨΓͷ੒ޭ཰ͰܭԿճೖ͔ͬͨ ֬཰ؔ਺ ύϥϝʔλʔ Qɿ ද ͕ग़Δ֬཰<ϑϦʔεϩʔͷ੒ޭ֬཰>
P ( X = x ) = nCrpx(1 p )n x Oɿηοτ͋ͨΓͷճ਺ ( x = 1 , 2 , · · · , n )

ೋ߲෼෍ Oճ࣮ࢪ දͷग़ͨճ਺Y ίΠϯ౤͛ දͷ֬཰Q

ೋ߲෼෍ # ύϥϝʔλʔ p = 0.7 trial_size =
10000 sample_size = 30 set.seed(71) # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒ gen_binom_var <-‐ function() { return(sum(rbern(sample_size, p))) } result <-‐ rdply(trial_size, gen_binom_var()) # ೋ߲෼෍ͷີ౓ؔ਺ dens <-‐ data.frame(y=dbinom(seq(sample_size), sample_size, 0.7))*trial_size # άϥϑඳը ggplot() + layer(data=resuylt, mapping=aes(x=V1), geom="bar", stat = "bin", binwidth=1, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=seq(sample_size)+.5, y=y), geom="line", stat="identity", position="identity",colour="red" ) + ggtitle("Bernoulli to Binomial.") ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒ 3ίʔυ

ϙϫιϯ෼෍ 1PJTTPO%JTUSJCVUJPO

ϙϫιϯ෼෍ ηοτճ਺ O ˠ 㱣 ɺ̍ճ͋ͨΓ੒ޭ཰Q ˠ
OQΛҰఆͱͨ࣌͠ͷ෼෍ ྫɿ๮ͷইͷνΣοΫΛߟ͑ΔɻDN෯Ͱইͷ༗ແΛ Λௐ΂ΔɻϕϧψʔΠࢼߦͷ܁Γฦ͠ ˠ ೋ߲෼෍ ˠແݶʹ෼ׂͯ͠ϙΞιϯ෼෍ɻ DN DN DN DN DN DN DN DN DN DN DN શ෦Ͱই͕̓ͭ DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN 㱣 ʹ෼ׂ

֬཰ؔ਺ ύϥϝʔλʔ ͜ͷྫͷ৔߹ɺ๮Ұຊ͋ͨΓͷইͷฏۉతͳ਺ Еɿฏۉੜىճ਺ P ( X = x )
= e x x ! ( x = 1 , 2 , · · · , 1) ЕOQͷؔ܎͕อͨΕ͍ͯΔ ๮̍ຊ͋ͨΓͷইͷ਺ ͷ෼෍ˠϙΞιϯ෼෍ ্ݶ͕ͳ͘ͳͬͨʂ ϙΞιϯ෼෍

Oճ࣮ࢪ ϙΞιϯ෼෍ ίΠϯ౤͛ දͷ֬཰Q OQҰఆͰOˠ㱣

ϙΞιϯ෼෍ trial_size = 5000; width <-‐ 1; # ΋ͱ΋ͱͷ໰୊ઃఆ
p = 0.7; n = 10; np <-‐ p*n # n→∞ɺp→0ɺnp=Ұఆ n = 100000; p <-‐ np/n # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒ gen_binom_var <-‐ function() { return(sum(rbern(n, p))) } result <-‐ rdply(trial_size, gen_binom_var()) # ϙΞιϯ෼෍ͷີ౓ؔ਺ dens <-‐ data.frame(y=dpois(seq(20), np))*trial_size # άϥϑඳը ggplot() + layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=seq(20)+.5, y=y), geom="line", stat="identity", position="identity", colour="red" ) + ggtitle("Bernoulli to Poisson.") ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒ 3ίʔυ

ਖ਼ن෼෍ /PSNBM%JTUSJCVUJPO

ਖ਼ن෼෍ ೋ߲෼෍ͷOΛେ͖͘͢Δͱਖ਼ن෼෍ͰۙࣅͰ͖Δɻ ֬཰ີ౓ؔ਺ ύϥϝʔλʔ Мɿඪ४తภࠩ Жɿฏۉ஋ f(x) = 1 p
2⇡ 2 exp ⇢ 1 2 (x µ) 2 2 ( 1 < x < 1) ྫɿࢮ͵·ͰϑϦʔεϩʔ౤͛ͯܭԿճೖ͔ͬͨ

ਖ਼ن෼෍ ϙΞιϯ෼෍ͱͷҧ͍ DN DN DN DN DN DN DN DN
DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN DN ʜ DN DN DN DN DN DN DN DN DN DN ˠ㱣DN DN DN DN DN DN DN DN DN DN Ԇʑͱ৳͹͢

Oճ࣮ࢪ ίΠϯ౤͛ දͷ֬཰Q Oˠ㱣 ਖ਼ن෼෍

ਖ਼ن෼෍ # ύϥϝʔλʔ n <-‐ 10000; p <-‐ 0.7;
trial_size = 10000 width=10 # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒ gen_binom_var <-‐ function() { return(sum(rbern(n, p))) } result <-‐ rdply(trial_size, gen_binom_var()) # ਖ਼ن෼෍ͷີ౓ؔ਺ dens <-‐ data.frame(y=dnorm(seq(6800,7200), mean=n*p, sd=sqrt(n*p*(1-‐p)))*trial_size*width) # άϥϑඳը ggplot() + layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=seq(6800,7200), y=y), geom="line", stat="identity", position="identity", colour="red") + ggtitle("Bernoulli to Normal.") 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒

ඪ४ਖ਼ن෼෍ 4UBOEBSE/PSNBM%JTUSJCVUJPO

ඪ४ਖ਼ن෼෍ ਖ਼ن෼෍Λฏۉ̌ɺඪ४ภࠩ̍ʹͨ͠΋ͷ ֬཰ີ౓ؔ਺ ύϥϝʔλʔ ͳ͠ ( 1 < x <
1) f(x) = 1 p 2⇡ exp ⇢ 1 2 x 2

Oճ࣮ࢪ ίΠϯ౤͛ දͷ֬཰Q Oˠ㱣 ඪ४ਖ਼ن෼෍ ฏۉΛʹͣΒ͢ ඪ४ภࠩΛʹ ॖΊΔ

ඪ४ਖ਼ن෼෍ # ύϥϝʔλʔ n <-‐ 10000; p <-‐ 0.7
trial_size = 30000 width=0.18 # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒ gen_binom_var <-‐ function() { return(sum(rbern(n, p))) } result <-‐ rdply(trial_size, gen_binom_var()) m <-‐ mean(result$V1); sd <-‐ sd(result$V1); result <-‐ (result -‐ m)/sd # ඪ४ਖ਼ن෼෍ͷີ౓ؔ਺ dens <-‐ data.frame(y=dnorm(seq(-‐4,4,0.05), mean=0, sd=1)*trial_size*width) # άϥϑඳը ggplot() + layer(data=result, mapping=aes(x=V1), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=seq(-‐4,4,0.05), y=y), geom="line", stat="identity", position=“identity", colour="red" ) + ggtitle("Bernoulli to Standard Normal.") 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒

ΧΠೋ৐෼෍ $IJ4RVBSF%JTUSJCVUJPO

ΧΠೋ৐෼෍ ඪ४ਖ਼ن෼෍ʹै͏֬཰ม਺Λ̎৐ͯ͠Lݸ ଍͠߹Θͤͨ ɹɹɹɹɹɹɹɹɹɹɹɹɹɹɹɹɹ͕ै͏෼෍ ֬཰ີ౓ؔ਺ ύϥϝʔλʔ Lɿࣗ༝౓ ̎৐ͨ͠ඪ४ਖ਼ن෼෍ͷ֬཰ม਺ͷ ɹɹɹɹɹ଍͠߹Θͤͨ਺ f
( x, k ) = (1 / 2)k/ 2 ( k/ 2) xk/ 2 1 e x/ 2 (0  x ) Xi Z = X2 1 + · · · + X2 k

ฏۉΛʹͣΒ͢ Oճ࣮ࢪ ίΠϯ౤͛ දͷ֬཰Q Oˠ㱣 ΧΠೋ৐෼෍ ฏۉΛʹͣΒ͢ ඪ४ภࠩΛʹ ॖΊΔ ฏۉΛʹͣΒ͢
Lͷͱ͖ ೋ৐ͨ͠ͱ͜Ζ ʹੵΜͰ͍͘ Lͷͱ͖ Lͭ αϯϓϦϯά͢Δ Lͭͷ௕͞Λ ଍ͨ͠΋ͷΛ ϓϩοτ͢Δ

ΧΠೋ৐෼෍ # ύϥϝʔλʔ p <-‐ 0.7; n <-‐ 1000;
trial_size <-‐ 100000; width <-‐ 0.3; df <-‐ 3 # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒(3·Θ͠) gen_binom_var <-‐ function() { return(sum(rbern(n, p))) } gen_chisq_var <-‐ function() { result <-‐ rdply(trial_size, gen_binom_var()) return(((result$V1 -‐ mean(result$V1))/sd(result$V1))**2) } # ࣗ༝౓dfͷ෼͚ͩੜ੒͢Δ result <-‐ rlply(df, gen_chisq_var(),.progress = "text") res <-‐ data.frame(x=result[[1]] + result[[2]] + result[[3]]) # ΧΠೋ৐෼෍ͷີ౓ؔ਺(ࣗ༝౓=3) xx <-‐ seq(0,20,0.1) dens <-‐ data.frame(y=dchisq(x=xx, df=df)*trial_size*width) # άϥϑඳը ggplot() + layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=xx, y=y), geom="line", stat="identity", position="identity", colour="blue" ) + ggtitle("Bernoulli to Chisquare") 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒

ࢦ਺෼෍ &YQPOFOUJBM%JTUSJCVUJPO

ࢦ਺෼෍ ϙΞιϯ෼෍͕୯Ґ௕ʢ΋͘͠͸୯Ґ࣌ؒʣʹฏۉ Կճੜى͢Δ͔ɺΛද͢෼෍Ͱ͕͋ͬͨɺಉ͡ࣄ৅ Λҧ͏֯౓Ͱଊ͑௚ͨ͠ͷ͕ࢦ਺෼෍ɻ DN DN DN DN DN DN
DN DN DN DN DN શ෦Ͱই͕̓ͭ ϙΞιϯ෼෍తߟ͑ํ ࢦ਺෼෍తߟ͑ํ ࣍ൃੜ͢Δ·ͰʹɺͲͷ͘Β͍ͷڑ཭ ࣌ؒ ͔

ࢦ਺෼෍ ֬཰ີ౓ؔ਺ ύϥϝʔλʔ Еɿ୯Ґ௕ʢ࣌ؒʣʹฏۉతʹੜى͢Δճ਺ f ( x, ) = ⇢
e x ( x 0) 0 ( x < 0)

Oճ࣮ࢪ ίΠϯ౤͛ දͷ֬཰Q OQҰఆͰOˠ㱣 ࢦ਺෼෍ ฏۉΛʹͣΒ͢ ִؒͷ௕͞Λϓϩοτ

trial_size = 7000; width <-‐ .01; # ΋ͱ΋ͱͷ໰୊ઃఆ
p = 0.7; n = 10; np <-‐ p*n; # n→∞ɺp→0ɺnp=Ұఆ n = 10000; p <-‐ np/n # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒ gen_exp_var <-‐ function() { cnt <-‐ 0 while (TRUE) { cnt <-‐ cnt + 1 if (rbern(1, p)==1){ return(cnt) # 1͕ग़ͨΒͦΕ͕Կճ໨͔Λฦ͢ } } } data <-‐ data.frame(x=rdply(trial_size, gen_exp_var())/n) names(data) <-‐ c("n", "x") # ࢦ਺෼෍ͷີ౓ؔ਺ dens <-‐ data.frame(y=dexp(seq(0, 1.5, 0.1), np)*trial_size*width) ggplot() + layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=seq(0, 1.5, 0.1), y=y), geom="line", stat="identity", position="identity", colour="red" ) + ggtitle("Bernoulli to Exponential.") 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒ ࢦ਺෼෍

ΨϯϚ෼෍ (BNNB%JTUSJCVUJPO

ΨϯϚ෼෍ ֬཰ີ౓ؔ਺ ύϥϝʔλʔ ୯Ґ௕ʢ࣌ؒʣ͋ͨΓЌճൃੜ͢Δࢦ਺෼෍ͷ Ћݸͷ࿨ͷ֬཰෼෍Λߟ͑ͨ࣌ɺΨϯϚ෼෍ͱͳΔ f(x, ↵, ) = ↵
(↵) x ↵ 1 exp( x) (0  x < 1) ↵ X i=1 Xi ⇠ (↵, ) Xi ⇠ Exp( )

Oճ࣮ࢪ ίΠϯ౤͛ දͷ֬཰Q OQҰఆͰOˠ㱣 ฏۉΛʹͣΒ͢ ΨϯϚ෼෍ ִؒͷ௕͞Λϓϩοτ Lͭͷ௕͞Λ ଍ͨ͠΋ͷΛ ϓϩοτ͢Δ
LͭαϯϓϦϯά

trial_size = 7000; width <-‐ .035; # ΋ͱ΋ͱͷ໰୊ઃఆ
p = 0.7; n = 10; np <-‐ p*n; # n→∞ɺp→0ɺnp=Ұఆ n = 10000; p <-‐ np/n; alpha <-‐ 5 # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒ get_interval <-‐ function(){ cnt <-‐ 0 while (TRUE) { cnt <-‐ cnt + 1 if (rbern(1, p)==1){ return(cnt) } } } gen_exp_var <-‐ function() { data <-‐ data.frame(x=rdply(trial_size, get_interval())/n) names(data) <-‐ c("n", "x") return(data) } result <-‐ rlply(alpha, gen_exp_var()) data <-‐ data.frame(x=result[[1]]$x + result[[2]]$x + result[[3]]$x + result[[4]]$x + result[[5]]$x) # ΨϯϚ෼෍ͷີ౓ؔ਺ dens <-‐ data.frame(y=dgamma(seq(0, 3,.01), shape=alpha, rate=np)*trial_size*width) ggplot() + layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=seq(0,3,.01), y=y), geom="line", stat="identity", position="identity", colour="red" ) + ggtitle("Bernoulli to Gamma") 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒ ΨϯϚ෼෍

ٯΨϯϚ෼෍ *OWFSTF(BNNB%JTUSJCVUJPO

ٯΨϯϚ෼෍ ֬཰ີ౓ؔ਺ ύϥϝʔλʔ f(x, ↵, ) = ↵ (↵) x
(↵+1) exp ✓ x ◆ (0  x < 1) ୯Ґ௕ʢ࣌ؒʣ͋ͨΓЌճൃੜ͢Δࢦ਺෼෍ͷЋݸͷ࿨ͷ ֬཰෼෍Λߟ͑ͨ࣌ɺΨϯϚ෼෍ͱͳΔ Xi ⇠ Exp( ) Z = ↵ X i=1 Xi ⇠ (↵, ) 1/Z ⇠ IG(↵, ) ͜ͷ࣌ɺ;ͷٯ਺͕ै͏෼෍ΛٯΨϯϚ෼෍ͱݴ͏ɻ

ฏۉΛʹͣΒ͢ Oճ࣮ࢪ ίΠϯ౤͛ දͷ֬཰Q OQҰఆͰOˠ㱣 ฏۉΛʹͣΒ͢ ִؒͷ௕͞Λϓϩοτ ٯΨϯϚ෼෍ Lͭͷ௕͞Λ ଍ͨ͠΋ͷΛ
ϓϩοτ͢Δ Lͭ αϯϓϦϯά શͯͷσʔλΛ ZYͰม׵ [ [ X X

trial_size = 7000; width <-‐ .; # ΋ͱ΋ͱͷ໰୊ઃఆ
p = 0.7; n = 10; np <-‐ p*n; # n→∞ɺp→0ɺnp=Ұఆ n = 10000; p <-‐ np/n; alpha <-‐ 5 # ϕϧψʔΠ෼෍ʹै͏ཚ਺ੜ੒ get_interval <-‐ function(){ cnt <-‐ 0 while (TRUE) { cnt <-‐ cnt + 1 if (rbern(1, p)==1){ return(cnt) } } } gen_exp_var <-‐ function() { data <-‐ data.frame(x=rdply(trial_size, get_interval())/n) names(data) <-‐ c("n", "x") return(data) } result <-‐ rlply(alpha, gen_exp_var()) data <-‐ data.frame(x=1/(result[[1]]$x + result[[2]]$x + result[[3]]$x + result[[4]]$x + result[[5]]$x)) # ٯΨϯϚ෼෍ͷີ౓ؔ਺ dens <-‐ data.frame(y=dinvgamma(seq(0, 23,.01), shape=5, rate=1/np)*trial_size*width) ggplot() + layer(data=data, mapping=aes(x=x), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=seq(0,3,.01), y=y), geom="line", stat="identity", position="identity", colour="red" ) + ggtitle("Bernoulli to Inversegamma") 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒ ٯΨϯϚ෼෍

ඪ४Ұ༷෼෍ 4UBOEBSE6OJGPSN%JTUSJCVUJPO

ඪ४Ұ༷෼෍ ֬཰ີ౓ؔ਺ ύϥϝʔλʔ ͳ͠ ͔ΒͷؒͰ౳֬཰Ͱൃੜ͢ΔΑ͏ͳݱ৅Λ ද͢෼෍ f(x) = ⇢ 1
(0  x  1) 0 (otherwise)

ඪ४Ұ༷෼෍ ίΠϯ౤͛ දͷ֬཰Q Rճ࣮ࢪ ࢼߦ݁ՌΛ֤ܻʹׂΓ౰ͯ ද ཪ ද ཪ ද
ཪ ճ໨ ճ໨ Rճ໨ ˞΋ͬͱޮ཰ͷྑ͍΍Γํ΋͋Δͱࢥ͍·͕͢෼͔Γ΍͢͞ͷͨΊʜ Z = x1(1 / 2)1 + x2(1 / 2)2 + · · · + xq(1 / 2)q

width <-‐ 0.02 p <-‐ 0.5; sample_size <-‐
1000 trial_size <-‐ 100000 gen_unif_rand <-‐ function() { # sample_sizeܻͷ2ਐগ਺ΛϕϧψʔΠ෼෍ʹ # ै͏ཚ਺͔Βੜ੒ return (sum(rbern(sample_size, p) * (rep(1/2, sample_size) ** seq(sample_size)))) } gen_rand <-‐ function(){ return( rdply(trial_size, gen_unif_rand()) ) } system.time(res <-‐ gen_rand()) ggplot() + layer(data=res, mapping=aes(x=V1), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + ggtitle("Bernoulli to Standard Uniform") 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒ ඪ४Ұ༷෼෍

Ұ༷෼෍ 4UBOEBSE6OJGPSN%JTUSJCVUJPO

֬཰ີ౓ؔ਺ ύϥϝʔλʔ BɿԼݶ Cɿ্ݶ f(x, a, b) = ⇢ (b
a) 1 (a  x  b) 0 (otherwise) B͔ΒCͷؒͰ͸౳֬཰Ͱൃੜ͢ΔΑ͏ͳݱ৅Λ ද͢෼෍ Ұ༷෼෍

ίΠϯ౤͛ දͷ֬཰Q Rճ࣮ࢪ Ұ༷෼෍ ฏߦҠಈͱҾ͖Ԇ͹͠

a <-‐ 5 b <-‐ 8; width <-‐
0.05 p <-‐ 0.5 sample_size <-‐ 1000 trial_size <-‐ 500000 gen_unif_rand <-‐ function() { # sample_sizeܻͷ2ਐগ਺ΛϕϧψʔΠ෼෍ʹ # ै͏ཚ਺͔Βੜ੒ return (sum(rbern(sample_size, p) * (rep(1/2, sample_size) ** seq(sample_size)))) } gen_rand <-‐ function(){ return( rdply(trial_size, gen_unif_rand()) ) } system.time(res <-‐ gen_rand()) res$V1 <-‐ res$V1 * (b-‐a) + a ggplot() + layer(data=res, mapping=aes(x=V1), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + ggtitle("Bernoulli to Uniform") + xlim(4,9) 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒ Ұ༷෼෍

ϕʔλ෼෍ #FUB%JTUSJCVUJPO

ϕʔλ෼෍ ֬཰ີ౓ؔ਺ ύϥϝʔλʔ Ћɿฒ΂ସ͑ͨ͋ͱͷαϯϓϦϯά͢Δϙδγϣϯ Ќɿಉ͘͡ϙδγϣϯ͕ͩɺޙΖ͔Β਺͑ͨͱ͖ͷ΋ͷ f ( x, ↵, )
= 1 B ( ↵, )x ↵ 1(1 x ) 1 (0 < x < 1) ಠཱʹඪ४Ұ༷෼෍ʹै͏ Ћ Ќ ݸͷ֬཰ม਺ Λେ͖͞ͷॱʹฒ΂ସ͑ͨͱ͖ɺখ͍͞ํ͔Β Ћ൪Ίʢେ͖͍ํ͔Β͸Ќ൪໨ʣͷ֬཰ม਺9ͷ෼ ෍͕ϕʔλ෼෍# Ћ Ќ ͱͳΔɻ ˠ# ɾ ͸ϕʔλؔ਺

ίΠϯ౤͛ දͷ֬཰Q Rճ࣮ࢪ ϕʔλ෼෍ େ͖͞ͷॱʹฒ΂ͯ ̎ͭΊͷ஋Λϓϩοτ ͭαϯϓϦϯά

width <-‐ 0.03; p <-‐ 0.5 digits_length <-‐ 30;
set_size <-‐ 3 trial_size <-‐ 30000 gen_unif_rand <-‐ function() { # digits_lengthܻͷ2ਐগ਺ΛϕϧψʔΠ෼෍ # ʹै͏ཚ਺͔Βੜ੒ return (sum(rbern(digits_length, p) * (rep(1/2, digits_length) ** seq(digits_length)))) } gen_rand <-‐ function(){ return( rdply(set_size, gen_unif_rand())$V1 ) } unif_dataset <-‐ rlply(trial_size, gen_rand, .progress='text') p <-‐ ceiling(set_size * 0.5); q <-‐ set_size -‐ p + 1 get_nth_data <-‐ function(a){ return(a[order(a)][p]) } disp_data <-‐ data.frame(lapply(unif_dataset, get_nth_data)) names(disp_data) <-‐ seq(length(disp_data)); disp_data <-‐ data.frame(t(disp_data)) names(disp_data) <-‐ "V1" x_range <-‐ seq(0, 1, 0.001) dens <-‐ data.frame(y=dbeta(x_range, p, q)*trial_size*width) ggplot() + layer(data=disp_data, mapping=aes(x=V1), geom="bar", stat = "bin", binwidth=width, fill="#6666ee", color="gray" ) + layer(data=dens, mapping=aes(x=x_range, y=y), geom="line", stat="identity", position="identity", colour="red" ) + ggtitle("Bernoulli to Beta") 3ίʔυ ੺ͷۂઢɿཧ࿦తͳ֬཰ ώετάϥϜɿཚ਺͔Βੜ੒ ϕʔλ෼෍

To be continued…

ࢀߟ ͜ͷεϥΠυͷ3ͷίʔυ IUUQTHJTUHJUIVCDPNNBUTVLFOBFEBEBFD ޻ՊܥͷͨΊͷ౷ܭ֓࿦ IUUQXXXBNB[PODPKQEQ ࢦ਺෼෍ͱϙΞιϯ෼෍ͷ͍͚ͳ͍ؔ܎ IUUQXXXTMJEFTIBSFOFUUFSBNPOBHJTT ϕʔλ෼෍ͷੑ࣭ IUUQXXXLXBOTFJBDKQIT[TVHBLVDUPVLFJCFUB CFUBIUN

全ての確率はコイン投げに通ず

全ての確率はコイン投げに通ず

More Decks by kenmatsu4

Other Decks in Science

Featured

Transcript