world

Q1

factorial() is R’s equipped function of calculates factorial of argument.

Q1 requires me create same function.

For example,

factorial(41)
## [1] 3.345253e+49

I create two tyep of function to calculate factorial.

First type is this.

f0a <- function(n){ ## function of factorial
    a <- c(1:n) ## create vector which size is n
    b <- prod(a)    ## multiply all elements in vector(a)
    return(b)
}

check this

f0a(41)
## [1] 3.345253e+49

Second type is this.

f0b <- function(n){ ## function of factorial
    x <- n
    for (i in 1:(n-1)){ ## needs n-1 multiply not n
        x <- x*(n-1)
        n <- n-1    ## update the n 
        }
    return(x)
}

also

f0b(41)
## [1] 3.345253e+49
f0b(0)
## [1] 0

Incedentally, factorial of 0 is 1 by the definition.

So I have to insert this.

f1 <- function(n){  ## true function of factorial
    if(n==1){   ## by the definition
        return(1)
    }else{
        return(f0b(n))  ## call my factorial function
    }
}
f1(41)
## [1] 3.345253e+49
f1(0)
## [1] 0

f1() was completed as factorial function.

Q2

I create function of calculate the number of combinations.

Formula is below, and it is implemented as choose().

\[_nc_k=\frac{n!}{k!(n-k)!}\]

f2 <- function(n,k){
    if(n==k){   ## this is take all
        return(1)   ## its only one pattern
    }else{
        if(k==0){
            return(1)
        }else{
            a <- f1(n)/f1(n-k)
            b <- a/f1(k)
            return(b)
        }
    }
}

Strictly speaking,I should suppose k larger than n(when the number of combination is 0).

But person who knows f2()’s role will not input larger k.

Check this.

choose(10,6)
## [1] 210
f2(10,6)
## [1] 210

Q3

\[f(x)=_nc_x p^x (1-p)^{n-x}\]

According to below formula, I create function to calculate the PMF of binomial distributions.

aruguments is threee, the number of trials(n) and success(x), probability of success(p).

f3 <- function(n,x,p){
    a <- f2(n,x)
    b <- p^x
    c <- (1-p)^(n-x)
    d <- a*b*c  ## I had taken a long time to notice order of dbinom is (x,n,p) not (n,x,p)
    return(d)   ## therefore this function is uselessly polite
}

Check this.

dbinom(6,10,0.46)
## [1] 0.169177
f3(10,6,0.46)
## [1] 0.169177

Q4

Basically, (under the holds constant p,) Q4 requires probability mass as function of the number of success(n).

f4 <- function(n,p){    ## if coding 0.5 instead p, p is held as constant.
    book1 <- rep(NA,n+1)    ## prepare null vector to input answer
    hit <- 0    ## the number of success is hit
    if(hit > n){
        return(book1)   ## when hit beyonds n, calculates are not executed
    }else{
        for(i in 1:(n+1)){  ## answer vector is larger than n because of nC0=1 pattern
            book1[i] <- f3(n,hit,p) ## inputting answer according to order
            hit <- hit + 1
        }
    }
    return(book1)
}

Check this.

dbinom(c(0.4),4,0.46)
## [1] 0
f4(4,0.46)
## [1] 0.08503056 0.28973376 0.37021536 0.21024576 0.04477456

Eventually, I need to draw graph.

For the sake of simplicity, function is separated getting answer vector function and histogram drawing function.

If integrates f4a and zu, get error.

Error: evaluation has nested too deeply. Or infinite recursion

zu <- function(n){  ### drawing graph
    a <- length(n)
    barplot(n, col = "gray",axes=TRUE,
    xlab = "The number of success", ylab = "Probability", main ="")
    axis(1,0:a,labels = 0:a)
}

Function is complete once. so I drawing probability mass of n=2,5,10,50 by two different p.

n=2, p=0.23

zu(f4(2,0.23))

n=5,p=0.23

zu(f4(5,0.23))

n=10,p=0.23

zu(f4(10,0.23))

n=50,p=0.23

zu(f4(50,0.23))

n=2,p=0.6826

zu(f4(2,0.6826))

n=5,p=0.6826

zu(f4(5,0.6826))

n=10,p=0.6826

zu(f4(10,0.6826))

n=50,p=0.6826

zu(f4(50,0.6826))

end