update with beta distribution description

08-Intro-to-Bayes.Rmd

@@ -1,4 +1,4 @@
-# (PART) 貝葉斯統計 {-}
+# (PART) 貝葉斯統計 Introduction to Bayesian Statistics {-}
 
 # 貝葉斯統計入門 {#intro-Bayes}
 
@@ -507,6 +507,12 @@ pi <- Vectorize(function(theta) dbeta(theta, 2,8))
 curve(pi, xlab=~ theta, ylab="Density", main="Beta prior: a=2, b=8",frame=FALSE, lwd=2)
 ```
 
+
+你會發現，beta分佈的圖形特徵由它的兩個超參數 (hyperparameter) `a, b` 決定。當相對地 `a` 比較大的時候，beta分佈的概率多傾向於較靠近橫軸右邊，也就是1的位置（較高概率），相對地 `b` 比較大的時候，beta分佈的曲線多傾向於靠近橫軸左邊，也就是0的位置（較低概率）。如果 `a` 和 `b`同時都在變大的話，beta分佈的曲線就變得比較“瘦”。這樣決定概率分佈形狀的參數，又被叫做**形狀參數 (shape parameters)**
+
+
+
+
 我們定義 $a>0$ 時[伽馬方程](https://zh.wikipedia.org/wiki/%CE%93%E5%87%BD%E6%95%B0)爲
 
 $$\Gamma(a)=\int_0^\infty x^{a-1}e^{-ax}\text{ d}x$$

12-Bayesian-stats.Rmd

@@ -2685,13 +2685,13 @@ Dat <- list(
 
 # fit use rjags package
 
-inits <- function (){
-  list (beta0=rnorm(1), beta1=runif(1) )
-}
+# inits <- function (){
+#   list (beta0=rnorm(1), beta1=runif(1) )
+# }
 
 post.jags <- jags(
   data = Dat,
-  inits = inits,
+  # inits = inits,
   parameters.to.save = c("beta0", "beta1", "theta[6]"),
   n.iter = 5000,
   model.file = paste(bugpath,
@@ -2866,7 +2866,7 @@ jagsModel <- jags.model(
                                   sep = ""),
                      data = Dat,
                      n.chains = 2, 
-                     inits = inits,
+                     # inits = inits,
                      quiet = TRUE)
 # Step 2 update 10000 iterations
 
@@ -2904,7 +2904,7 @@ Dat <- list(
 )
 post.jags <- jags(
   data = Dat,
-  inits = inits,
+  # inits = inits,
   parameters.to.save = c("beta0", "beta1", "theta[6]"),
   n.iter = 20000,
   model.file = paste(bugpath,
@@ -2960,17 +2960,19 @@ model{
 ```
 
 
-```{r OpenBUGSPractical0424, message=TRUE, warning=FALSE, echo=TRUE}
+```{r OpenBUGSPractical0424, message=TRUE, warning=FALSE, echo=TRUE, eval=TRUE, cache=TRUE}
+
 # R2JAGS codes:
 Dat <- list(
      y = c(1, 3, 6, 8, 11, 15, 17, 19), 
      n = c(20, 20, 20, 20, 20, 20, 20, 20), 
      x = c(30, 32, 34, 36, 38, 40, 42, 44), 
      N = 8 
 )
