正規線形モデルにおけるパラメータ（）に対する共役事前分布は，正規逆カイ2乗分布（normal inverse chi-square distribution）であることが知られています（Hagan and Forster 2004: pp.305-307, 繁桝 1985: pp.176-177）．
正規逆カイ2乗分布の乱数を生成する方法はいくつかあるようですが（ギブスサンプリングで生成するのが王道かもしれません？）（←削除：McAlinn先生からTwitterにて教えていただいた（2021年1月3日）のですが，この記事で述べる方法も「ギブスサンプリング」と呼べるとのことです．と
の条件付き分布を交互に生成するものだけを「ギブスサンプリング」と呼ぶと私は勘違いしていました．）
，
と分解して，前者の周辺分布が逆カイ2乗分布，後者の条件付き分布が多変量正規分布に従うことを利用するのが，素直な方法のようにも思えます．
以下のRコードは，その素直な方法で正規逆カイ2乗分布の乱数を生成しています．また，の1変量周辺分布のヒストグラムを描いて，それが解析的な答え（位置-尺度型のt分布）になっていることを見た目で確認しています．
（の分散共分散行列が
となるようにパラメタライゼーションしています．また，
です．）
set.seed(20220103)
b0 <- c(2, 9, 1)
S0 <- matrix(c(3,1,1,1,5,1,1,1,7), nrow=3, ncol=3)
l0 <- 7.5;
nu0 <- 10;
nsim <- 100000;Ch0 = chol(S0);
p <- length(b0)result <- matrix(nrow = nsim, ncol = 4);
for(i in 1:nsim){
tau <- rchisq(1,nu0) /(nu0 * l0);
b <- b0 + (Ch0/sqrt(tau)) %*% rnorm(p);
result[i, 1] <- tau;
result[i, 2:(p+1)] <- t(b);
};scale <- matrix(nrow = p, ncol = 1)
par(mfrow = c(2,2))
for(i in 1:p){
scale[i] <- sqrt(S0[i, i])*sqrt(l0);
hist(result[,i + 1], prob = TRUE)
curve((1/scale[i])*dt((x -b0[i])/scale[i], nu0), add = TRUE)
}
上記のコードは事前分布の乱数を生成するものです．事後分布は，ハイパーパラメータがちょっと複雑なものになると思います（試していません）．
（...そもそも論として，乱数を生成するのであれば，共役事前分布を仮定する実用上の意義がないように思います...．）
参考文献
O'Hagan, A. and Forster, J. (1994:1st ed., 2004:2nd ed.) Kendall's Advanced Theory of Statistics. Vol.2B. Bayesian Inference.
