################################################################################
#Beta distribution examples #
################################################################################
################################################################################
#Beta distribution
a <- 0.5
b <- 2
#Use as ad-hoc way to get an upper limit for the plot
ylim.max <- max(dbeta(x = 0.01, shape1 = a, shape2 = b), dbeta(x = 0.99, shape1 = a, shape2 = b))
par(xaxs = "i") #Remove the extra 4% around x-axis limits that are put there by default
curve(expr = dbeta(x, shape1 = a, shape2 = b), xlim = c(0,1), ylab = "f(x)", xlab = "x",
ylim = c(0, ylim.max), col = "red",
main = substitute(paste("Beta distribution with ", alpha == a, " and ", beta == b), list(a = a, b = b)))
abline(h = 0)
#Density to integrate - because dbeta() already exists,
# do not really need to program in like this. However,
# this will be useful to see for later.
beta.density <- function(x, a, b) {
#1/beta(a = a, b = b) * x^(a-1) * (1-x)^(b-1) #Could program in density directly, but may be more likely to lead to a numerical error
dbeta(x = x, shape1 = a, shape2 = b)
}
#Integrates to 1 over entire region
integrate(f = beta.density, lower = 0, upper = 1, a = a, b = b)
#Integrate over a specific region
lower.bound <- 0.1
upper.bound <- 0.5
integrate(f = beta.density, lower = lower.bound, upper = upper.bound, a = a, b = b)
#Check
pbeta(q = upper.bound, shape1 = a, shape2 = b) - pbeta(q = lower.bound, shape1 = a, shape2 = b)
#Simple lines on plot showing limits
#segments(x0 = lower.bound, y0 = 0, x1 = lower.bound, y1 = dbeta(x = lower.bound, shape1 = a, shape2 = b),
# lty = "dashed")
#segments(x0 = upper.bound, y0 = 0, x1 = upper.bound, y1 = dbeta(x = upper.bound, shape1 = a, shape2 = b),
# lty = "dashed")
#Shaded area on plot - notice that duplicate x values at the beginning and end of vector in polygon().
# This is needed to correspond to the lower y limits
x.values <- seq(from = lower.bound, to = upper.bound, by = 0.001)
polygon(x = c(lower.bound, x.values, upper.bound), border = NA,
y = c(0, dbeta(x.values, shape1 = a, shape2 = b),0), col = "blue")
#Can use border = "blue" too, but f(x) is then covered by blue
#Show how to save results
ans <- integrate(f = beta.density, lower = lower.bound, upper = upper.bound, a = a, b = b)
ans
names(ans)
ans$value
ans$subdivisions
ans$message
#E(X)
beta.density.mu <- function(x, a, b) {
x * dbeta(x = x, shape1 = a, shape2 = b)
}
mu.x <- integrate(f = beta.density.mu, lower = 0, upper = 1, a = a, b = b)
mu.x
a/(a+b)
#Var(X)
beta.density.var <- function(x, a, b, mu) {
(x - mu)^2 * dbeta(x = x, shape1 = a, shape2 = b)
}
var.x <- integrate(f = beta.density.var, lower = 0, upper = 1, a = a, b = b, mu = mu.x$value)
var.x
a*b/ ((a+b)^2 * (a+b+1))
################################################################################
# Trapezoidal rule
library(pracma)
x.values2 <- seq(from = lower.bound, to = upper.bound, by = 0.2)
trapz(x = x.values2, y = dbeta(x = x.values2, shape1 = a, shape2 = b))
#For two trapezoids:
#All parts together: area = 0.5 * width * (height1 + height2)
0.5 * 0.2 * (dbeta(x = 0.1, shape1 = a, shape2 = b) + dbeta(x = 0.3, shape1 = a, shape2 = b)) +
0.5 * 0.2 * (dbeta(x = 0.3, shape1 = a, shape2 = b) + dbeta(x = 0.5, shape1 = a, shape2 = b))
0.5 * 0.2 * dbeta(x = 0.1, shape1 = a, shape2 = b) + 1 * 0.2 * dbeta(x = 0.3, shape1 = a, shape2 = b) +
0.5 * 0.2 * dbeta(x = 0.5, shape1 = a, shape2 = b)
#Alternative ways to write this that will be helpful shortly: SUM(w_i * f(x_i))
w1 <- 0.5 * 0.2
w2 <- 1 * 0.2
w3 <- 0.5 * 0.2
sum(w1*dbeta(x = 0.1, shape1 = a, shape2 = b) + w2*dbeta(x = 0.3, shape1 = a, shape2 = b) +
w3*dbeta(x = 0.5, shape1 = a, shape2 = b))
curve(expr = dbeta(x, shape1 = a, shape2 = b), xlim = c(0,1), ylab = "f(x)", xlab = "x",
ylim = c(0, ylim.max), col = "red",
main = substitute(paste("Beta distribution with ", alpha == a, " and ", beta == b), list(a = a, b = b)))
abline(h = 0)
#Draw trapezoids
segments(x0 = x.values2, y0 = 0, x1 = x.values2, y1 = dbeta(x = x.values2, shape1 = a, shape2 = b),
lty = "dashed")
numb.val <- length(x.values2)
segments(x0 = x.values2[1:(numb.val-1)], y0 = dbeta(x = x.values2[1:(numb.val-1)], shape1 = a, shape2 = b),
x1 = x.values2[2:numb.val], y1 = dbeta(x = x.values2[2:numb.val], shape1 = a, shape2 = b),
lty = "dashed")
################################################################################
# Gauss-Legendre quadrature
options(width = 60)
library(pracma)
beta.density2 <- function(x, alpha, beta) {
dbeta(x = x, shape1 = alpha, shape2 = beta)
}
save.int <- quadgk(f = beta.density2, a = lower.bound, b = upper.bound, alpha = a, beta = b)
save.int
#3 point
save.xw <- gaussLegendre(n = 3, a = lower.bound, b = upper.bound)
save.xw$x
save.xw$w
sum(save.xw$w * dbeta(x = save.xw$x, shape1 = a, shape2 = b))
#7 point
save.xw <- gaussLegendre(n = 7, a = lower.bound, b = upper.bound)
save.xw$x
save.xw$w
sum(save.xw$w * dbeta(x = save.xw$x, shape1 = a, shape2 = b))
#Another way to do same as above
library(gaussquad)
leg.rule <- legendre.quadrature.rules(n = 3)
leg.rule #x's and w's are not same as above because it assumes a range of -1 to 1
#It appears that legendre.quadrature() transforms f(x) to correspond to the -1 to 1 limits.
legendre.quadrature(functn = beta.density2, rule = leg.rule[[3]], lower = lower.bound, upper = upper.bound,
alpha = a, beta = b)
#Relate save.xw to leg.rule
save.xw <- gaussLegendre(n = 3, a = lower.bound, b = upper.bound) #Re-do to make sure n = 3
#Using the u-substitution formula:
(2*0.14508-(0.1+0.5))/0.4
(2*save.xw$x - (lower.bound + upper.bound))/(upper.bound - lower.bound)
#Take into account different sized interval
save.xw$w * 2/(upper.bound - lower.bound)
################################################################################
# Other quadrature
#This will not work because these polynomials are defined over -Inf to +Inf
# Would need to transform our integration to be over this range
save.xw.hermite <- gaussHermite(n = 3)
save.xw.hermite$x
save.xw.hermite$w
sum(save.xw.hermite$w * dbeta(x = save.xw.hermite$x, shape1 = a, shape2 = b))
#