MCMCでパラメータθの推定がしたい。

前提・実現したいこと

自殺における相対リスクのベイズ推定という論文のRのコードを実行しようとしてもエラーが起きます。

発生している問題・エラーメッセージ

if (u[i] <= exp(logr)) { でエラー: 
   TRUE/FALSE が必要なところが欠損値です 
 追加情報:  警告メッセージ: 
 rgamma(1, shape = m * alpha[i - 1] + a1, rate = sumthe + b2) で: 
   NA が生成されました

該当のソースコード

data <- read.csv("hyougo.csv",header=T)
library(coda)
d <- data$d
e <- data$e
m <- nrow(data)
a1 <- 1
b1 <- 1
a2 <- 0.1
b2 <- 1
fa <- function(alpha,beta,theta,a,b){
  slthe <- sum(log(theta))
  logf <- (m*alpha)*log(beta) + (alpha-1)*slthe - b*alpha - m*lgamma(alpha)
  return(logf)
}
chain <- 3
alphst <- c(10,75,150)
betast <- c(10,75,150)
burn <- 10000
thin <- 10
Niter <- 102000
reall <- array(0, dim=c((Niter*chain), (m+3))) 
for (t in 1:chain){
  alpha <- numeric(Niter)
  beta <- numeric(Niter)
  theta <- array(0, dim=c(Niter, m))
  alpha[1] <- alphst[t]
  beta[1] <- betast[t]
  sigma <- 1.9
  u <- runif(Niter)
  k <- numeric(Niter)
  for (i in 2:Niter){
    for (j in 1:m){
      theta[i,j] <- rgamma(1, shape=d[j] + alpha[i-1], rate=e[j] + beta[i-1])
    }
    sumthe <- sum(theta[i,])
    beta[i] <- rgamma(1, shape=m*alpha[i-1] + a1, rate=sumthe + b2)
    y <- rnorm(1, alpha[i-1], sigma)
    logr <- fa(y, beta[i],theta[i,],a2,b2) - fa(alpha[i-1],beta[i],theta[i,],a2,b2)
    if(u[i] <= exp(logr)){
      alpha[i] <- y
      k[i] <- 1
    }
    else{
      alpha[i] <- alpha[i-1]
    }
  }
  result <- cbind(alpha=alpha,beta=beta,theta=theta,k=k)
  reall[((t-1)*Niter+1):(t*Niter),] <- result
}

補足情報（FW/ツールのバージョンなど）

データは
d　　　　e
2.2 2.4
3.4 3.4
5.7 5
7.4 11
9.3 12
15 13.4
18.5 12.8
28.6 25.8
45.3 39.8
78.4 66
です。
論文のURL:https://www.kansai-u.ac.jp/step/img/public/20140331/0331-1.pdf