#----------------------------Questions:1------------------------------------#
###################### =n <- 10
x <- 4
betafun <- function(n,x,theta) {
postfun <- (gamma(n+2) * theta^x * (1-theta)^(n-x))/(gamma(x+1)*gamma(n-x+1))
return(postfun)
}
theta <- seq(0,1, length=100)
theta <- theta[c(-1,-100)]
beta <- c(1:length(theta))
for (j in 1:length(theta)) {
beta[j] <- betafun(n,x,theta[j])
}
plot(theta,beta)
###################### d
Mybeta <- function(alpha,beta,n,x) {
mybate <- function(alpha,beta,n,x,theta) {
posterior_PDF <- (gamma(alpha+beta+n) * theta^(alpha+x-1) * (1-theta)^(beta+n-x-1))/(gamma(alpha+x)*gamma(beta+n-x))
return(posterior_PDF)
}
theta <- seq(0,1,length = 100)
theta <- theta[c(-1,-100)]
mybateDate <- c(1:length(theta))
for (m in 1:length(theta)) {
theta1 <- theta[m]
mybateDate[m] <- mybate(alpha,beta,n,x,theta1)
}
plot(x = theta,y = mybateDate)
}
###################### e
Mybeta(alpha = 2, beta=2, n=10, x=6)
Mybeta(alpha = 4, beta=2, n=10, x=6)
Mybeta(alpha = 2, beta=4, n=10, x=6)
Mybeta(alpha = 20, beta=20, n=10, x=6)
mymix <- function(w,x,a1,b1,a2,b2) {
mixbeta <- w* x^(a1-1) * (1-x)^(b1-1)/beta(a1,b1) + (1-w)* x^(a2-1) * (1-x)^(b2-1)/beta(a2,b2)
return(mixbeta)
}
###################### question c
Mymixplot <- function(w,a1,b1,a2,b2) {
mymix <- function(w,x,a1,b1,a2,b2) {
mixbeta <- w* x^(a1-1) * (1-x)^(b1-1)/beta(a1,b1) + (1-w)* x^(a2-1) * (1-x)^(b2-1)/beta(a2,b2)
return(mixbeta)
}
x <- seq(0,1,length = 100)
x <- x[c(-1,-100)]
mixdata <- c(1:length(x))
for (k in 1:length(x)) {
mixdata[k] <- mymix(w,x[k],a1,b1,a2,b2)
}
plot(x,mixdata)
}
# plot
Mymixplot(w=0.3, a1=2,b1=4,a2=4,b2=2)
Mymixplot(w=0.5, a1=2,b1=4,a2=4,b2=2)
Mymixplot(w=0.7, a1=2,b1=4,a2=4,b2=2)
Mymixplot(w=0.9, a1=2,b1=4,a2=4,b2=2)
# posterior
mixPostPlot <- function(w,n,x,a1,b1,a2,b2) {
mybateOne <- function(a1,b1,n,x,theta) {
posterior_PDFOne <- (gamma(a1+b1+n) * theta^(a1+x-1) * (1-theta)^(b1+n-x-1))/(gamma(a1+x)*gamma(b1+n-x))
return(posterior_PDFOne)
}
mybateTwo <- function(a2,b2,n,x,theta) {
posterior_PDFTwo <- (gamma(a2+b2+n) * theta^(a2+x-1) * (1-theta)^(b2+n-x-1))/(gamma(a2+x)*gamma(b2+n-x))
return(posterior_PDFTwo)
}
theta <- seq(0,1,length = 100)
theta <- theta[c(-1,-100)]
mixdata <- c(1:length(theta))
for (k in 1:length(theta)) {
mixdata[k] <- w*mybateOne(a1,b1,n,x,theta[k]) + (1-w)*mybateTwo(a1,b1,n,x,theta[k])
}
plot(theta,mixdata)
}
mixPostPlot(0.5,10,6,2,4,2,4)
# prior
x <- seq(0,100,by = 1)
y <- dbinom(x,100,0.4)
plot(x,y)