人口モンテカルロの実装

Question

私は、単純なモデル（関数model()を参照）のために、this紙（79ページ図3参照）に記述されているPopulation Monte CarloアルゴリズムをPythonを使って1つのパラメータで実装しようとしています。残念ながら、このアルゴリズムはうまくいかず、何が間違っているのか分かりません。以下の私の実装を参照してください。実際の関数はabc()です。他のすべての関数はヘルパー関数として見ることができ、正常に動作するようです。

アルゴリズムが働くかどうかをチェックするには、モデルの唯一のパラメータをparam = 8に設定して観測データを生成します。したがって、ABCアルゴリズムの結果の後部は、8付近に中心を置く必要があります。なぜ私は思っています。

ご意見やご感想をお寄せいただきありがとうございます。私は以来、PMCアルゴリズムの神託の実装を探していたとして、私はこの質問につまずい

# imports 

from math import exp 
from math import log 
from math import sqrt 
import numpy as np 
import random 
from scipy.stats import norm 


# globals 
N = 300    # sample size 
N_PARTICLE = 300  # number of particles 
ITERS = 5   # number of decreasing thresholds 
M = 10    # number of words to remember 
MEAN = 7    # prior mean of parameter 
SD = 2    # prior sd of parameter 


def model(param): 
    recall_prob_all = 1/(1 + np.exp(M - param)) 
    recall_prob_one_item = np.exp(np.log(recall_prob_all)/float(M)) 
    return sum([1 if random.random() < recall_prob_one_item else 0 for item in range(M)]) 

## example 
print "Output of model function: \n" + str(model(10)) + "\n" 

# generate data from model 
def generate(param): 
    out = np.empty(N) 
    for i in range(N): 
    out[i] = model(param) 
    return out 

## example 

print "Output of generate function: \n" + str(generate(10)) + "\n" 


# distance function (sum of squared error) 
def distance(obsData,simData): 
    out = 0.0 
    for i in range(len(obsData)): 
    out += (obsData[i] - simData[i]) * (obsData[i] - simData[i]) 
    return out 

## example 

print "Output of distance function: \n" + str(distance([1,2,3],[4,5,6])) + "\n" 


# sample new particles based on weights 
def sample(particles, weights): 
    return np.random.choice(particles, 1, p=weights) 

## example 

print "Output of sample function: \n" + str(sample([1,2,3],[0.1,0.1,0.8])) + "\n" 


# perturbance function 
def perturb(variance): 
    return np.random.normal(0,sqrt(variance),1)[0] 

## example 

print "Output of perturb function: \n" + str(perturb(1)) + "\n" 

# compute new weight 
def computeWeight(prevWeights,prevParticles,prevVariance,currentParticle): 
    denom = 0.0 
    proposal = norm(currentParticle, sqrt(prevVariance)) 
    prior = norm(MEAN,SD) 
    for i in range(len(prevParticles)): 
    denom += prevWeights[i] * proposal.pdf(prevParticles[i]) 
    return prior.pdf(currentParticle)/denom 


## example 

prevWeights = [0.2,0.3,0.5] 
prevParticles = [1,2,3] 
prevVariance = 1 
currentParticle = 2.5 
print "Output of computeWeight function: \n" + str(computeWeight(prevWeights,prevParticles,prevVariance,currentParticle)) + "\n" 


# normalize weights 
def normalize(weights): 
    return weights/np.sum(weights) 


## example 

print "Output of normalize function: \n" + str(normalize([3.,5.,9.])) + "\n" 


# sampling from prior distribution 
def rprior(): 
    return np.random.normal(MEAN,SD,1)[0] 

## example 

print "Output of rprior function: \n" + str(rprior()) + "\n" 


# ABC using Population Monte Carlo sampling 
def abc(obsData,eps): 
    draw = 0 
    Distance = 1e9 
    variance = np.empty(ITERS) 
    simData = np.empty(N) 
    particles = np.empty([ITERS,N_PARTICLE]) 
    weights = np.empty([ITERS,N_PARTICLE]) 

    for t in range(ITERS): 
    if t == 0: 
     for i in range(N_PARTICLE): 
     while(Distance > eps[t]): 
      draw = rprior() 
      simData = generate(draw) 
      Distance = distance(obsData,simData) 

     Distance = 1e9 
     particles[t][i] = draw 
     weights[t][i] = 1./N_PARTICLE 

     variance[t] = 2 * np.var(particles[t]) 
     continue 


    for i in range(N_PARTICLE): 
     while(Distance > eps[t]): 
     draw = sample(particles[t-1],weights[t-1]) 
     draw += perturb(variance[t-1]) 
     simData = generate(draw) 
     Distance = distance(obsData,simData) 

     Distance = 1e9 
     particles[t][i] = draw 
     weights[t][i] = computeWeight(weights[t-1],particles[t-1],variance[t-1],particles[t][i]) 


    weights[t] = normalize(weights[t]) 
    variance[t] = 2 * np.var(particles[t]) 

    return particles[ITERS-1] 



true_param = 9 
obsData = generate(true_param) 
eps = [15000,10000,8000,6000,3000] 
posterior = abc(obsData,eps) 
#print posterior 

Answer 1

は、非常に偶然、私は私自身の研究にこの正確な論文で技術を適用する過程で、現在です。

あなたが得た結果を投稿できますか？私の推測では、1）距離関数（および/または類似性の閾値）の貧弱な選択を使用している、または2）十分な粒子を使用していないということです。私はここで間違っているかもしれませんが（私はサンプルの統計にあまり精通していませんが）、あなたの距離関数は暗黙のうちにあなたのランダムの順序が重要であることを私に示唆しています。私は実際にコンバージェンスのプロパティに影響を与えるかどうかを判断するためにこれ以上検討する必要がありますが、単純に平均値または中央値をサンプルの統計値として使用しないでください。

私の距離関数としての標本平均と[0.5、0.3、0.1]のεを用いた3回の反復の絶対差を使って、1000個の粒子と真のパラメータ値を持つコードを実行しました。私の推定された事後分布のピークは、各反復ごとにそうであるように8に近づいているように見えますが、母集団の分散が減少します。それでもなお顕著な右バイアスが存在することに注意してください。これはモデルの非対称性のためです（パラメータ値8以下では8回以上の成功は得られませんが、8より大きいパラメータ値はすべて右につながります）分布の偏り）。ここで

が私の結果をプロットします：MCの仕組みに慣れていない人があなたのコードを通過するため