全問共通

#==========================================================
# File: 03.HW03.R
# Working directory: /Users/Jaehyun/Dropbox/Graduate School
#                    /Political Methodology I/2014(Yanai)
# Title: Research methods in Poltical Science - Homework 03
# Date: 2014-10-15
# Modified : 2014-10-19
# Author : Jaehyun Song
#==========================================================

# Call "ggplot2" due to creating graphs
require(ggplot2)
## Loading required package: ggplot2

問題1

  1. 関数の名前はJay.factorialと命名し、引数はxと指定する。
  2. 戻り値はreturnであり、掛け算のために最初の段階で1を割り当てる。
  3. if関数で条件付き分岐をする。
    1. 引数が0より小さい場合はエラーメッセージを出力する。(if(x < 0)の部分)
    2. 引数が0の場合は1を返する。(else if(x == 0)の部分)
    3. 引数が0より大きい場合はべき乗の演算(\(n! = n(n-1)(n-2)...1\))を行い、その結果を返す。(elseの部分)
#=================
# Probelm 1
#=================

#To define a function, called "Jay.factorial", which perfome a function same as "factorial" in R
Jay.factorial <- function(x){ 
  # Argument:x (integer)
  # Return value:result
  
  result <- 1 # To declare the return value

  if(x < 0){       # The case which the argument is less than zero
    #Display an error message
    print("The argument must be greater than ZERO.") 
  }
  
  else if(x == 0){ #If x equals to zero
    return(1) # Return 1
  }
  
  else{            # If the argument is greater than zero
    for(i in x:1){
      result <- result * i # Multiply "i" and the previous value, "result"
    }
    return(result) #Return the result
  }
}

最後にRの内蔵関数(factorial)と自作関数の結果を比較する。

#Comparison between the user defiend function, Jay.factorial, and factorial in R
factorial(10) == Jay.factorial(10)
## [1] TRUE

問題2

  1. 組み合わせを演算するユーザー定義関数をJay.chooseと命名し、因数はnkにする。これはn個のものからk個のものを選ぶ事を意味する。
  2. if関数で条件付き分岐をする。
    1. knより大きい場合、0を返す。(if(n < k)の部分)
    2. それ以外の場合は組合せの演算を行う。(elseの部分) 組合せの公式は\[\begin{pmatrix} n \\ k \end{pmatrix}=\frac { n! }{ k!(n-k)! }\]である。べき乗の演算はJay.factorial()で処理する。
    3. 結果はresultに書き込まれ、最後にresultの値を返す。(return(result)の部分)
#=================
# Probelm 2
#=================

#To define a function, called "Jay.choose", which perfome a function same as "choose" in R
Jay.choose <- function(n, k){ 
  # Arguments: n, k
  # Retrun value:result
  
  if(n < k){  # if k is greater than n
    return(0) # return "0"
  }
  
  else{       #If n is greater than k
    # Calculate combination
    result <- Jay.factorial(n) / (Jay.factorial(k) * Jay.factorial(n-k))
    
    return(result) # Return the result
  }
}

最後にRの内蔵関数(choose)と自作関数の結果を比較する。

#Comparison between the user defiend function, Jay.factorial, and factorial in R
choose(10, 5) == Jay.choose(10, 5)
## [1] TRUE

問題3

  1. 二項分布の確率質量を求めるユーザー定義関数をJay.dbinomと命名し、因数はxnpとする。xは成功回数であり、nは試行回数、pは成功確率である。
  2. 確率質量を求めるために\[f\left( x \right) =\begin{pmatrix} n \\ x \end{pmatrix}{ p }^{ x }{ (1-p) }^{ n-x }\quad (x=0,1,2,...,n)\]の演算を行う。組合せは問題2のJay.chooseを利用して演算する。
  3. 最後に結果をresultに書き込み、返す。
#=================
# Probelm 3
#=================

#To define a function, called "Jay.dbinom", which perfome a function same as "dbinom" in R
Jay.dbinom <- function(x, n, p){
  # Arguments: x - number of success (integer)
  #            n - number of bernoulli trials (integer)
  #            p - probability of success (integer)
  # Return Value: result

    # Calculate the binomial probability
    result <- Jay.choose(n, x) * p^x * (1-p)^(n-x)
  
  return(result)
}

最後にRの内蔵関数(dbinom)と自作関数の結果を比較する。

#Comparison between the user defiend function, Jay.factorial, and factorial in R
dbinom(1, 10, .5) == Jay.dbinom(1, 10, .5) # This result can be ignored
## [1] FALSE

比較の結果、内臓関数と自作関数の値が一致しないことが確認できる。しかし、これは丸めの問題である可能性が高いため、個別的に検証する。

