标签: bayesian

我一直在研究一个三级 RStan 模型，其中重复的宽带测量（年份 ID = yrid）嵌套在地方当局（LA ID = 铺设）内，最终嵌套在区域内（区域 ID = rnid）。（记录的）因变量是速度，（记录的）预测变量是人口密度 (pd) 和超高速宽带渗透率 (sfbb)。目前在地方当局和区域级别（2 级和 3 级）有随机拦截。

如何扩展模型以在 1 级或 2 级具有随机斜率？

这是数据的子集、RStan 模型和整个 R 代码。任何帮助将非常感激。

 library(rstan)
 ###Data
  yrid = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
     21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,
     41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,
     61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,
     81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,
     101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,
     121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,
     141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,
     161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,
     181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199)
  laid <- c(1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,
     6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10,
     11,11,11,11,12,12,12,12,13,13,13,13,14,14,14,14,15,15,15,15,
     16,16,16,16,17,17,17,17,18,18,18,18,19,19,19,19,20,20,20,20,
     21,21,21,21,22,22,22,22,23,23,23,23,24,24,24,24,25,25,25,25,
     26,26,26,26,27,27,27,27,28,28,28,28,29,29,29,29,30,30,30,30,
     31,31,31,31,32,32,32,32,33,33,33,33,34,34,34,34,35,35,35,35,
     36,36,36,36,37,37,37,37,38,38,38,38,39,39,39,39,40,40,40,40,
     41,41,41,41,42,42,42,42,43,43,43,43,44,44,44,44,45,45,45,45,
     46,46,46,46,47,47,47,47,48,48,48,48,49,49,49,49,50,50,50)
 rnid <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,
     3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,
     4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
     4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
     5,5,5,5,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
     6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,
     7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,8,8)
