読者です 読者をやめる 読者になる 読者になる

Data Science by R and Python

統計学を、広く、深く、わかりやすく。

Rでベータ関数を描いた

特に、何か面白いことをしたわけではないですが、
ベータ関数を描いたので、コード載せとく。

logf = function(x){
	rr = 0
	for(i in 1:x){
		rr = rr + log(i)
	}
	return(rr)
}

#beta distibutionの点の値を求める関数
d_beta = function(theta,n=129,y=118){
	rr = logf(n+1) - logf(y) - logf(n-y) + y*log(theta) + (n-y)*log(1-theta)
	mle = y/n
	print(paste("Theta's MLE is",mle,sep=" "))
	return(exp(rr))
}

theta_mle = 118/129
mle_likeli = purpose(theta_mle)

x = seq(0,1,0.0001)
y = d_beta(x)
plot(x,d_beta(x),type="l",lty=1,xlab=expression(theta),ylab=expression(paste("p(",theta,"|",y[1],",...,",y[129],")")))
lines(c(theta_mle, theta_mle), c(mle_likeli, -1), col = "red4", lwd = 1.5,lty=2)
text(theta_mle-0.15, 2.3, expression(paste(hat(theta)[MLE])), adj = 0, col = "red4", cex = 1.8)
arrows(x0 = theta_mle-0.1, y0 = 1.5, x1 = theta_mle-0.05, y1 =  -0.05, col = "red4", lwd = 1,angle=15,length=0.1)
mtext(expression(paste("0.915")), side = 1, line = 2., at = theta_mle, adj = 0.5, col = "red4", cex = 1.0, padj = -2)

#背景が必要な場合
library(png)
ima = readPNG("~/Desktop/code.png")
lim <- par()
rasterImage(ima, lim$usr[1], lim$usr[3], lim$usr[2], lim$usr[4])
grid()
points(x,purpose(x),type="l",lty=1,col="black")
#segments(theta_mle, 0, theta_mle, mle_likeli,lty=2)
#segments(0, mle_likeli, theta_mle, mle_likeli,lty=2)