Введите код в функцию
apply_fun <- function() {
x <- c(11, 5, 2, -5, 7, 2, -11, 9, -5, -5, -4, 17, 2, -10, -11, -10,-4, 2, 1, 13)
a <- 0.1
b <- 0.1
c <- 0
d <- 100^2
M <- 1e3
sample <- array(NA, dim=c(M,2))
mu <- mean(x)
sig2 <- var(x)
for( m in 1:M ){
mu <- rnorm(1, (length(x) + 1/d)^(-1) * (sum(x) + c/d),
sqrt( sig2/(length(x) + 1/d) ))
sig2 <- rigamma(1, .5*length(x)+a+.5,
.5*sum( (x-mu)^2 ) + 1/(2*d)*(mu-c)^2 + b )
sample[m,] <- c(mu,sig2)
}
plot( density( sample[,1] ))
plot( density( sample[,2] ))
}
и затем используйте lapply
для каждого начального значения
output <- lapply(c(314159, 523626, 626789), function(x) {set.seed(x);apply_fun()})
, где rigamma
равно
rigamma = function(n, a, b) return(1/rgamma(n, shape = a, rate = b))