dbinom(1, 10, .5)
## [1] 0.009766
Jay.dbinom(1, 10, .5)
## [1] 0.009766

以上の結果より、内臓関数と自作関数は少なくとも小数点9桁までは同じ数値を返すことが確認できる。


問題4

最初に二項分布の確率質量を書き込むためのデータフレームを生成する。ここでは最終的に使うデータフレームfinal.data.frameとループ文内で使用するtemp.data.frameを生成した。temp.data.frameは第1行のみを使い上書きするため欠損値(NA)で生成したが、final.data.framerbind()コマンドを用いて結合させるためNULLで設定する。

#=================
# Probelm 4
#=================

# Creat an empty data frame called "temp.data.frame" (Temporary data)
temp.data.frame <- data.frame(p=NA, n=NA, x=NA, result=NA)
# Creat an empty data frame called "final.data.frame" (Full data)
final.data.frame <- data.frame(p=NULL, n=NULL, x=NULL, result=NULL)

成功確率p, 試行回数n, 成功回数x、それぞれに対してループを使うため3つのfor文を用いる。これらは独立しているわけではなく、三重構造になる。3つ目のfor文の中では臨時データフレーム, temp.data.framep, n, x, 確率質量を書き込み、次はrbind()コマンドを使用してfinal.data.frameの最後業に追加する。確率質量は問題3のユーザー定義関数であるJay.dbninom()を利用する。

for(p in c(0.5, 0.8)){         # Loop for p = 0.5 and p = 0.8
  for(i in c(2, 3, 5, 10)){    # Loop for n = 2, 3, 5, 10
    for(j in 0:i){             # Loop for x = 0 ~ n
      
      # Input p, n, x, and result into the temporary data frame
      temp.data.frame[1, ] <- c(paste("p =", p), paste("n =", i), j, Jay.dbinom(j, i, p))
      
      # Combining final.data.frame with temp.data.frame
      final.data.frame <- rbind(final.data.frame, temp.data.frame) 
    }
  }  
}

この過程によって生成されたfinal.data.frameは中身は以下のようである。

# Display first six data in "final.data.frame"
head(final.data.frame)
##         p     n x result
## 1 p = 0.5 n = 2 0   0.25
## 2 p = 0.5 n = 2 1    0.5
## 3 p = 0.5 n = 2 2   0.25
## 4 p = 0.5 n = 3 0  0.125
## 5 p = 0.5 n = 3 1  0.375
## 6 p = 0.5 n = 3 2  0.375
# Display last six data in "final.data.frame"
tail(final.data.frame)
##          p      n  x       result
## 43 p = 0.8 n = 10  5 0.0264241152
## 44 p = 0.8 n = 10  6  0.088080384
## 45 p = 0.8 n = 10  7  0.201326592
## 46 p = 0.8 n = 10  8  0.301989888
## 47 p = 0.8 n = 10  9  0.268435456
## 48 p = 0.8 n = 10 10 0.1073741824

final.data.frameの中身は数値データではなく文字データとして書き込まれているため、このままggplotを用いて4つのグラフを表示させるとn=10, n=2, n=3, n=5の順番で表示される。これを避けるためにfinal.data.framen列のデータをfactor型へ変換し、順序付けする。

# Factoring the column, "n", in order to order the plots 
final.data.frame$n <- factor(final.data.frame$n, 
                             levels = c("n = 2", "n = 3", "n = 5", "n = 10"))

結果グラフを表示する。作図をする前に2つの図において共通的に使われるオプションを一つのオブジェクトに収める。そのためにリスト型のオブジェクト(prob4.opts)を生成し、その中にオプションを書き込む。

# Save option components into an object called "prob4.opts"
prob4.opts <- list(geom_bar(aes(x = as.integer(x), y = as.numeric(result)), stat="identity"), 
                   facet_wrap(~n), ylab("Probability"), xlab("Number of success"))

成功確率(p)に分けて作図するため、同一データフレーム内の異なる行のデータ(1 ~ 24行:p = 0.5, 24 ~ 48行:p = 0.8)を用いる。

# Display the plot of "p = 0.5" with dividing by "n"
ggplot(data=final.data.frame[1:24, ]) + prob4.opts + ggtitle("In The Case of p = 0.5")

plot of chunk prob4-6

以上はp = 0.5の場合の二項分布の確率質量である。

続いて、p = 0.8の場合の確率質量の棒グラフを表示する。

# Display the plot of "p = 0.8" with dividing by "n"
ggplot(data=final.data.frame[25:48, ]) + prob4.opts + ggtitle("In The Case of p = 0.8")

plot of chunk prob4-7