bayesian hw 3.qmd

# bayesian hw 3 {.unnumbered}

---
title: Bayesian hw 3
subtitle: chapter 3
author: Congyuan He
date: now

format: 
  html:
    toc: true
    toc-location: left
    toc-depth: 3
    toc-expand: 1
    number-depth: 3
---

## 第一题

![](images/clipboard-1853463496.png)

## a

临床试验的例子，

安慰剂服从$Y～N(\mu,\sigma^2) ,i=1,..,n$ ,

实验药服从$N～N(\mu + \delta,\sigma^2),i=1+n,..,n+m$

## b

![](images/clipboard-420108849.png)

## c

![](images/clipboard-1203345143.png)

```{r}
# 先验参数
n <- m <- 50
mu <- 10
delta <- 1
sigma <- 2

# 模拟数据
y1 <- rnorm(n,mu,sigma)
y2 <- rnorm(m,mu+delta,sigma)
y <- c(y1,y2)

# Gibbs
set.seed(2024)
S<-10000
PHI<-matrix(nrow=S,ncol=3)
PHI[1,]<-phi<-c(0, 0, 4) # sigma2

## Gibbs sampling algorithm
for(s in 2:S) {

# generate a new mu value from its full conditional
mu_n<- n*100*100*mean(y1)/(4+n*100*100)
s2_n<- (100*100*4)/(4+n*100*100)
phi[1]<-rnorm(1, mu_n, sqrt(s2_n) )

# generate a new delta value from its full conditional
# delta_m<-  ( m*100*100*mean(y2)/(4+m*100*100))-( n*100*100*mean(y1)/(4+n*100*100))
delta_m<-  (100*100*sum(y2-mu_n)/(4+m*100*100))

s2_m<- (100*100*4)/(4+m*100*100)
phi[2]<-rnorm(1, delta_m, sqrt(s2_m) )

# generate a new sigma^2 value from its full conditional
a <- 0.01+(m+n)*0.5
b <-0.01+(m+n)*0.5*((sum((y-mean(y))^2)/(m+n))+((mean(y)-(n*mu+m*(mu+delta))/(m+n))^2))
phi[3]<- 1/rgamma(1, a, b)
PHI[s,]<-phi         }


```

```{r}
plot(density(PHI[,1]),main='mu')
plot(density(PHI[,2]),main='delta')
plot(density(PHI[,3]),main='sigma2')
```

## d

```{r}
library('rjags')
library("coda")


model_code <- "
model {

  for (i in 1:n) {
    Y[i] ~ dnorm(mu, 1/sigma2)
  }
  
  for (i in (n+1):(n+m)) {
    Y[i] ~ dnorm(mu + delta, 1/sigma2)
  }

  mu ~ dnorm(0, 1/10000)
  delta ~ dnorm(0, 1/10000)
  tau ~ dgamma(0.01, 0.01)
  sigma2 =  1/tau
}
"

n <- 50
m <- 50
mu_n <- 10
delta <- 1
sigma <- 2

set.seed(2024)
Y <- c(rnorm(n, mu_n, sigma), rnorm(m, mu_n + delta, sigma))

data_list <- list(Y = Y, n = n, m = m)



model <- jags.model(textConnection(model_code), data = data_list, n.chains = 3, n.adapt = 5000)

# 运行 MCMC 采样
update(model, 5000) # Burn-in阶段
samples <- coda.samples(model, variable.names = c("mu", "delta", "sigma2"), n.iter = 5000)

plot(samples)

```

## e

```{r}
# gelman.diag(samples)
```

```{r}
# 自相关
acf(as.matrix(samples)[, "mu"])
acf(as.matrix(samples)[, "delta"])
acf(as.matrix(samples)[, "sigma2"])
```

```{r}
# 有效样本量
effectiveSize(samples)
```

## 第二题

![](images/clipboard-1198587299.png)

## a
首先共轭
beta分布对于概率$\theta$来说是常规的先验分布设置手段。因为其足够灵活并且在0，1之间变化; 其次，这样设置还有现实意义的解释，即beta分布中，a=命中多少个，b=未能命中的个数; 最后，添加了一个exp(m)这样一个调节的量，表示在常规命中率的基础上引入变化
exp（m）等于试验总次数取exp保持其为正，exp(m)q/(a+b)=成功的概率刚好等于q，表示平时命中率用到了样本信息

```{r}

set.seed(2024)
theta_i <- c(0.845, 0.847, 0.880, 0.674, 0.909, 0.898, 0.770, 0.801, 0.802, 0.875)
m_samples <- rnorm(1000, mean = 0, sd = 10)

probabilities <- numeric(10000)

for (i in 1:1000) {
  m <- m_samples[i] 

  for (j in 1:10) {
    shape1 <- exp(m) * theta_i[j]       
    shape2 <- exp(m) * (1 - theta_i[j]) 
    probability <- rbeta(1, shape1, shape2)
    probabilities[(i - 1) * 10 + j] <- probability
  }
}


hist(probabilities,probability = TRUE)
lines(density(probabilities), col = "red", lwd = 2) 

```

## b

m引入了不确定性
exp(m)先验样本量

## c

后验：

$$
beta(y+exp(m)*q_{i} , exp(m)*(1-q_i)+n-y)
$$

## d

