Блог

Если я хочу получить Гиббса проще:

  library(mvtnorm)

 #Pseudo Data
 #Sample Size
 n = 50

 #The response variable
 Y = matrix(rnorm(n,20))

 #The Design matrix
 X = matrix(c(rep(1,n),rnorm(n,3),rnorm(n,10)),nrow=n)
  k = ncol(X)

 #Number of samples
 B = 1100

 #Variables to store the samples in 
  beta = matrix(NA,nrow=B,ncol=k)
  sigma = c(1,rep(NA,B))
  psi = rep(NA,B)

 #The Gibbs Sampler
   for(i in 1:B){

  #The least square estimators of beta and sigma
   V = solve(t(X)%*%X)
   bhat = V%*%t(X)%*%Y

  sigma.hat = t(Y-X%*%bhat)%*%(Y-X%*%bhat)/(n-k)


  #Sample beta from the full conditional 
  beta[i,] = rmvnorm(1,bhat,sigma[i]*V)

  #Sample sigma from the full conditional
  sigma[i 1] = 1/rgamma(1,(n-k)/2,(n-k)*sigma.hat/2)

 #Obtain the marginal posterior of psi
  psi[i] = (beta[i,2] beta[i,3])/sigma[i 1]
 }


 #Plot the traces
   dev.new()
   par(mfrow=c(2,2))
 plot(beta[,1],type='l',ylab=expression(beta[1]),main=expression("Plot of "*beta[1]))
 plot(beta[,2],type='l',ylab=expression(beta[2]),main=expression("Plot of "*beta[2]))
  plot(beta[,3],type='l',ylab=expression(beta[3]),main=expression("Plot of "*beta[2]))
plot(sigma,type='l',ylab=expression(sigma^2),main=expression("Plot of "*sigma^2))
  

Затем я получил график трассировки:

Мой вопрос в том, если я снова запущу, чтобы получить другой график beta_1 (rep. beta_2 ). Как отобразить эти два графика трассировки beta_1 (rep. beta_2 ) на одном и том же рисунке другим цветом? Я не знаю, как написать этот код?

Боюсь, вы не сможете вернуться к графику 1, как только перейдете к графику 4. Итак, сначала вам нужно сохранить ваши симуляции в списке, поэтому оберните ваше моделирование в функцию:

 lambda = function(){

   beta = matrix(NA,nrow=B,ncol=k)
   sigma = c(1,rep(NA,B))
   psi = rep(NA,B)

   for(i in 1:B){

   V = solve(t(X)%*%X)
   bhat = V%*%t(X)%*%Y

   sigma.hat = t(Y-X%*%bhat)%*%(Y-X%*%bhat)/(n-k)
 
   beta[i,] = rmvnorm(1,bhat,sigma[i]*V)
   sigma[i 1] = 1/rgamma(1,(n-k)/2,(n-k)*sigma.hat/2)
   psi[i] = (beta[i,2] beta[i,3])/sigma[i 1]
   }
   return(list(beta=beta,sigma=sigma,psi=psi))
 }
  

Запустите их с помощью replicate :

 results = replicate(2,lambda(),simplify=FALSE)
names(results) = paste0(1:length(results))
  

Теперь определите цвета и постройте:

 cols = c("#F0545480","#30475E80")
names(cols) = c("rep1","rep2")

par(mfrow=c(2,2))
for(i in 1:3){
   plot(results[[1]]$beta[,i],type='l',
   ylab = substitute(expression(beta[x]), list(x = i)),
   main = substitute(expression(beta[x]), list(x = i)),
   col = cols[1])
   lines(results[[2]]$beta[,i],col=cols[2])
}
plot(results[[1]]$sigma,type='l',ylab=expression(sigma^2),main=expression("Plot of "*sigma^2),col=cols[1])
lines(results[[2]]$sigma,type='l',ylab=expression(sigma^2),main=expression("Plot of "*sigma^2),col=cols[2])
  

У меня есть решение для вас, в котором я изменил ваш код, чтобы он работал как автономные функции. my_gibbs_test_data() создает данные, которые подаются в пробоотборник Гиббса my_gsamp() , и обе эти функции сами вызываются, gibbs_eval_test() которые принимают в качестве входных данных общее количество попыток / запусков случайных данных и выборку Гиббса, которая должна быть выполнена (наряду с входными данными для вышеупомянутых функций). Результатом всего этого являются данные из сэмплера в аккуратном формате, которые затем могут быть переданы ggplot() для построения графика.

Я продемонстрировал это ниже. (Я предполагаю, что код достаточно понятен для OP, учитывая, что это просто их собственный код, слегка переупакованный.)

 #Requisite libraries
