Question 1

I write R codes to calculate factorials.

cal_factorial <- function(n) { ## function to calculate the factorial
  ## Aarguments:
  ##  n = non-negative integer
  if (n > 0) {
    X <- seq(1, n, by=1) 
    ## calculate the permutation of choosing n items among n items
    Y <- prod(X)
  }
  if (n == 0) { ## denote 0!=1
    Y <- 1
  }
  if (n < 0) { ## print error message if the input value is wrong
    Y <- ("You can not calculate the factorial.")
  }
  return(Y)
}

Let’s use this function. When n=3, run:

cal_factorial(3)
## [1] 6

I got the result. Then, I try using the built-in function factorial().

factorial(3)
## [1] 6

I got the same result.


Question 2

I write R codes to calculate the number of possible combinations that can be obtained by taking a subset of elements from a larger set.

combination <- function(n, k) { ## function to calculate the number of possible combinations
  if (n >= k){
    ## use the function we made before
    ## calculate the numerator
    numerator <- cal_factorial(n)
    ## calculate the denominator
    denominator <- cal_factorial(k) * cal_factorial(n-k)
    Z <- numerator / denominator
  }
  if (n < k){ ## print error message if the input value is wrong
    Z <- ("You can not calculate the number of possible combinations")
  }
  return(Z)
}

Let’s use this function. When n=6, run:

combination(6, 3)
## [1] 20

I got the result. Then, I try using the built-in function choose().

choose(6, 3)
## [1] 20

I got the same result.


Question 3

I write R codes to calcilate the probability mass of binomial distributions.

binomial_dis <- function(x, n, p) { ## function to calculate the probability mass of binomial distributions
  ## Aarguments:
  ##  x = vector
  ##  n = the number of independent Bernoulli trials
  ##  p = the probability of success for each trial
  A <- combination(n, x)　## use the function we made before
  B <- p^x
  C <- (1 - p)^(n - x)
  D <- A * B * C
  return(D)
}

Let’s use this function. When x=2, n=5, p=0.5 run:

binomial_dis(x=2, n=5, p=0.5)
## [1] 0.3125

I got the result. Then, I try using the built-in function dbinom().

dbinom(x=2, size=5, p=0.5) 
## size means the number of independent Bernoulli trials
## [1] 0.3125

I got the same result.


Question 4

I use my own binomial PMF, display graphically the PMF of binomial distributions for n = 2, 5, 10, and 50. I think when substituting 0.3, 0.5, 0.7 for p.

n <- c(2, 5, 10, 50) ## substitute 2, 5, 10 and 50 for n
p <- c(0.3, 0.5, 0.7) ## substitute 0.3, 0.5 and 0.7 for p
par(oma = c(0, 0, 2, 0))
## write the titele on the upper part of screen
par(mfrow=c(2,2)) ## set the tetrameric screen for plots
for(i in 1:length(p)) { ## calculate the PMF of the binomial distribution when p=0.3, 0.5, 0.7 individually
  for(j in 1:length(n)) { ## calculate the PMF of the binomial distribution when n=2, 5, 10, 50 individually
    x <- 0:n[j] ## x = 0,1,2,...,n
    ## prepare vector to save the results
    distributions <- rep(NA, length(x))
    for(k in 1:length(x)) {
      ## use the function we made before
      ## save the k-th result in the k-th row
      distributions[k] <- binomial_dis(x[k], n[j], p[i])
    }
    ## make the plots
    barplot(distributions, names=c(0:n[j]),
            ylim=c(0, max(distributions)), xlab="x", ylab="f(x)",
            main=paste("n=", n[j]))
  }
  ## make the title of the plots
  mtext(text = paste("p=", p[i]),
        side=3, outer=T, cex=1.3)
}