2

特定の確率でバイナリ観測を生成するAR(1) プロセスをモデル化しようとしています。

パーティクル フィルター(順次重要度リサンプリング)を実装しましたが、その動作は奇妙に思えます。


いくつかのヘルパー関数を定義した後:

sigmoid2 = lambda x: 1/(1+np.exp(-x))
bern_likelihood = lambda x,y: ( (sigmoid2(x))**y ) * ( (1-sigmoid2(x))**(1-y) )
normal = lambda x, mu, sigma: (1/(sigma*np.sqrt(2*pi)))*np.exp((-(1/2)*((x-mu)**2))/(sigma**2))
newVal = lambda nseT:  np.random.normal(loc=0,scale=nseT,size=nParticles)

潜在変数を生成します。

thetaPath = []
x = 0​
AR_noise = .1
​ ​
nTrials = 800
for i in range(nTrials):
    x += np.random.normal(loc=0,scale=AR_noise)
    thetaPath.append(x)

thetaPath = np.array(thetaPath)
​

これは、シグモイド リンクを介して、観測されたバイナリ データを生成します。

frac_corr = sigmoid2(thetaPath)
trials = np.random.random(size=nTrials)<frac_corr

次に、ここに私の粒子フィルター コードを示します。

stL = 0
nParticles = 500

#AR_noise = 0.1
dummy = arange(nParticles)

x0 = np.random.normal(loc=0,scale=AR_noise,size=nParticles)   #initialise a set of particles
w0_unNorm = bern_likelihood(x0,trials[0])    #calculate their weights
w0 = w0_unNorm/np.sum(w0_unNorm)    #normalise their weights
xi = cp.deepcopy(x0)
wi = cp.deepcopy(w0)
particle_pths = np.zeros([nParticles,nTrials])
particle_pths[:,0] = cp.deepcopy(x0)
for i in range(1,nTrials):


    #This part does resampling
    if (1/(np.sum(wi**2)))<(nParticles/2):
        #sample with replacement
        resamp =np.random.choice(dummy,size=nParticles,replace=True,p=wi)

        stL += 1   #simple counter for number of times that resampled
        xi = cp.deepcopy(xi[resamp])
        wi = [1/nParticles]*nParticles


    oldxi = xi
    xi += newVal(AR_noise*2)
    wi *= bern_likelihood(xi,trials[i])*normal(xi,oldxi,AR_noise)/normal(xi,oldxi,AR_noise*2)
    wi /= np.sum(wi)

   particle_pths[:,i] = xi

ここで、私の提案は次のように定義されます

normal(xi,oldxi,AR_noise*2)

一方、実際のモデルは

normal(xi,oldxi,AR_noise)

粒子の数が多い場合も少ない場合も、基になる状態の推定値が離散的で不安定であることがわかります。両方の画像の加重平均の推定値: ここに画像の説明を入力

リサンプリングしていないときは、これは表示されません。 ここに画像の説明を入力


この動作の理由として何が考えられますか? それは予想通りですか、それともコードに誤りがありますか?

(提供されたコードは、関連するモジュールをインストールした後、箱から出してすぐに動作するはずです)

4

0 に答える 0