Beta Distribution

Author

Hubert Baechli

Loading some Packages

Code

library(fitdistrplus)

Lade nötiges Paket: MASS

Lade nötiges Paket: survival

Code

library(survival)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ dplyr::select() masks MASS::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Code

# set.seed(1)

Generating Betakurve Parameter

Code

beta_par <- function(mu, sigma) {
  mu <- max(mu, 0.0001)
  mu <- min(mu, 0.9999)
  sigma <- max(sigma, 0.00001)
  par <- c(sqrt(mu * (1 - mu))*0.99, sigma)
  sigma <- min(par)
  a_b <- (mu * (1 - mu)) / sigma^2 - 1  # calc alpha + beta
  a <- mu * a_b
  b <- (1 - mu) * a_b
  Output <- list(alpha = a,
                   beta = b)
  return(Output)
}

sigma <- 0.1
mu <- 0.4

beta_par(mu = mu, sigma = sigma)

$alpha
[1] 9.2

$beta
[1] 13.8

Generating random values

Code

gen_rbeta <- function(n = n,
                      mu = mu, 
                      sigma = sigma
                      ) {
  par <- beta_par(mu = mu, sigma = sigma)
  vec <- rbeta(n = 1000, shape1 = par$alpha, shape2 = par$beta)

  return(vec)
}

Testing for different Parameters

Code

n <- 1000
m <- 10

data01 <- data.frame(ID = seq(1,n))  

# Weitere Spalten hinzufügen
for (i in seq(0, 10)) {
  temp <- round(gen_rbeta(n = n, mu = i / m, sigma = 0.1), 3)
  temp <- sort(temp)
  name <- sprintf("Set%03d", i)
  data01[[name]] <- temp
}
summary(data01$Set000)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0

Code

data001 <- pivot_longer(data01,
                     cols = starts_with("Set"),
                     names_to = "Set",
                     values_to = "Values")

ggplot(data = data001, aes(x = ID, y = Values, color = Set))  +
  geom_line() +
  ggtitle("") +
  xlab("number") +
  labs(color = "Data Sets") +
  theme_light() +
  theme()

Code

data02 <- data.frame(ID = seq(1,n))  

# Weitere Spalten hinzufügen
for (i in seq(0, 10)) {
  temp <- round(gen_rbeta(n = n, mu = i / m, sigma = 0.5), 3)
  temp <- sort(temp)
  name <- sprintf("Set%03d", i)
  data02[[name]] <- temp
}
data002 <- pivot_longer(data02,
                     cols = starts_with("Set"),
                     names_to = "Set",
                     values_to = "Values")

ggplot(data = data002, aes(x = ID, y = Values, color = Set))  +
  geom_line() +
  ggtitle("") +
  xlab("number") +
  labs(color = "Data Sets") +
  theme_light() +
  theme()

Code

data03 <- data.frame(ID = seq(1,n))  

# Weitere Spalten hinzufügen
for (i in seq(0, 10)) {
  temp <- round(gen_rbeta(n = n, mu = 0.5, sigma = i / m), 3)
  temp <- sort(temp)
  name <- sprintf("Set%03d", i)
  data03[[name]] <- temp
}
data003 <- pivot_longer(data03,
                     cols = starts_with("Set"),
                     names_to = "Set",
                     values_to = "Values")

ggplot(data = data003, aes(x = ID, y = Values, color = Set))  +
  geom_line() +
  ggtitle("") +
  xlab("number") +
  labs(color = "Data Sets") +
  theme_light() +
  theme()

Code

data04 <- data.frame(ID = seq(1,n))  

# Weitere Spalten hinzufügen
for (i in seq(0, 10)) {
  temp <- round(gen_rbeta(n = n, mu = 0.1, sigma = i / m), 3)
  temp <- sort(temp)
  name <- sprintf("Set%03d", i)
  data04[[name]] <- temp
}
data04 <- pivot_longer(data04,
                     cols = starts_with("Set"),
                     names_to = "Set",
                     values_to = "Values")

ggplot(data = data04, aes(x = ID, y = Values, color = Set))  +
  geom_line() +
  ggtitle("") +
  xlab("number") +
  labs(color = "Data Sets") +
  theme_light() +
  theme()

Fit Betaparameter

Code

datah <- data03[, "Set002"]
beta_par(mu = 0.5, sigma = 2/10)

$alpha
[1] 2.625

$beta
[1] 2.625

Code

hist(datah)

Code

fit <- fitdist(datah, "beta", start = list(shape1 = 0.5, shape2 = 0.2))
fit

Fitting of the distribution ' beta ' by maximum likelihood 
Parameters:
       estimate Std. Error
shape1 2.882962  0.1243341
shape2 2.796242  0.1202756