#”a”と”b”の集団を作る
Na<-1;Nb<-11;k<-4
Ns<-Na+Nb
A<-c(rep("a",Na),rep("b",Nb))
A
#k倍する
B<-rep(A,k)
#Ns個抜き取る
sample(B,Ns)
Niter<-1000
#Niter回、繰り返して
Numa<-rep(0,Niter)
for(i in 1:Niter){
S<-sample(B,Ns)
Numa[i]<-length(which(S=="a"))
}
#"a"の抜き取られ数の分布をみる
hist(Numa)
#"a"が２個抜き取られた確率は？
length(which(Numa==2))/Niter
#n!を使って計算してみる
gamma(12+1)/(gamma(2+1)*gamma(10+1))*gamma(36+1)/(gamma(2+1)*gamma(34+1))*gamma(4+1)*gamma(44+1)/gamma(48+1)
#"a"の抜き取り数は、0,1,2,3,4のいずれか
#その確率を計算して足し合わせてみる
P0<-gamma(12+1)/(gamma(0+1)*gamma(12+1))*gamma(36+1)/(gamma(4+1)*gamma(32+1))*gamma(4+1)*gamma(44+1)/gamma(48+1)
P1<-gamma(12+1)/(gamma(1+1)*gamma(11+1))*gamma(36+1)/(gamma(3+1)*gamma(33+1))*gamma(4+1)*gamma(44+1)/gamma(48+1)
P2<-gamma(12+1)/(gamma(2+1)*gamma(10+1))*gamma(36+1)/(gamma(2+1)*gamma(34+1))*gamma(4+1)*gamma(44+1)/gamma(48+1)
P3<-gamma(12+1)/(gamma(3+1)*gamma(9+1))*gamma(36+1)/(gamma(1+1)*gamma(35+1))*gamma(4+1)*gamma(44+1)/gamma(48+1)
P4<-gamma(12+1)/(gamma(4+1)*gamma(8+1))*gamma(36+1)/(gamma(0+1)*gamma(36+1))*gamma(4+1)*gamma(44+1)/gamma(48+1)
#足して１になる
P0+P1+P2+P3+P4