+
 post.jags <- jags(
   data = Dat,
-  inits = inits,
+  # inits = inits,
   parameters.to.save = c("ED95", "OR", "P35", "beta0", "beta1"),
   n.iter = 20000,
   model.file = paste(bugpath,
@@ -2986,6 +2988,8 @@ post.jags <- jags(
 print(post.jags)
 ```
 
+
+
 所以，藥物每增加劑量1mg，有療效的比值比是OR = 1.46 (95%CrI: 1.31, 1.63)。能夠達到95%患者都有療效的劑量是 45.02 mg (95% CrI: 42.7, 48.2 mg)。如果給予患者藥物劑量爲 35 mg，患者的疼痛能夠得到緩解(有療效)的概率是 32.3% (95% CrI: 22.7%, 42.3%)。跟着看到這裏的你是不是也覺得貝葉斯的結果和過程能夠更加豐富地回答我們想知道的問題呢？ 
 
 ```{r R-OpenBUGSPractical0425, cache=FALSE, fig.width=3.5, fig.height=5, fig.cap='Posterior density plots of ED95, OR, and P35.', fig.align='center', out.width='80%', message=TRUE, warning=FALSE}
@@ -3216,7 +3220,31 @@ list(alpha = 10, beta = 0, gamma = -5, logsigma2 = 5)
 一開始我們先採集1000個事後概率分布樣本。模型中的四個未知參數$\alpha, \beta, \gamma, \sigma^2$的1000次事後MCMC採樣的歷史痕跡圖繪制如圖 \@ref(fig:BayesianChapter0505)。我們可以看到模型收斂的速度很快。刨除前50次採樣的圖 (Fig. \@ref(fig:BayesianChapter0506))，可以對採樣過程看得更加清楚。像圖 \@ref(fig:BayesianChapter0506) 這樣粗粗的有點像毛毛蟲一樣的歷史痕跡圖通常象徵已經達到理想的收斂。
 
 
-```{r BayesianChapter0504, cache=TRUE, message=TRUE, warning=FALSE, echo=TRUE}
+<!-- ```{r BayesianChapter0504, cache=TRUE, message=TRUE, warning=FALSE, echo=TRUE} -->
+<!-- library(BRugs) -->
+<!-- # Step 1 check model -->
+<!-- modelCheck(paste(bugpath, "/backupfiles/gambia-model.txt", sep = ""))  -->
+<!-- # Load the data  -->
+<!-- modelData(paste(bugpath, "/backupfiles/gambia-data.txt", sep = ""))      -->
+<!-- # compile the model with two separate chains -->
+<!-- modelCompile(numChains = 2)  -->
+<!-- # generate initial values  -->
+<!-- # the choice is arbitrary -->
+<!-- initlist <- list(alpha = 0, beta = 1, gamma = 5, logsigma2 = 1) -->
+<!-- modelInits(bugsData(initlist)) -->
+<!-- initlist1 <- list(alpha = 10, beta = 0, gamma = -5, logsigma2 = 5) -->
+<!-- modelInits(bugsData(initlist1)) -->
+
+<!-- # Set monitors on nodes of interest#### SPECIFY, WHICH PARAMETERS TO TRACE: -->
+<!-- parameters <- c("alpha", "beta", "gamma", "sigma2") -->
+<!-- samplesSet(parameters) -->
+
+<!-- # Generate 1000 iterations -->
+<!-- modelUpdate(1000) -->
+<!-- ``` -->
+
+
+```{r BayesianChapter0505, cache=FALSE, fig.width=7, fig.height=8, fig.cap='History plots for iterations 1-1000 for the Gambia example.', fig.align='center', out.width='80%', message=FALSE, warning=FALSE, echo=FALSE}
 library(BRugs)
 # Step 1 check model
 modelCheck(paste(bugpath, "/backupfiles/gambia-model.txt", sep = "")) 
@@ -3237,10 +3265,6 @@ samplesSet(parameters)
 
 # Generate 1000 iterations
 modelUpdate(1000)
-```
-
-
-```{r BayesianChapter0505, cache=FALSE, fig.width=7, fig.height=8, fig.cap='History plots for iterations 1-1000 for the Gambia example.', fig.align='center', out.width='80%', message=FALSE, warning=FALSE, echo=FALSE}
 #### PUT THE SAMPLED VALUES IN ONE R DATA FRAME:
 chain <- data.frame(alpha  = samplesSample("alpha"),
                     beta  = samplesSample("beta"),

backupfiles/logistic-reg-model-centred-stat.txt

@@ -7,8 +7,10 @@ model{
   }
   # priors 
   beta0 ~ dunif(-100, 100)
-  beta1 ~ dunif(-100, 100)
+  beta1 ~ dnorm(1, 0.000000001)
+
+  # generated values
   OR <- exp(beta1) # odds ratio of positive response per 1 mg increase in dose
-  ED95 <- (logit(0.95) - beta0)/beta1 + mean(x[]) # dose that gives 95% of maximal response
+  ED95 <- mean(x[]) + ((logit(0.95) - beta0) / (beta1))  # dose that gives 95% of maximal response
   logit(P35) <- beta0 + beta1 * (35 - mean(x[]))
-}
+}