# Define x range
x <- seq(0, 1, length.out = 1000)

# Compute Beta densities
beta_11 <- dbeta(x, 1, 1)
beta_21 <- dbeta(x, 2, 1)
beta_33 <- dbeta(x, 3, 3)
beta_36 <- dbeta(x, 3, 6)
beta_48 <- dbeta(x, 4, 8)

# Plot the distributions
plot(x, beta_11, type = "l", lwd = 2,
     xlab = "x", ylab = "Density",
     main = "Beta Distribution PDFs",
     ylim = c(0, max(beta_21)))

lines(x, beta_21, lwd = 2)
lines(x, beta_33, lwd = 2)
lines(x, beta_36, lwd = 2)
lines(x, beta_48, lwd = 2)

# Add legend
legend("topright",
       legend = c("Beta(1,1)", "Beta(2,1)", "Beta(3,3)",
                  "Beta(3,6)", "Beta(4,8)"),
       lwd = 2)

# Define x range
x <- seq(0, 1, length.out = 1000)

# Compute Beta densities
beta_11 <- dbeta(x, 1, 1)
beta_21 <- dbeta(x, 2, 1)
beta_33 <- dbeta(x, 3, 3)
beta_36 <- dbeta(x, 3, 6)
beta_48 <- dbeta(x, 4, 8)

# Plot the first distribution
plot(x, beta_11, type = "l", lwd = 2, col = "black",
     xlab = "x", ylab = "Density",
     main = "Beta Distribution PDFs",
     ylim = c(0, max(beta_21)))

# Add remaining distributions with different colors
lines(x, beta_21, lwd = 2, col = "blue")
lines(x, beta_33, lwd = 2, col = "red")
lines(x, beta_36, lwd = 2, col = "darkgreen")
lines(x, beta_48, lwd = 2, col = "purple")

# Add legend
legend("topright",
       legend = c("Beta(1,1)", "Beta(2,1)", "Beta(3,3)",
                  "Beta(3,6)", "Beta(4,8)"),
       col = c("black", "blue", "red", "darkgreen", "purple"),
       lwd = 2)

mu <- 6
sigma <- 0.5
N <- 500000
# Generate N random samples from a log-normal distribution
sl <- rlnorm(N, mu, sigma)
ggplot(tibble(samples = sl), aes(samples)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 50) +
ggtitle("Log-normal distribution\n") +
coord_cartesian(xlim = c(0, 2000))

# Generate N random samples from a normal distribution,
# and then exponentiate them
sn <- exp(rnorm(N, mu, sigma))
ggplot(tibble(samples = sn), aes(samples)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 50) +
ggtitle("Exponentiated samples from\na normal distribution") +
coord_cartesian(xlim = c(0, 2000))