特定の確率でバイナリ観測を生成する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)
粒子の数が多い場合も少ない場合も、基になる状態の推定値が離散的で不安定であることがわかります。両方の画像の加重平均の推定値:
この動作の理由として何が考えられますか? それは予想通りですか、それともコードに誤りがありますか?
(提供されたコードは、関連するモジュールをインストールした後、箱から出してすぐに動作するはずです)