1

PyMC を使用して大数の法則の非常に単純な例を実装しようとしています。目標は、さまざまなサイズのサンプルの多くのサンプル平均を生成することです。たとえば、以下のコードでは、5 つのサンプル (samples_to_average = 5) のグループを繰り返し取得し、それらの平均を計算してから、結果のトレースの 95% CI を見つけています。

以下のコードは実行されますが、samples_to_average をリストに変更して、1 回のパスでさまざまなサンプル サイズの範囲の信頼区間を計算できるようにしたいと考えています。 

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 に渡す必要があると思います (globals / default パラメーターを使用して sample_average を指定するのではなく)。

getConfidenceInterval() で作成された変数を sample_average() に渡すにはどうすればよいですか? あるいは、samples_to_average の異なる値を使用して複数のモデルを評価できる別の方法は何ですか? 可能であれば、グローバルを避けたいと思います。

4

1 に答える 1

2

質問に答える前に、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 を使用してサンプル平均を 2D 配列に集約し、Numpy のパーセンタイル関数の軸オプションを使用して、各列に沿ってパーセンタイルを計算しました。

于 2013-11-09T16:56:16.677 に答える