確率分布の関数

• こちらで二項分布を使っている

choose((v - 1),(Nout - 1))*(0.1*p)^(v - Nout)*(1 - 0.1*p)^(Nout - 1)*(1 - 0.1*p)

• (v-1)から(Nout-1)を取り出す場合の数　ｘ　0.1pの確率で(v-Nout)回とって、(1-0.1p)の確率で(Nout-1)回とる、という計算のこと
• このように、全部でA回の試行をしたとき、確率pの事象がB回起こって、A-B回は起こらなかったときの確率分布が２項分布で、Rでは次のように計算する

dbinom(x=B,size=A,prob=p)

• ここで"binom"は二項分布を表す文字、"d"は確率『密度(density)』を表す接頭辞
• "binom"を付けた関数には以下がある
  • "dbinom()"
  • "pbinom()"
  • "qbinom()"
  • "rbinom()"
• 正規分布("norm"に対しては同様に"dnorm()""pnorm()""qnorm()""rnorm()"がある。以下色々な確率分布について、同様)
• さて、２項分布の関数"xbinom()"に戻ろう
  • A,B,pを与えたときに、その生起確率は"dbinom(x=B,size=A,prob=p)"が与える。これは、確率密度分布(distribution)の"d"

A<-10
B<-3
pr<-0.4
dbinom(x=B,size=A,prob=pr)

#Bは0,1,2,...,Aをとるからそう置き換えて
B<-0:A
dbinom(x=B,size=A,prob=pr)
# その和は１
sum(dbinom(x=0:A,size=A,prob=pr))

#dbinomの累積は
cumsum(dbinom(x=B,size=A,prob=pr))

#累積密度分布の関数は実は用意してあって
pbinom(q=B,size=A,prob=pr)

#じゃあ、累積確率pがこのくらい、とわかっているときにBの値が知りたいとしたら。。。
qbinom(p=0.5,size=A,prob=pr)

# pの値をいろいろ振ると
p<-seq(from=0,to=1,by=0.01)
qbinom(p=p,size=A,prob=pr)

plot(p,qbinom(p=p,size=A,prob=pr))

# 最後にrbinom()
# これは、prob=prで起きる事象がA回繰り返されたときに、何回起きるか(Bの値がいくつになるか)を、ランダムに試してみる関数
ntrial<-1000 # A回繰り返しを１０００回試してみた
rbinom(n=ntrial,size=A,prob=pr)
hist(rbinom(n=ntrial,size=A,prob=pr))

• おまけ

# 検討するパラメタとその範囲を決めるための指定
maxNout.num<-10 # 検討する最大アウト数
maxBatter.num<-20 # 検討する最大バッター数
kizami<-0.1 # 検討するヒット率(非アウト率)を0-1の間のいくつ刻みで考えるか

# 検討するパラメタの値のリストを作る
Nout.nums<-1:maxNout.num # 検討するアウト数のリスト
Batter.nums<-1:maxBatter.num # 検討するバッター数のリスト
hitps<-seq(from=0,to=1,by=kizami) # 検討するヒット率のリスト

# パラメタの数だけループを作って、すべての場合を検討する
# パラメタ数次元のarrayに格納する
result<-result2<-array(0,c(length(Nout.nums),length(Batter.nums),length(hitps)))
# resultの次元を確認
dim(result)
for(i in 1:length(Nout.nums)){
	Nout.num<-Nout.nums[i]
	for(j in 1:length(Batter.nums)){
		Batter.num<-Batter.nums[j]
		for(k in 1:length(hitps)){
			hitp<-hitps[k]
			if(Nout.num<=Batter.num){
				result[i,j,k]<-choose(Batter.num-1,Nout.num-1)*(1-hitp)^(Nout.num-1)*hitp^(Batter.num-(Nout.num-1))*(1-hitp)
				result2[i,j,k]<-dbinom(x=Nout.num-1,size=Batter.num-1,prob=(1-hitp))*(1-hitp)
			}
		}
	}
}

apply(result2[,,3],1,sum)

persp(result2[,,3])