标签: bayesian

我正在尝试研究基于"做贝叶斯数据分析:R,JAGS和斯坦(2015)的教程"一书中的贝叶斯分析.

在本书中,有一些例子.所以,我试图在R中复制这个例子.但是,在这个例子中我收到了一条错误信息.

具体而言,这是示例数据.

data
   y        s
1  1 Reginald
2  0 Reginald
3  1 Reginald
4  1 Reginald
5  1 Reginald
6  1 Reginald
7  1 Reginald
8  0 Reginald
9  0     Tony
10 0     Tony
11 1     Tony
12 0     Tony
13 0     Tony
14 1     Tony
15 0     Tony

y<-data$y
s<-as.numeric(data$s)
Ntotal=length(y)
Nsubj=length(unique(s))

dataList=list(y=y, s=s, Ntotal=Ntotal, Nsubj=Nsubj)

另外,这是我的模特.

modelString=" 
model{
  for(i in 1:Ntotal){
    y[i] ~ dbern(theta[s[i]])
  }
  for(s in 1:Nsubj){
    theta[s] ~ dbeta(2,2)
  }
}
"
writeLines(modelString, con="TEMPmodel.txt") …

我正在寻找一个R可用于根据计数数据训练狄利克雷先验的包。我正在询问一位正在使用 的同事R，而我自己并没有使用它，所以我不太清楚如何寻找软件包。搜索起来有点困难，因为“R”是一个非特定的搜索字符串。CRAN上好像没有什么，但是还有其他地方可以看吗？

我正在使用 OpenBUGS 和 R 包研究二项式混合模型R2OpenBUGS。我已经成功构建了更简单的模型，但是一旦我添加另一个级别来进行不完善的检测，我就会不断收到错误variable X is not defined in model or in data set。我尝试了许多不同的方法，包括更改数据结构以及将数据直接输入 OpenBUGS。我发布此内容是希望其他人有此错误的经验，并且也许知道为什么 OpenBUGS 无法识别变量 X，尽管据我所知它已明确定义。

我也收到了错误expected the collection operator c error pos 8- 这不是我以前收到的错误，但我也同样感到困惑。

模型和数据模拟函数均来自 Kery 的《生态学家 WinBUGS 简介》(2010)。我要注意的是，这里的数据集是代替我自己的数据的，这是相似的。

我包括构建数据集和模型的函数。抱歉长度。

# Simulate data: 200 sites, 3 sampling rounds, 3 factors of the level 'trt', 
# and continuous covariate 'X'

data.fn <- function(nsite = 180, nrep = 3, xmin = -1, xmax = 1, alpha.vec = c(0.01,0.2,0.4,1.1,0.01,0.2), beta0 = 1, beta1 …

我正在使用贝叶斯推理做课堂作业。为此，我使用MCMCregress来自 的函数MCMCpack。当我想要获得残差时，问题就出现了，因为该函数不提供残差，所以我必须“手动”计算它们（在 R 中）。

我的模型是：

Model <- MCMCregress(Y~X1+X2+X3+X4+X5, data=DATA)

其中X1和X5是连续的，而X2和X3是X4二分的。模型输出为我提供了每个变量的估计：

(Intercept) = 1.90, 
X1 = -0.02, 
X2 = -0.05, 
X3 = 0.32,
X4 = 0.61,
X5 = -0.003,
sigma2 = 0.54

我想我必须这样做：

1.90 - 0.02*X1 - 0.05*X2 + 0.32*X3 ...

但我知道我在 R 代码中遗漏了一些重要的东西，所以我想知道哪个是正确的 R 代码以获得残差。

这是一个可重现的示例（尽管它与原始数据不对应）：

Y  <- c(0.2,0.8,1,4.3,5,3.5,3.2,1,3.3,1,2,4,3.6,3,5,4.3,3.2,4,1,2)
X1 <- c(17,13,15,NA,12,24,15,NA,12,22,14,12,18,NA,10,13,12,11,26,10)
X2 <- c(0,0,0,1,0,0,1,1,1,1,0,1,0,0,0,0,0,1,NA,NA)
X3 <- c(0,0,1,1,1,0,1,0,0,0,0,1,0,NA,0,NA,NA,0,1,0)
X4 <- c(1,0,1,0,0,1,0,0,NA,0,NA,NA,0,0,1,0,1,1,0,1)
X5 <- c(2.46,4.56,32.1,NA,NA,NA,NA,NA,NA,3.76,5.67,4.56,NA,
        17.32,12.2,4.56,7.2,1.2,NA,9.2)

X2 …

出于比较目的，我想利用 PyMC3 之外的后验密度函数。

对于我的研究项目，我想了解 PyMC3 与我自己定制的代码相比的性能如何。因此，我需要将它与我们自己的内部采样器和似然函数进行比较。

我想我想出了如何调用内部PyMC3后验的方法，但是感觉很别扭，想知道是否有更好的方法。现在我正在手动转换变量，而我应该能够将参数字典传递给 pymc 并获得后验密度。这可能以直接的方式进行吗？

非常感谢！

演示代码：

import numpy as np
import pymc3 as pm
import scipy.stats as st

# Simple data, with sigma = 4. We want to estimate sigma
sigma_inject = 4.0
data = np.random.randn(10) * sigma_inject

# Prior interval for sigma
a, b = 0.0, 20.0

# Build PyMC model
with pm.Model() as model:
    sigma = pm.Uniform('sigma', a, b)      # Prior uniform between 0.0 and 20.0
    likelihood = pm.Normal('data', 0.0, sd=sigma, observed=data)

# …

我一直在关注 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)

    # …

