Binary Probit (rbprobitGibbs) 精度矩阵的 Gibbs 采样器(Albert 和 Chib)

Gibbs Sampler (Albert and Chib) for Binary Probit (rbprobitGibbs) A precision matrix

目前,我正在研究 RStudio 帮助文件中的上述内容,其中包含以下示例:

##
## rbprobitGibbs example
##
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}

set.seed(66)    
simbprobit = function(X,beta) {
    ##  function to simulate from binary probit including x variable
    y=ifelse((X%*%beta+rnorm(nrow(X)))<0,0,1)
    list(X=X,y=y,beta=beta)
}

nobs=200
X=cbind(rep(1,nobs),runif(nobs),runif(nobs))
beta=c(0,1,-1)
nvar=ncol(X)
simout=simbprobit(X,beta)

Data1=list(X=simout$X,y=simout$y)
Mcmc1=list(R=R,keep=1)

out=rbprobitGibbs(Data=Data1,Mcmc=Mcmc1)

summary(out$betadraw,tvalues=beta)

if(0){
    ## plotting example
    plot(out$betadraw,tvalues=beta)
}

当我单步执行代码时,我没有看到任何地方设置了 A 矩阵。只有当我到达这条线时:

out=rbprobitGibbs(Data=Data1,Mcmc=Mcmc1)

我看到输出中显示了 A 矩阵,我知道它必须是 k * k 矩阵,其中 betabar 是 k * 1 矩阵。

Prior Parms: 
betabar
# [1] 0 0 0
A
#      [,1] [,2] [,3]
# [1,] 0.01 0.00 0.00
# [2,] 0.00 0.01 0.00
# [3,] 0.00 0.00 0.01

所以我可以理解A是如何得到它的维度的;但是,我不清楚 A 中的值是如何设置为 0.01 的。我试图弄清楚如何允许用户调用 rbprobitGibbs 函数通过 A 将精度设置为他们喜欢的任何值。我可以看到 A 的输出位置,但它的值是如何基于某些输入的?有没有人有什么建议? TIA。

更新:

这是生成的输出,但据我所知,无论是否使用 prior = list(rep(0,3), .2*diag(3)),它都是相同的:

> out
$betadraw
           [,1]      [,2]       [,3]
 [1,] 0.3565099 0.6369436 -0.9859025
 [2,] 0.4705437 0.7211755 -1.1955608
 [3,] 0.1478930 0.6538157 -0.6989660
 [4,] 0.4118663 0.7910846 -1.3919411
 [5,] 0.0385419 0.9421720 -0.7359932
 [6,] 0.1091359 0.7991905 -0.7731041
 [7,] 0.4072556 0.5183280 -0.7993501
 [8,] 0.3869478 0.8116237 -1.2831395
 [9,] 0.8893555 0.5448905 -1.8526630
[10,] 0.3165972 0.6484716 -0.9857531
attr(,"class")
[1] "bayesm.mat" "mcmc"      
attr(,"mcpar")
[1]  1 10  1

它通过先验精度矩阵的缩放常数得到这个因子。在源代码中,您会注意到,如果您不提供先验精度,那么它将生成一个方 k 矩阵并将其乘以 .1。这里没什么好看的。 bayesm 中所有各种函数的这些缩放参数可以在 ./bayesm/R/bayesmConstants.R file 中找到。

    if (is.null(Prior$A)) {
        A = BayesmConstant.A * diag(nvar)
    }

如果你愿意,你可以提供你自己的常量,比如说 .2,你可以这样做,prior = list(rep(0,k), .2*diag(k)),甚至可以在先验中引入一些相关信息。

参加派对已经很晚了,但我 运行 遇到了同样的问题并且刚刚弄明白了。为了更改 A 矩阵和先验矩阵,您还必须命名它们,因为所有其他输入变量都已命名。

例如你的代码应该是,

rbprobitGibbs(Data=Data1, Prior=list(betabar=betabar1, A=A1), Mcmc=Mcmc1)

如果这样做,您就可以为 betabar 和 A 设置自己的值。