library(mvtnorm)
library(ggplot2)
library(dplyr)
library(tidyr)
library(magrittr)


#----
#Random data generator for the Gibbs sampler
#Returns the data as a list; [[1]] is X, [[2]] is Y
my_gibbs_test_data <- function(myn = 50, myB = 1100, yot = 20, xot1 = 3, xot2 = 10){
  #Generating test data
  Y <- matrix(rnorm(myn, yot))
  
  X <- matrix(
    c(rep(1, myn), rnorm(myn, xot1), rnorm(myn, xot2)), 
    nrow = myn)
  
  return(list(X, Y))
}
#t2 <- my_gibbs_test_data()
#----


#----
#Gibbs sampler
#This function spits out the data as a data.frame
#Easier for plotting with ggplot() down the line
my_gsamp <- function(myX, myY, myn, myB){
  
  #
  myk <- ncol(myX)
  
  #Least sq. estimators of beta and sigma
  myV <- base::solve(t(myX) %*% myX)
  b_hat <- myV %*% t(myX) %*% myY
  sig_hat <- as.numeric( ( (t(myY - myX %*% b_hat) %*% (myY - myX %*% b_hat)) / (myn - myk) ) )
  
  
  #Sample storage variables
  #Initializing all of these with same column lengths
  #The last rows of mybeta and mypsi will remain NA after
  #the calculations
  mybeta <- matrix(nrow = myB   1, ncol = myk)
  #
  mysigma <- c(1, rep(NA, myB))
  #mysigma <- c(1, rep(NA, myB-1))
  mypsi <- rep(NA, myB   1)
  
  #Main for loop
  for(i in 1:myB){
    #i <- 1
    mybeta[i, ] <- mvtnorm::rmvnorm(1, b_hat, mysigma[i] * myV)
    mysigma[i 1] <- 1/stats::rgamma(1, (myn - myk)/2, (myn - myk) * sig_hat/2)
    mypsi[i] <- (mybeta[i, 2]   mybeta[i, 3]) / mysigma[i 1]
  }
  
  mydat <- data.frame(mybeta, mysigma, mypsi)
  names(mydat) <- c(paste0("b", 1:ncol(mybeta)), "sig", "psi")
  #mydat$index <- seq(nrow(mydat))
  
  return(mydat)
}
#tst <- my_gsamp(t2[[1]], t2[[2]], n, B)
#----


#----
#Function that replicates the Gibbs sampling locreps number of times
gibbs_eval_test <- function(locreps = 4, locn = 50, locB = 1100, locyot = 20, locxot1 = 3, locxot2 = 10){
  
  #All that this function does is calls my_gibbs_test_data() and my_gsamp()
  #as many number of times as specified by the user via locreps
  
  #locn, locB, locyot, locxot1, and locxot2 are the parameters that control the test data
  #and the Gibbs sampler itself
  
  #Note: this function returns the data in the long format
  
  for(i in 1:locreps){
    #i <- 1
    #locreps = 1
    #locn = 50
    #locB = 1100
    #locyot = 20
    #locxot1 = 3
    #locxot2 = 10
    
    loctestdat <- my_gibbs_test_data(locn, locB, locyot, locxot1, locxot2)
    
    if(i == 1){
      outdat <- my_gsamp(loctestdat[[1]], loctestdat[[2]], locn, locB)
      
      #Index for plotting
      outdat$index <- seq_along(outdat$b1)
      #Pivoting longer
      outdat %<>% pivot_longer(cols = -index, names_to = "meas_var", values_to = "meas_val")
      
      outdat$run <- rep(i, nrow(outdat))
      
    } else{
      locdat <- my_gsamp(loctestdat[[1]], loctestdat[[2]], locn, locB)
      
      #Index for plotting
      locdat$index <- seq_along(locdat$b1)
      #Pivoting longer
      locdat %<>% pivot_longer(cols = -index, names_to = "meas_var", values_to = "meas_val")
      
      locdat$run <- rep(i, nrow(locdat))
      
      #Appending to outdat
      outdat <- dplyr::bind_rows(outdat, locdat)
    }
    
  }
  
  return(outdat)
  
}
#----


#----
#Testing
t4 <- gibbs_eval_test()



#Plotting
t4 %>% 
  filter(meas_var != "psi") %>%
  ggplot(aes(color = as.factor(run)))   
  geom_line(aes(y = meas_val, x = index))   
  facet_wrap(~meas_var, scales = "free")

#----
  

Приведенный выше пример графика содержит данные из 4 повторений процесса выборки Гиббса (с использованием встроенных значений по умолчанию). ИМХО, это решение обеспечивает большую гибкость, а также позволяет легко визуализировать.