### Gamma (two plots), F (two parameters: two plots) and t(one paramter: one plot)
### mean, variance, skewness, excess kurtosis


### Gamma Distribution
# Plot 1: Fix shape and changing rate
curve(expr = dgamma(x = x, shape = 2, rate = 2),
      xlab = "", ylab = "", main = "Gamma PDFs with r = 2",
      lwd = 2, col = 1, xlim = c(0, 2),
      ylim = c(0, 8))

for (i in 1:5) {
  curve(expr = dgamma(x = x, shape = 2, 
                      rate = c(3, 5, 10, 15, 20)[i]),
        lwd = 2, col = (2:6)[i], add = TRUE)
}

legend(x = "topright", legend = c("Rate = 2", "Rate = 3",
                                  "Rate = 5", "Rate = 10",
                                  "Rate = 15", "Rate = 20"),
       lwd = 2, col = 1:6)

# When fixing shape (r) and changing rate (lambda) for Gamma Distribution,
# When the rate (lambda) increases, based on the formula, 
# the mean decreases (shift to the left),
# the variance decreases
# the skewness does not change
# the excess kurtosis does not change


# Plot 2: Fix rate and changing shape
curve(expr = dgamma(x = x, shape = 2, rate = 2),
      xlab = "", ylab = "", main = "Gamma PDFs with lambda = 2",
      lwd = 2, col = 1, xlim = c(0, 20),
      ylim = c(0, 0.8))

for (i in 1:5) {
  curve(expr = dgamma(x = x, shape = c(3, 5, 10, 15, 20)[i], 
                      rate = 2),
        lwd = 2, col = (2:6)[i], add = TRUE)
}
legend(x = "topright", legend = c("Shape = 2", "Shape = 3",
                                  "Shape = 5", "Shape = 10",
                                  "Shape = 15", "Shape = 20"),
       lwd = 2, col = 1:6)

# When fixing rate (lambda) and changing shape (r) for Gamma Distribution,
# When the shape (r) increases, based on the formula, 
# the mean increases (shift to the right),
# the variance increases
# the skewness decreases
# the excess kurtosis decreases

### F Distribution
# Plot 1: Fix df2 and changing df1
curve(expr = df(x = x, df1 = 3, df2 = 10),
      xlab = "", ylab = "", main = "F PDFs with df2 = 10",
      lwd = 2, col = 1, xlim = c(0, 4),
      ylim = c(0, 1))

for (i in 1:5) {
  curve(expr = df(x = x, df1  = c(5, 15, 30, 60, 80)[i], 
                  df2=10),
        lwd = 2, col = (2:6)[i], add = TRUE)
}  
legend(x = "topright", legend = c("df1 = 3", "df1 = 5",
                                  "df1 = 15", "df1 = 30",
                                  "df1 = 60", "df1 = 80"),
       lwd = 2, col = 1:6)

# When fixing df2 and changing df1 for F distribution,
# When the df1 increases,
# According to the formula: m = df1, n= df2, we can test diffrent values to see the relationship
# mean <- n/(n-2) when n > 2
# variance <- (2*(n^2)*(m+n-2))/(m*(n-2)^2*(n-4)) when n > 4
# skewness <- ((2*m+n-2)*sqrt(8*(n-4)))/((n=6)*sqrt(m*(m+n-2))) when n > 6
# excess kurtosis <- 12*(m*(5*n-22)*(m+n-2)+(n-4)*(n-2)^2) / (m*(n-6)*(n-8)*(m+n-2)) when n > 8

# the mean does not change,
# the variance decreases
# the skewness decreases
# the excess kurtosis decreases

# Plot 2: Fix df1 and changing df2
curve(expr = df(x = x, df1 = 2, df2 = 3),
      xlab = "", ylab = "", main = "F PDFs with df1 = 2",
      lwd = 2, col = 1, xlim = c(0, 4),
      ylim = c(0, 1))

for (i in 1:5) {
  curve(expr = df(x = x, df1  = 2, 
                  df2=c(5, 15, 30, 60, 80)[i]),
        lwd = 2, col = (2:6)[i], add = TRUE)
} 
legend(x = "topright", legend = c("df2 = 3", "df2 = 5",
                                  "df2 = 15", "df2 = 30",
                                  "df2 = 60", "df2 = 80"),
       lwd = 2, col = 1:6)

# When fixing df1 and changing df2,
# When the df2 increases,
# According to the formula: m = df1, n= df2, we can test different values to see the relationship 
# mean <- n/(n-2) when n > 2
# variance <- (2*(n^2)*(m+n-2))/(m*(n-2)^2*(n-4)) when n > 4
# skewness <- ((2*m+n-2)*sqrt(8*(n-4)))/((n=6)*sqrt(m*(m+n-2))) when n > 6
# excess kurtosis <- 12*(m*(5*n-22)*(m+n-2)+(n-4)*(n-2)^2) / (m*(n-6)*(n-8)*(m+n-2)) when n > 8

# the mean decreases (shift to the left),
# the variance decreases
# the skewness increases
# the excess kurtosis decreases


### T distribution
curve(expr = dt(x = x, df = 3),
      xlab = "", ylab = "", main = "T PDFs with changing df",
      lwd = 2, col = 1, xlim = c(-4, 4),
      ylim = c(0, 0.5))

for (i in 1:5) {
  curve(expr = dt(x = x, df  = c(5, 15, 30, 60, 80)[i]),
        lwd = 2, col = (2:6)[i], add = TRUE)
}  
legend(x = "topright", legend = c("df = 3", "df = 5",
                                  "df = 15", "df = 30",
                                  "df = 60", "df = 80"),
       lwd = 2, col = 1:6)

# When changing df for t distribution,
# When the df increases, based on the formula, 
# the mean does not change, when df > 1
# the variance decreases, when df > 2
# the skewness does not change, when df > 3
# the excess kurtosis decreases, when df > 4