```{r}
### data input
q <- c(0.845,0.847,0.880,0.674,0.909,0.898,0.770,0.801,0.801,0.875)
y <- c(64,72,55,27,75,24,28,66,40,13)
n <- c(75,95,63,39,83,26,41,82,54,16)

### prior
mu_m <- 0
tau_m_2 <- 10

### Metropolis + Gibbs sampling
set.seed(1)
S<-10000
PHI<-matrix(nrow=S,ncol=length(n)+1)
PHI[1,]<-phi<-c(q, 0)

# help function for Metropolis step
log_m <- function(theta,m){
  like  <- sum(dbeta(theta,exp(m)*q,exp(m)*(1-q),log=TRUE))
  prior <- dnorm(m,mu_m,sqrt(tau_m_2),log=TRUE)
  return(like+prior)}

# proposal candidate standard deviation
can_sd <- 2

# monitoring acceptance rate
acs <- 0

for(s in 2:S){
  
  # Metropolis for m
  can    <- rnorm(1,phi[11],can_sd)
  logR   <- log_m(phi[1:10],can)-
    log_m(phi[1:10],phi[11])
  if(log(runif(1))<logR){
    phi[11] <- can; acs <- acs+1 
  }
  
  # Gibbs for theta
  for (i in 1:length(n) ) {
    phi[i] <- rbeta(1, exp(phi[11])*q[i]+y[i],exp(phi[11])*(1-q[i])+n[i]-y[i] )
  }
  
  PHI[s,] <- phi
}

acs/S

# 目标接受率
# 
# •	0.2 到 0.5：
# 这是 Metropolis 算法的常见推荐范围。一般来说，较低的接受率可能意味着步长太大，
# 导致大量候选值被拒绝；较高的接受率可能意味着步长太小，导致缺乏足够的探索。适中的接受率通常是最理想的。
# 
# 如何调整接受率？
# 
# •	调整候选分布的标准差（can_sd）：
# 如果接受率较低，可以尝试增大 can_sd，从而让候选值更有可能接受。
# 如果接受率较高，可以尝试减小 can_sd，以避免候选值过于接近当前值，从而增加采样的多样性。

# 在Metropolis-Hastings算法中，can_sd 是用来生成候选值（candidate values）的标准差，它决定了每次提议的新候选点（can）与当前点（phi[11]）之间的距离。
# 这是一个非常重要的超参数，它影响了算法的收敛速度和效率。具体来说，它决定了提议的新参数值的“步长”，即每次更新时的幅度

apply(PHI[5000:10000, ], 2, 
      function(x) quantile(x,c(0.05,0.5,0.95)))


```

## e

```{r}
library(rjags)
library(coda)

model_string <- "
model {


  for (i in 1:n_players) {

    clutch_makes[i] ~ dbin(theta[i], clutch_attempts[i])

    theta[i] ~ dbeta(alpha[i], beta[i])

    alpha[i] <- exp(m) * theta_i[i]
    beta[i] <- exp(m) * (1 - theta_i[i]) 
  }
   m ~ dnorm(0,1/10)
}
"

theta_i <- c(0.845, 0.847, 0.880, 0.674, 0.909, 0.898, 0.770, 0.801, 0.802, 0.875)
clutch_makes <- c(64, 72, 55, 27, 75, 24, 28, 66, 40, 13)
clutch_attempts <- c(75, 95, 63, 39, 83, 26, 41, 82, 54, 16)
n_players <- length(theta_i)


data_jags <- list(
  theta_i = theta_i,         
  clutch_makes = clutch_makes,
  clutch_attempts = clutch_attempts,
  n_players = n_players
)



model <- jags.model(textConnection(model_string), data = data_jags)

update(model, 10000)

n_samples <- 10000
samples <- coda.samples(model, variable.names = c('theta','m'), n.iter = n_samples)

# plot(samples)
summary(samples)



```

## f

自己的代码需要推导后验分布或者用更复杂的M-H方法

## g


```{r}
effectiveSize(samples)
```

## 第三题

![](images/clipboard-1294109230.png)

## a

```{r}
library(MASS)
data(galaxies)
?galaxies
Y <- galaxies

n <- length(Y)

hist(Y,breaks=25)


```

$\mu=20000$

$\sigma=5000$

$k=15$

## b

```{r}
library(rjags)
library(coda)

model_string <- "
model {
  for (i in 1:s) {
    Y[i] ~ dt(mu, tau, k)
  }
  mu ~ dnorm(0, 1/10000^2)
  tau ~ dgamma(0.01, 0.01)
  k ~ dunif(1, 30)
}
"

data_jags <- list(Y =galaxies, s = length(galaxies))

# 设置初始值
init_values <- list(
  list(mu = 0, tau = 100, k = 15)
)

# 创建模型
model <- jags.model(textConnection(model_string), data = data_jags, inits = init_values, n.chains = 1)

# 更新模型
update(model, 10000)

# 采样
n_samples <- 10000
samples <- coda.samples(model, variable.names = c("mu", "tau", "k"), n.iter = n_samples)

# 绘制轨迹图
plot(samples)
```

## c

```{r}

posterior_means <- apply(as.matrix(samples), 2, mean)

# 从后验均值计算t分布
mu_posterior <- posterior_means["mu"]
sigma_posterior <- 1/sqrt(posterior_means["tau"])
k_posterior <- posterior_means["k"]

# 生成t分布的随机样本
t_dist_samples <- rt(10000, df = k_posterior, ncp = mu_posterior)

# 绘制观察数据和t分布的对比图
hist(galaxies, probability = TRUE, col = rgb(0, 0, 1, 0.5), main = "Comparison of Data and t-distribution", xlim = c(min(galaxies), max(galaxies)))
lines(density(t_dist_samples), col = "red", lwd = 2,xlim=c(10000,35000))


```