pd <- c(7.59262,7.59875,7.6027,7.60375,7.5301,7.53444,7.53604,7.54136,8.378,8.3936,
   8.40061,8.41183,7.36682,7.36992,7.37607,7.38268,7.20065,7.2011,7.20162,7.20578,
   7.78846,7.79947,7.80743,7.81992,7.71797,7.72011,7.72396,7.73026,7.66336,7.66561,
   7.66744,7.66833,7.66973,7.67587,7.68327,7.69321,7.4326,7.43449,7.43762,7.44167,
   7.43053,7.43053,7.43189,7.43396,8.33459,8.34315,8.34548,8.35036,7.15921,7.16325,
   7.16379,7.16943,7.4898,7.48869,7.48689,7.48796,7.61918,7.62046,7.62075,7.62261,
   6.55763,6.56541,6.57438,6.58286,6.27777,6.27833,6.28133,6.28339,6.80184,6.8045,
   6.80572,6.81113,7.31315,7.32324,7.32804,7.33446,7.24893,7.24843,7.24744,7.24993,
   7.80751,7.81927,7.83475,7.84514,7.80045,7.80147,7.80543,7.80792,7.74119,7.74253,
   7.74323,7.74457,7.6027,7.6042,7.60564,7.60852,8.29695,8.30721,8.31356,8.32186,
   8.07527,8.09465,8.11516,8.13795,8.06994,8.07091,8.07347,8.07788,8.19141,8.19883,
   8.20841,8.21603,7.05652,7.05893,7.06613,7.07089,7.85991,7.86511,7.8699,7.87721,
   8.18894,8.19332,8.19572,8.20125,7.26382,7.26669,7.2701,7.27351,6.32972,6.33505, …

我是 R 中 MCMCglmm 包的新手，一般来说，我对 glm 模型还是很陌生。我有一个物种特征数据集，以及它们是否被引入到它们的原生范围之外。

我想测试是否可以通过任何物种特征来解释被引入（作为二进制 0/1 响应变量）。我还想纠正物种之间的系统发育。

有人告诉我，对于二元响应，我可以使用 family =“threshold”，我应该将残差方差固定为 1。但是我在处理之前所需的其他参数时遇到了一些问题。

我已经为随机效应指定了 R 值，但是如果我指定 RI 还必须指定 G 并且我不清楚如何决定这个参数的值。我试过设置默认值，但收到错误消息：

Error in MCMCglmm(fixed, random = ~species, data = data2, family = "threshold",  : 
prior$G has the wrong number of structures

我已经阅读了帮助小插曲和课程，但还没有找到一个二进制响应的例子，我不清楚如何决定先验值。这是我到目前为止：

fixed=Intro_binary ~ Trait1+ Trait2 + Trait3 
Ainv=inverseA(redTree1)$Ainv

binary_model = MCMCglmm(fixed, random=~species, data = data, family = "threshold", ginverse=list(species=Ainv),
 prior = list( 
    G = list(),    #not sure about the parameters for random effects.
    R = list(V = 1, fix …

我一直在关注 PyMC3 的高斯混合模型示例：https : //github.com/pymc-devs/pymc3/blob/master/pymc3/examples/gaussian_mixture_model.ipynb 并让它与人工数据集很好地工作。

我已经用一个真实的数据集尝试过它，我正在努力让它给出合理的结果：

关于我应该缩小/扩大/改变哪些参数以获得更好的拟合的任何想法？痕迹似乎很稳定。这是我从示例中调整的模型片段：

model = pm.Model()
with model:
    # cluster sizes
    a = pm.constant(np.array([1., 1., 1.]))
    p = pm.Dirichlet('p', a=a, shape=k)
    # ensure all clusters have some points
    p_min_potential = pm.Potential('p_min_potential', tt.switch(tt.min(p) < .1, -np.inf, 0))


    # cluster centers
    means = pm.Normal('means', mu=[0, 1.5, 3], sd=1, shape=k)
    # break symmetry
    order_means_potential = pm.Potential('order_means_potential',
                                     tt.switch(means[1]-means[0] < 0, -np.inf, 0)
                                     + tt.switch(means[2]-means[1] < 0, -np.inf, 0))

    # measurement error
    sd = pm.Uniform('sd', lower=0, upper=2, shape=k)

    # …

我在使用基于 Metropolis-Hastings 算法的 MCMC 技术对后验分布进行采样（因此是贝叶斯方法）方面很新。为此，我在 R 中使用了 mcmc 库。我的分布是多维的。为了检查这个metro算法是否适用于多元分布，我在一个多维student-t分布（包mvtnorm，函数dmvt）上成功地做到了。现在我想对我的多元分布（2 个变量 x 和 y）应用同样的东西，但它不起作用；我收到一个错误：X[, 1] 中的错误：维数不正确

这是我的代码：

library(mcmc)
library(mvtnorm)
my.seed <- 123

logprior<-function(X,...)
{
      ifelse( (-50.0 <= X[,1] & X[,1]<=50.0) & (-50.0 <= X[,2] & X[,2]<=50.0), return(0), return(-Inf))
}

logpost<-function(X,...)
{
      log.like <- log( exp(-((X[,1]^2 + X[,2]^2 - 4)/10 )^2) * sin(4*atan(X[,2]/X[,1])) )
      log.prior<-logprior(X)
      log.post<-log.like + log.prior # if flat prior, the posterior distribution is the likelihood one
      return (log.post)
}

x <- seq(-5,5,0.15)
y <- seq(-5,5,0.15)
X<-cbind(x,y)

#out <- metrop(function(X) …

runjagsR的软件包非常棒。并行能力和使用extend.jags功能的能力让我的生活变得更好。然而，有时，在我运行模型后，我意识到老化阶段应该更长。如何从我的run.jags输出中修剪额外的样本，以便我可以重新估计我的参数分布并检查收敛？

jags.object <- run.jags(model, n.chains=3, data=data, monitor =c('a','b'), sample=10000)

是否有可能在 rstanarm 中提取用于 MCMC 采样的 stan 代码？

我想将我自己的模型参数化和先前选择与 rstanarm 中使用的参数化进行比较。

文献称MCMC中的metropolis-hasting算法是上世纪发展起来的最重要的算法之一，具有革命性。文献还说，正是 MCMC 的这种发展给了贝叶斯统计第二次诞生。

我了解 MCMC 的作用 - 它提供了一种从任何复杂概率分布中抽取样本的有效方法。

我也知道贝叶斯推理是什么——它是计算参数的完整后验分布的过程。

我很难在这里连接点：MCMC 在贝叶斯推理过程中的哪一步发挥作用？为什么 MCMC 如此重要以至于人们说是 MCMC 给了贝叶斯统计第二次出生？

想知道你们中是否有人知道为什么 JAGS 会告诉我这里的初始值与尺寸不匹配。

我正在尝试拟合一个空间明确的捕获-重新捕获模型，在该模型中我在每个时间步估计鱼的位置 (x,y)。对于 T=21 时间步长，有 M=64 个人。这是在 array 中估计的s，它通过 i=M 和 t=T 从每个坐标的两个正态分布中循环——x,y。使维度s= (64,2,21)。

我对这个数组的初始值是合适栖息地内的合理位置，是一个维度为 64、2、21 的数组。

然而，JAGS 给了我错误，Error in setParameters(init.values[[i]], i) : RUNTIME ERROR: Dimension mismatch in values supplied for s. 如果我只是不初始化它，我会得到同样的错误，但对于z维度为 64,21的状态矩阵。如果我也不提供 的值z，我会得到错误Error in node y[1,1,2] Node inconsistent with parents，y我的观察数组在哪里，维度 64、7、21（第二个元素是 #of 检测门）。

非常感谢任何帮助。完整代码见下文。

在SourceForge交叉发布，并附有初始值s。不太确定如何在这里发布它。

sink("mod.txt")
cat("
model{


#######################################
# lamda0 = baseline detection rate, gate dependent and …

我正在尝试将 PYMC3 用于贝叶斯模型，我想在新的看不见的数据上反复训练我的模型。我想我每次看到数据时都需要用先前训练过的模型的后验更新先验，类似于这里https://docs.pymc.io/notebooks/updating_priors.html 的实现方式。他们使用以下函数从样本中找到 KDE，并使用对 from_posterior 的调用替换模型中参数的每个原始定义。

def from_posterior(param, samples):
    smin, smax = np.min(samples), np.max(samples)
    width = smax - smin
    x = np.linspace(smin, smax, 100)
    y = stats.gaussian_kde(samples)(x)

    # what was never sampled should have a small probability but not 0,
    # so we'll extend the domain and use linear approximation of density on it
    x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
    y = np.concatenate([[0], y, [0]])
    return Interpolated(param, x, y)

这是我的原始模型。 …

我正在使用 gemtc 包对类似于以下的数据集进行贝叶斯网络元分析：

df <- data.frame(study = c("A", "A", "B", "B", "C", "C", "D", "D", "E", "E", "F", "F", 
                           "G", "G", "H", "H", "I", "I", "J", "J", "K", "K"),
                                 treatment = c("A", "B", "B", "C", "B", "C", "A", "B", "B", 
                                                "C", "B", "C", "A", "B", "B", "C", "B", "C", 
                                                "A", "C", "B", "C"),
                                 responders = c(1, 5, 0, 0, 3, 1, 0, 2, 0, 2, 0, 2, 0,
                                                0, 1, 2, 0, 0, 2, 9, 1, 1),
                                 sampleSize = c(32, …