我使用以下行来更新每次试验中的 beta 发行版并给出 arm 推荐（我使用 scipy.stats.beta）：

self.prior = (1.0,1.0)
def get_recommendation(self):
    sampled_theta = []
    for i in range(self.arms):
        #Construct beta distribution for posterior
        dist = beta(self.prior[0]+self.successes[i],
                    self.prior[1]+self.trials[i]-self.successes[i])
        #Draw sample from beta distribution
        sampled_theta += [ dist.rvs() ]
    # Return the index of the sample with the largest value
    return sampled_theta.index( max(sampled_theta) )

但目前，它只适用于奖励是二元的（成功或失败）。我想修改它，使其适用于非二元奖励。（例如奖励：2300、2000，...）。我怎么做？

背景：

通常，R 给出众所周知的分布的分位数。在这些分位数中，较低的 2.5% 到较高的 97.5% 覆盖了这些分布下 95% 的面积。

问题：

假设我有一个 F 分布（df1 = 10，df2 = 90）。在 R 中，如何确定此分布下的 95% 区域，使得该 95% 仅覆盖高密度区域，而不是 R 通常给出的 95%（请参阅下面的 R 代码）？

注意：显然，最高密度是“模式”（下图中的虚线）。所以我想，人们必须从“模式”转向尾部。

这是我的 R 代码：

curve(df(x, 10, 90), 0, 3, ylab = 'Density', xlab = 'F value', lwd = 3)

Mode = ( (10 - 2) / 10 ) * ( 90 / (90 + 2) )

abline(v = Mode, lty = 2)

CI = qf( c(.025, .975), 10, …

我正在努力实现一个模型，其中狄利克雷变量的集中因子依赖于另一个变量。

情况如下：

系统由于组件故障而失败（共有三个组件，每次测试/观察时只有一个组件失败）。

组件发生故障的概率取决于温度。

这是该情况的（已注释的）简短实现：

import numpy as np
import pymc3 as pm
import theano.tensor as tt


# Temperature data : 3 cold temperatures and 3 warm temperatures
T_data = np.array([10, 12, 14, 80, 90, 95])

# Data of failures of 3 components : [0,0,1] means component 3 failed
F_data = np.array([[0, 0, 1], \
       [0, 0, 1], \
       [0, 0, 1], \
       [1, 0, 0], \
       [1, 0, 0], \
       [1, 0, 0]])

n_component = 3

# When temperature …

我发现到处都发布了同样的问题，而且我似乎无法找到任何适合我的数据的解决方案，我想知道我是否正在尝试将我的数据拟合到一个模型太复杂了。

我正在尝试将我的数据拟合到 MCMCglmm 包中的多项逻辑回归模型。我查看了许多不同的文档、教程和 MCMCglmm 手册本身，主要是Florian Jaeger 的教程，它非常出色且全面。然而，我在他选择先验的 G 结构和 R 结构的值时迷失了方向，并且不断收到此错误消息

Error in priorformat(if (NOpriorG) { : \n  V is the wrong dimension for some prior$G/prior$R elements\n

特别是，我不确定n在这两个数据中应该给我的数据赋予什么值，但有些不透明的错误消息表明这是一个问题V

这是我的数据的（相关）子集：

CG_imm       locuteur   enquete       loc_age   loc_sexe      left    liquid   right    articulation_C1    voice_C1   NC_C1    NC_right    voice_right    right2          pos    logfreq      realization\nabordable    44ajs1     Nantes        79         M            bl      l        p        stop               V          Vstop    stop        NV             NVstop          adj    NA           2\nadmettre     91adb1     Brunoy        54         M            tR      R        E …