ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2016-11-13 複素数とガンマ関数・ベータ分布

[][][][]複素数とガンマ関数・ベータ分布

library(fAsianOptions)

cgamma(1i)

my.cbeta <- function(a,b){
	ret <- cgamma(a)*cgamma(b)/cgamma(a+b)
	return(ret)
}

my.dcbeta <- function(x,a,b){
	x^(a-1) * (1-x)^(b-1) / my.cbeta(a,b)
}

x <- seq(from=-1,to=2,length=100)
a <- -0.5
b <- -0.5
y <- my.dcbeta(x,a,b)

plot(x,y)

2014-06-25 この先、どうなる?

ryamada222014-06-25

[][][]この先、どうなる?

  • 今、N回の観測をして、○がn1,×がn2回だったとする
  • この後、m回繰り返したら、m回のうち、何回成功するんだろう?
  • まず観測データから、○確率pの確率密度関数をベータ関数で表してPr(p) = ¥beta(n1+1,n2+1)=¥frac{p^{n1}(1-p)^{n2}}{B(n1+1,n2+1)}、ただし、B(a,b)はベータ関数
  • 今、○確率がpのときにm回のうちk回成功するのは¥begin{pmatrix}m¥¥k¥end{pmatrix} p^k (1-p)^{m-k}なので
  • 結局¥int_0^1¥begin{pmatrix}m¥¥k¥end{pmatrix} p^k (1-p)^{m-k} ¥times Pr(p) dp
  • これを式変形すると
  • ¥begin{pmatrix}m¥¥k¥end{pmatrix}¥frac{B(k+n1+1,m-k+n2+1)}{B(n1+1,n2+1)}¥int_0^1 ¥frac{p^{k+n1}(1-p)^{m-k+n2}}{B(k+n1+1,m-k+n2+1)}となるが、積分は1なので
  • ¥begin{pmatrix}m¥¥k¥end{pmatrix}¥frac{B(k+n1+1,m-k+n2+1)}{B(n1+1,n2+1)}
  • これって、なんかの名前が付いているはずだと思うのだが…はて?

f:id:ryamada22:20140623153624j:image

# n1の○、n2の×の下で、後m回
my.partial.binom <- function(m,n1,n2){
# ○は0,1,...,mになりえる
	M <- 0:m
	a <- lchoose(m,M)
	b <- lbeta(n1+1,n2+1)
	B <- lbeta(M+n1+1,M[length(M):1]+n2+1)
	exp(a+B-b)
}

plot(my.partial.binom(10,2,3))
plot(my.partial.binom(13,2,3))
# 既観察n1,n2
n1 <- 2
n2 <- 4
# 既観察とこれからの観察の総観察数N
N <- (n1+n2):30
# 追加試行の打ち○の回数別の確率のデータ格納
X <- matrix(0,length(N),max(N)+1)
# 別の格納庫
Ps <- list()
# ○率p
p <- seq(from=0,to=1,length=100)
# Nの値ごとに、○率pがどれくらいの確率になるかの格納
ps <- matrix(0,length(N),length(p))
for(i in 1:length(N)){
	tmp <- my.partial.binom(N[i]-(n1+n2),n1,n2)
	X[i,1:length(tmp)] <- tmp
	Ps[[i]] <- tmp
	for(j in 1:length(tmp)){
		ps[i,] <- ps[i,] + Ps[[i]][j] * dbeta(p,n1+(j-1)+1,n2+N[i]-(n1+n2)-(j-1)+1)
	}
	#print(sum(tmp*(0:(length(tmp)-1))))/(N-(n1+n2))
}

matplot(0:max(N),t(X),type="l")
# Nが変わっても○率の分布は変わらない
matplot(p,t(ps),type="l")
# Nが変わっても、○率の期待値は変わらない
for(i in 1:length(N)){
	print(sum(ps[i,] * p))
}