smc - py - prot

gistfile1.py

def smc(M, steps, prior, potential, proposal_kernel, correction_steps=1, ess_ratio=2):
    """
    smc(M, steps, prior, potential, proposal_kernel, correction_steps=1, ess_ratio=2)
    
    Creates a particle approximation of a probability distribution using the Sequential
    Monte Carlo (SMC) algorithm with Markov Chain Monte Carlo (MCMC) correction steps using
    the Metropolis-Hastings algorithm.
    
    Parameters
    ----------
    M : int
        The number of particles to use for the approximation.
    steps : int
        Number of intermediary distributions used for approximating the target
        distribution. Each intermediary distribution will be used to perform an
        update of the current particle approximation.
    prior : scipy.stats.rv_continuous
        Prior distribution, used to sample an initial set of particles
        and later used for the acceptance of the MCMC proposals.
    potential : (np.ndarray, int) -> np.ndarray
        Function representing the sequence of potentials used to
        approximate the target distribution. The function should accept
        as first argument the array of particles for which to evaluate
        the potential and as a second argument the index of the
        potential to evaluate. It should return an array of evaluations of
        the potential.
    proposal_kernel : (np.ndarray, int) -> np.ndarray
        Function sampling proposals for the Metropolis-Hastings agorithm.
        The function should accept the current array of particles as its first
        argument and the current step and return an array of proposal particles. 
    correction_steps : int, optional
        The number of correction steps to apply to the set of particles for each
        update of the approximation.
    ess_ratio : int, optional
        Specifies the lower bound on the effective sample size. A resampling step
        is then performed if `ESS < M/ess_ratio`.
    """
    
    # initial set of particles is created using a Monte Carlo approximation of the prior distribution
    x = prior.rvs(size=M)
    weights = np.full(M,1/M)
    
     
    for n in range(steps+1):
        weights = np.log(weights)
        
        proposals = proposal_kernel(x)
        proposal_weights = potential(proposals, n)
        
        current_pot = potential(x, n)
        
        np.add(weights, current_pot - potential(x,n-1), out=weights)
        
        for j in range(correction_steps):
            acceptance = stats.uniform.rvs(size=M)
            np.log(acceptance, out=acceptance)

            potential_ratio = proposal_weights - current_pot
            prior_ratio = prior.logpdf(proposals) - prior.logpdf(x)
            
            accepted = acceptance < (potential_ratio + prior_ratio)
            x[accepted] = proposals[accepted]
            
            if j < correction_steps-1:
                not_accepted = np.logical_not(accepted, out=accepted)
                proposal_weights[not_accepted] = current_pot[not_accepted]

                current_pot = proposal_weights
                
                proposals = proposal_kernel(x)
                proposal_weights = potential(proposals, n)
    
        np.exp(weights, out=weights)
        np.divide(weights, weights.sum(), out=weights)
        
        ess = 1 / np.dot(weights, weights)
        print(ess)
        if ess < M/ess_ratio:
            x = np.random.choice(x, size=M, p=weights, replace=True)
            weights = np.full(M,1/M)
        
    return x,weights

#x,w = smc(1000, 11, prior, vpotential, gaussian_proposal, correction_steps=5)
plt.hist(x,weights=w)
np.dot(x,w)