将参数传递给确定性变量，pymc(Passing parameters to determinist

我试图实现使用PyMC大数定律的一个很简单的例子。我们的目标是产生不同大小的样品的许多示例平均值。例如，在下面的代码，我正在5个样品多次服用组（samples_to_average = 5），计算它们的平均，和然后找到所得轨迹的95％CI。

下面的代码运行，但我想要做的就是修改samples_to_average是一个清单，这样我就可以计算置信区间在一个循环中的一系列不同的样本量。

import scipy.misc
import numpy as np
import pymc as mc

samples_to_average = 5 
list_of_samples = mc.DiscreteUniform("response", lower=1, upper=10, size=1000)

@mc.deterministic
def sample_average(x=list_of_samples, n=samples_to_average):
    samples = int(n)
    selected = x[0:samples] 
    total = np.sum(selected) 
    sample_average = float(total) / samples 
    return sample_average 

def getConfidenceInterval():   
    responseModel = mc.Model([samples_to_average, list_of_samples, sample_average])
    mapRes = mc.MAP(responseModel)
    mapRes.fit() 
    mcmc = mc.MCMC(responseModel)
    mcmc.sample( 10000, 5000)
    upper = np.percentile(mcmc.trace('sample_average')[:],95)
    lower = np.percentile(mcmc.trace('sample_average')[:],5)
    return (lower, upper)     


print getConfidenceInterval()

大多数例子我使用确定性的装饰中使用全局变量随机看到。然而，为了实现我的目标，我想我需要做的是建立在getConfidenceInterval（）随机变量（正确的长度），并通过这sample_average（而不是使用全局/默认参数提供sample_average）。

如何能在getConfidenceInterval创建的变量（）传递到sample_average（），或者是什么，我可以使用samples_to_average的不同值评估多个模型另一种方式？我想如果可能的话，以避免全局变量。

解决你的问题之前，我想简化sample_average的方式是这么写的，这是更紧凑，更容易理解。

sample_average = mc.Lambda('sample_average', lambda x=list_of_samples, n=samples_to_average: np.mean(x[:n]))

现在，你可以概括这个到samples_to_average是一组参数的情况下：

samples_to_average = np.arange(5, 25, 5)

sample_average = mc.Lambda('sample_average', lambda x=list_of_samples, n=samples_to_average: [np.mean(x[:t]) for t in n])

该getConfidenceInterval功能也将不得不改变，如下图所示：

def getConfidenceInterval():
    responseModel = mc.Model([samples_to_average, list_of_samples, sample_average])
    mapRes = mc.MAP(responseModel)
    mapRes.fit()
    mcmc = mc.MCMC(responseModel)
    mcmc.sample( 10000, 5000)
    average = np.vstack((t for t in mcmc.trace('sample_average')))
    upper = np.percentile(average, 95, axis = 0)
    lower = np.percentile(average, 5, axis = 0)
    return (lower, upper)

我用vstack到样品平均聚集成一个二维数组，然后用于在numpy的的百分功能轴线选项计算沿每一列百分。