world
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.
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
\[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
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