Markov Chain montecarlo integration and infinite while loop

Question

I'm implementing a Markov Chain Montecarlo with metropolis and barkes alphas for numerical integration. I've created a class called MCMCIntegrator(). I've loaded it with some attributes, one of then is the pdf of the function (a lambda) we're trying to integrate called g.

import numpy as np
import scipy.stats as st


class MCMCIntegrator:

    def __init__(self):

        self.g = lambda x: st.gamma.pdf(x, 0, 1, scale=1 / 1.23452676)*np.abs(np.cos(1.123454156))
        self.size = 10000
        self.std = 0.6
        self.real_int = 0.06496359

There are other methods in this class, the size is the size of the sample that the class must generate, std is the standard deviation of the Normal Kernel, which you will see in a few seconds. The real_int is the value of the integral from 1 to 2 of the function we're integrating. I've generated it with a R script. Now, to the problem.

 def _chain(self, method=None):

        """
            Markov chain heat-up with burn-in

        :param method: Metrpolis or barker alpha
        :return: np.array containing the sample
        """
        old = 0
        sample = np.zeros(int(self.size * 1.5))
        i = 0

        if method:
            def alpha(a, b): return min(1, self.g(b) / self.g(a))

        else:
            def alpha(a, b): return self.g(b) / (self.g(a) + self.g(b))

        while i != len(sample):
            new = st.norm(loc=old, scale=self.std).rvs()
            new = abs(new)
            al = alpha(old, new)
            u = st.uniform.rvs()

            if al > u:
                sample[i] = new
                old = new
                i += 1

        return np.array(sample)

Below this method is an integrate() method that calculates the proportion of numbers in the [1, 2] interval:

    def integrate(self, method=None):
        """
            Integration step

        """

        sample = self._chain(method=method)
        
        # discarding 30% of the sample for the burn-in
        ind = int(len(sample)*0.3)
        sample = sample[ind:]
        setattr(self, "sample", sample)

        sample = [1 if 1 < v < 2 else 0 for v in sample]
        return np.mean(sample)

This is the main function:

def main():

    print("-- RESULTS --".center(20), end='\n')
    mcmc = MCMCIntegrator()
    print(f"\t{mcmc.integrate()}", end='\n')
    print(f"\t{np.abs(mcmc.integrate() - mcmc.real_int) / mcmc.real_int}")


if __name__ == "__main__":
    main()

I'm stuck in an infinity while loop and I have no idea why this is happening.

Answer 1

Couple things... You are hung up in the chain method because the alpha computation is returning NaN, because g() is returning NaN. Take a look at the print statements I inserted into your code and run it...

tips:

1. Make g() a class function just like chain.

2. Test g() on some test values, something is clearly amiss

3. Don't conditionally define a function like alpha. Wayyy confusing and error prone/tough to troubleshoot. Just pass alpha what it needs and you also can make it a class function alpha(a, b, method=None)
4. Take a look at where you are updating i in the `_chain' function.... You realize you are risking a long looping process because you only update i conditionally!
5. You are set up for disaster with your use of numpy array. You may have a bunch of trailing zeros after your actual data because you are over-writing the large zeros list. You do NOT need numpy array here, just use a python empty list and append new values to it, either zero or one...based on whatever.

Add a couple print statements in when you are troubleshooting (or unit test your functions). Try my adds to your function below... it is what I used to figure out what was going on

def _chain(self, method=None, verbose=True):

    """
        Markov chain heat-up with burn-in

    :param method: Metrpolis or barker alpha
    :return: np.array containing the sample
    """
    old = 0
    sample = np.zeros(int(self.size * 1.5))
    i = 0

    if method:
        def alpha(a, b): return min(1, self.g(b) / self.g(a))

    else:
        def alpha(a, b): 
            if verbose: print(f'g(a): {self.g(a)}, g(b): {self.g(b)}')
            return self.g(b) / (self.g(a) + self.g(b))

    while i != len(sample):
        new = st.norm(loc=old, scale=self.std).rvs()
        new = abs(new)
        al = alpha(old, new)
        u = st.uniform.rvs()
        if verbose: print(f'old: {old:.3f} new: {new:.3f} alpha: {al:.3f} u: {u:.3f}')
        if al > u:
            sample[i] = new
            old = new
            i += 1              # do you really want to conditionally update i?
        sys.exit(-1)            # to breakout of infinite loop...


    return np.array(sample)