Метрополис-Гастингс для гауссовой смеси и определения параметров - PullRequest
Новая
       24

Метрополис-Гастингс для гауссовой смеси и определения параметров

0 голосов
/ 16 марта 2019

Итак, у меня есть смесь двух гауссов, и я хочу определить средства (сигмы и веса известны).Для этого я пытаюсь реализовать алгоритм Метрополиса-Хастинга как метод MCMC.Я пытаюсь воспроизвести метод с этого сайта https://python4mpia.github.io/fitting_data/Metropolis-Hastings.html, но я не уверен, что я делаю правильный путь, потому что он очень медленный, и я не смог увидеть хорошие результаты.

Я думаю, что проблема заключается в определении логарифмического правдоподобия на каждой итерации. Я видел на этих слайдах https://www.ceremade.dauphine.fr/~xian/BCS/Bmix.pdf (слайд 296), что он должен быть вычисляем в O (п).Также я пытаюсь выяснить, может ли этот связанный вопрос помочь мне Ускорить Метрополис - Гастингс в Python , но не нашел ответов.

Мой вопрос, есть лиспособ значительно ускорить код ниже, особенно есть ли вероятность, что логарифмическая вероятность может быть вычислена быстрее? 

import matplotlib.pyplot as plt
import numpy as np
import math 

# define posterior distribution
def posterior(m1, m2, x):
    #mu, sigma = 5, 2.0 # mean and standard deviation that could be chosen
    num1 = math.exp( - ( x - m1 )**2 / ( 2.0 * sigma1 **2 ) )
    den1 = math.sqrt( 2 * math.pi * sigma1 **2)
    num2 = math.exp( - ( x - m2 )**2 / ( 2.0 * sigma2 **2 ) )
    den2 = math.sqrt( 2 * math.pi * sigma2 **2)
    return  num1 / den1+ num2/den2


# Generate toy Data
N = 10000
s = 10
x_0 = 0
mu1= 5
sigma1=1.0
mu2=1
sigma2=1.0

# Metroplis Hasting for sampling 
def Metroplis_sampling(N, s, x_0):

        p = posterior(mu1, mu2,  x_0)

        samples = []
        x = x_0

        for i in range(N):
            xn = x + np.random.normal(size=1)
            pn = posterior(mu1,  mu2, xn)
            if pn >= p:
                p = pn
                x = xn
            else:
                u = np.random.rand()
                if u < pn/p:
                    p = pn
                    x = xn
            if i % s == 0:
                samples.append(x)

        samples = np.array(samples)

        return samples

Samples = Metroplis_sampling(N, s, x_0)


#Define probabilities for samples
def probabilities(mus, samples):
    mu_1, mu_2 = mus[0], mus[1]
    probas = np.asarray( [math.log(posterior(mu_1, mu_2,  x)) for x in samples] )
    # return log likelihood.
    return np.sum(probas)



# initial guess for mu's as array.
guess = [4.5, 1.5]
# Prepare storing MCMC chain as array of arrays.
A = [guess]
# define stepsize of MCMC.
stepsizes = [0.005, 0.005]  # array of stepsizes
accepted  = 0.0
new_loglik = probabilities(guess, samples)
# Metropolis-Hastings with 10,000 iterations.
for n in range(10000):
    old_mus  = A[len(A)-1]  # old parameter values as array
    old_loglik = new_loglik
    # Suggest new candidate from Gaussian proposal distribution.
    new_mus = numpy.zeros([len(old_mus)])
    for i in range(len(old_mus)):
        # Use stepsize provided for every dimension.
        new_mus[i] = random.gauss(old_mus[i], stepsizes[i])
    new_loglik = probabilities(new_mus, samples)
    # Accept new candidate in Monte-Carlo fashing.
    if (new_loglik > old_loglik):
        A.append(new_mus)
        accepted = accepted + 1.0  # monitor acceptance
    else:
        u = random.uniform(0.0,1.0)
        if (u < math.exp(new_loglik - old_loglik)):
            A.append(new_mus)
            accepted = accepted + 1.0  # monitor acceptance
        else:
            A.append(old_mus)

print ("Acceptance rate = "+str(accepted/10000.0))
...