Home

Python Multiprocessing Programming for MCMC

 

A sample program was written in Python, using multiprocessing, so that multiple chains in MCMC were run concurrently. A model for one variable normal distribution was employed, that is, it was assumed that data s were sampled from a normal distribution

Hence, we have

 is called a likelihood function and denoted by , when it is considered to be a function of parameters  and with data s given and fixed. That is,

Hence, the posterior distribution is given by the following equation

Uniform distribution was used as a prior distribution

where,

 

For equation (1), the following algorithm, which is known as Metropolis-within-Gibbs, was employed.

 

Proposal distributions were chosen as follows

and were adjusted in an adapting stage before the main MCMC stage as follows (see the code in mainprg.py below)

 

        v, ck = scale_sgm(r_acpt_mean, ck)

        params['genMuSgm'] = params['genMuSgm'] * v

        v, ck = scale_sgm(r_acpt_sd, ck)

        params['genSDSgm'] = params['genSDSgm'] * v

 

params['genMuSgm'] and params['genSDSgm'] represent  and , r_acpt_mean and r_acpt_sd are acceptance rates for parameters and The function scale_sgm is given

def scale_sgm( acpt_r, ck ):

    v = 1.0

    if (acpt_r > 0.5):

        ck += 1

        v = 1.0 + 2.0 * (acpt_r - 0.4) / 0.6

    elif (acpt_r < 0.3):

        ck += 1

        v = 1.0 / (1.0 + 2.0 * (0.4 - acpt_r) / 0.4)

 

    return v, ck

 

The function scale_sgm adjusts  and , so that acceptance rates becomes from 0.3 to 0.5

 

Cycles of MCMC were executed according to the following steps

 

Step 0

Set  ,and for the first MCMC, the initial values  and  are set to be the mean and SD of the data. The initial values  and  are also given based on the SD of the data. For MCMC after the first MCMC, the initial values are set based on samples from the previous MCMC.

 

Step 1

As a candidate for , sample y as

Set the acceptance probability  as follows

Set 

      with probability  

             or

    with probability  

 

Step 2

As a candidate for , sample y as

Set acceptance probability  as follows

Set

      with probability

or

    with probability

 

Step 3

Set

If t reaches the maximum number of cycles, then exit. Otherwise, jump back to step 1.

 

 

Some tips on calculation are given below.

 

The following term in the equation in step1

is given by

because

 

When we calculate a ratio of probabilities

Values of probability may become 0, because of the range, which type double can represent. To consider this possibility, the following code is prepared.

 

v = 1.0

Lv = 0.0

for i in range(0, 325):

    v *= 0.1

    print("%d    v = %0.17e"%((i + 1), v))

 

Run this code, then we have the results in Figure 1.

Figure 1

 

Beyond the value , v becomes 0. To avoid to get value 0 for probability, calculations are based on log-transformed value. That is, the following transformation was employed.

In this case, we have

 

To restore the log-transformed value to the not-transformed one, the function exp may be used. But, as to the function exp, some caution should be needed. Consider the following code.

 

import math

 

v = 700.0

for i in range(0, 10):

    try:

        print("exp(%f) = %0.15e    exp(%f) = %0.15e"%(\

            v, math.exp(v), -v, math.exp(-v)))

        v += 3.0

    except Exception as e:

        print(e)

        print("v = %f (exp overflow)     exp(%f) = %0.15f"%(\

              v, -v, math.exp(-v)))

 

Run this code, then we get the results in Figure 2.

Figure 2

 

For a large value v, exp(v) causes overflow error.

Acceptance probability  need not become larger than 1 (see the code of the function mcmc_func in subprg.py below), so the following function exp0 was prepared:

 

def exp0( x ):

    if (x > 0.0):

        return 1.0

    else:

        return math.exp(x)

 

In step 2, log-transformed terms are given as follows.

For a candidate y such as

we have

In this case, the acceptance probability  is set to be -1, which means that actual acceptance probability is 0 (see the code of the function mcmc_func in subprg.py below).

 

The algorithm of the MCMC is implemented as the function mcmc_func in the following code file subprg.py.

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝  subprg.py  ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

import math

from RNGen import RNp2to191 as RN

 

def sqr( v ):

    return v * v

 

def exp0( x ):

    if (x > 0.0):

        return 1.0

    else:

        return math.exp(x)

 

def scale_sgm( acpt_r, ck ):

    v = 1.0

    if (acpt_r > 0.5):

        ck += 1

        v = 1.0 + 2.0 * (acpt_r - 0.4) / 0.6

    elif (acpt_r < 0.3):

        ck += 1

        v = 1.0 / (1.0 + 2.0 * (0.4 - acpt_r) / 0.4)

 

    return v, ck

   

 

def mcmc_func( params ):

    """

          params = {'jump':  'n':  'x':  'NSimu':

                    'init_mean':  'init_sd':  'genMuSgm': 'geSDSgm':  }

    """

    jump      = params['jump']

    n         = params['n']

    x         = params['x']

    NSimu     = params['NSimu']

    init_mean = params['init_mean']

    init_sd   = params['init_sd']

    genMuSgm  = params['genMuSgm']

    genSDSgm  = params['genSDSgm']

    sumX = 0.0

    ssumX = 0.0

    for i in range(0, n):

        sumX += x[i]

        ssumX += x[i] ** 2.0

    mean = [init_mean]

    sd = [init_sd]

    acpt_mean = 0

    acpt_sd = 0

 

    rn = RN()

    #

    #   Independent sequences can be obtained for different values of jump

    #

    for i in range(0, jump):

        rn.jump(100)

    #

    #       NSimu samples from MCMC

    #

    for t in range(0, NSimu):

        if (t % 1000) == 0:

            print("t/NSimu = %d/%d   for jump = %d"%(t, NSimu, jump))

 

        y = rn.normalMS(mean[t], genMuSgm)

        log_a = (-2.0 * sumX * (y - mean[t]) + n * (sqr(y) - sqr(mean[t]))) / (-2.0 * sqr(sd[t]))

        a = exp0(log_a)

        if (rn.uni() < a):

            mean.append(y)

            acpt_mean += 1

        else:

            mean.append(mean[t])

 

        y = rn.normalMS( sd[t], genSDSgm )

        a = 0.0

        if (y <= 0.0):

            a = -1.0

        else:

            log_a = n * math.log(sd[t] / y) \

                        + 0.5 * ((1.0 / sqr(sd[t])) - (1.0 / sqr(y))) \

                        * (ssumX - 2.0 * mean[t + 1] * sumX + n * sqr(mean[t + 1]))

            a = exp0(log_a)

        if (rn.uni() < a):

            sd.append(y)

            acpt_sd += 1

        else:

            sd.append(sd[t])

 

    return mean, acpt_mean, sd, acpt_sd

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

The following code makes Markov chains of child processes independent of each other.

 

    for i in range(0, jump):

        rn.jump(100)

 

As to the method jump, visit <this website>.

 

Sampling distributions from the posterior distributions are displayed as histograms (Figure 3 for ; Figure 4 for ).

Figure 3

 

Figure 4

 

Quartile values and 95% CIs are written in the output file (Figure 5).

Figure 5

 

The main program code is shown below.

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝  mainprg.py  ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

from multiprocessing.pool import Pool

import matplotlib.pyplot as plt

import math

from subprg import scale_sgm, mcmc_func

 

if __name__ == "__main__":

    #

    #       Prepare the input data file

    #

    s_in = input("Input data file = ")

    f_in = open(s_in, "r")

    #

    #       Set the contents of the input data file in the object data

    #

    data_f = f_in.readlines()

    f_in.close()

    #

    #       Prepare the output file

    #

    s_out = input("Output data file = ")

    f_out = open(s_out, "w")

    f_out.write("Input data file = " + s_in + "\n\n")

    #

    #           Set the data in list x

    #

    pos = 0

    while True:

        if data_f[pos].find("/", 0) == 0: break

        pos += 1

    n = 0

    x = []

    while True:

        pos += 1

        if data_f[pos].find("/", 0) == 0: break

        x.append(float(data_f[pos]))

        n += 1

       

    print("X...")

    i = 1

    for v in x:

        f_out.write("%5d:"%i + "  %f\n"%v)

        i+=1

 

    f_out.write("\n")

    sum = 0.0

    for v in x:

        sum += v

    mean = sum / n

    print("Mean = %f"%mean)

    f_out.write("Mean = %f\n"%mean)

    ssum = 0.0

    for v in x:

        ssum += (v - mean) ** 2.0

    var = ssum / n

    sd = math.sqrt(var)

    print("Variance = %f"%var + "     sd = %f"%sd)

    f_out.write("Variance = %f"%var + "     sd = %f\n"%sd)

    uvar = ssum / (n - 1.0)

    sd_uvar = math.sqrt(uvar)

    print("Unbiased Var. = %f"%uvar + "     sd(uvar) = %f"%sd_uvar)

    f_out.write("Unbiased Var. = %f"%uvar + "     sd(uvar) = %f\n"%sd_uvar)

 

    NSimu = 3000

    #

    #       Initial values of the parameters of MCMC

    #

    params = {'jump': 0,

              'n': n,

              'x': x,

              'NSimu': NSimu,

              'init_mean': mean,

              'init_sd': sd,

              'genMuSgm': 0.1 * sd,

              'genSDSgm': 0.1 * sd }

    #

    #           Adjust the parameters of MCMC

    #

    ck = -1

    n_adap = 0

    while True:

        print("ck = %d   n_adap = %d"%(ck, n_adap))

        L_mean, acpt_mean, L_sd, acpt_sd = mcmc_func( params )

        sum_mean = 0.0

        sum_sd = 0.0

        for t in range(1, NSimu + 1):

            sum_mean += L_mean[t]

            sum_sd += L_sd[t]

        params['init_mean'] = sum_mean / NSimu

        params['init_sd'] = sum_sd / NSimu

 

        r_acpt_mean = acpt_mean / float(NSimu)

        r_acpt_sd = acpt_sd / float(NSimu)

       

 

        ck = 0

        v, ck = scale_sgm(r_acpt_mean, ck)

        params['genMuSgm'] = params['genMuSgm'] * v

        v, ck = scale_sgm(r_acpt_sd, ck)

        params['genSDSgm'] = params['genSDSgm'] * v

 

        if (ck == 0):

            print("ck = %d"%ck)

            break

 

        n_adap +=1

  

    #

    #           The main stage of MCMC

    #

    pool = Pool()

 

    NSimu = 10000

    params['NSimu'] = NSimu

    n_mcmc = 4              #   The number of concurrent MCMCs

    ary_params = []

    for i in range(0, n_mcmc):

        temp = params.copy()

        temp['jump'] = i + 1

        ary_params.append(temp)

 

    #

    #       Run the processes of MCMCs

    #

    print("The main MCMCs started.")

    results = pool.map(mcmc_func, ary_params)

    print("The main MCMCs ended.")

    #

    #       Merge the samples from MCMCs

    #

    mrg_mean = []

    mrg_sd = []

    for t_mean, t_acpt_mean, t_sd, t_acpt_sd in results:

        print("t_acpt_mean = %d    t_acpt_sd = %d"%(t_acpt_mean, t_acpt_sd))

        mrg_mean += t_mean[1:NSimu+1]

        mrg_sd += t_sd[1:NSimu+1]

 

    #

    #           Calculate the quantiles

    #

    mrg_mean.sort()

    L_mean = mrg_mean[int(NSimu * n_mcmc * 0.025)]

    Q1_mean = mrg_mean[int(NSimu * n_mcmc * 0.25)]

    Med_mean = mrg_mean[int(NSimu * n_mcmc * 0.5)]

    Q3_mean = mrg_mean[int(NSimu * n_mcmc * 0.75)]

    U_mean = mrg_mean[int(NSimu * n_mcmc * 0.975)]

    print("\nQ1_mean = %f     Med_mean = %f    Q3_mean = %f"%(Q1_mean, \

                                                            Med_mean, Q3_mean))

    print("95%%CI = [%f, %f]"%(L_mean, U_mean))

    f_out.write("\nQ1_mean = %f     Med_mean = %f    Q3_mean = %f\n"%(Q1_mean, \

                                                            Med_mean, Q3_mean))

    f_out.write("95%%CI = [%f, %f]\n"%(L_mean, U_mean))

 

    plt.figure(figsize = (12, 9))

    plt.title("Posterior Distribution of the Mean")

    plt.hist(mrg_mean, bins = 50)

    plt.show()

   

    mrg_sd.sort()

    L_sd = mrg_sd[int(NSimu * n_mcmc * 0.025)]

    Q1_sd = mrg_sd[int(NSimu * n_mcmc * 0.25)]

    Med_sd = mrg_sd[int(NSimu * n_mcmc * 0.5)]

    Q3_sd = mrg_sd[int(NSimu * n_mcmc * 0.75)]

    U_sd = mrg_sd[int(NSimu * n_mcmc * 0.975)]

    print("\nQ1_sd = %f     Med_sdn = %f    Q3_sd = %f"%(Q1_sd, \

                                                            Med_sd, Q3_sd))

    print("95%%CI = [%f, %f]"%(L_sd, U_sd))

    f_out.write("\nQ1_sd = %f     Med_sd = %f    Q3_sd = %f\n"%(Q1_sd, \

                                                            Med_sd, Q3_sd))

    f_out.write("95%%CI = [%f, %f]\n"%(L_sd, U_sd))

   

    plt.figure(figsize = (12, 9))

    plt.title("Posterior Distribution of the SD")

    plt.hist(mrg_sd, bins = 50)

    plt.show()

 

    f_out.close()

    print(s_out + " was saved.")

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

Proposal distributions  and  of the MCMC are adjusted automatically by the following code.

    ck = -1

    n_adap = 0

    while True:

        print("ck = %d   n_adap = %d"%(ck, n_adap))

        L_mean, acpt_mean, L_sd, acpt_sd = mcmc_func( params )

        sum_mean = 0.0

        sum_sd = 0.0

        for t in range(1, NSimu + 1):

            sum_mean += L_mean[t]

            sum_sd += L_sd[t]

        params['init_mean'] = sum_mean / NSimu

        params['init_sd'] = sum_sd / NSimu

 

        r_acpt_mean = acpt_mean / float(NSimu)

        r_acpt_sd = acpt_sd / float(NSimu)

       

 

        ck = 0

        v, ck = scale_sgm(r_acpt_mean, ck)

        params['genMuSgm'] = params['genMuSgm'] * v

        v, ck = scale_sgm(r_acpt_sd, ck)

        params['genSDSgm'] = params['genSDSgm'] * v

 

        if (ck == 0):

            print("ck = %d"%ck)

            break

 

        n_adap +=1

 

 

Multiprocessing is implemented for the main MCMC using the class Pool as follows. Multiprocessing by the class Pool is explained at the website <Multiprocessing by the Class Pool of Python>.

 

#

    #           The main stage of MCMC

    #

    pool = Pool()

 

    NSimu = 10000

    params['NSimu'] = NSimu

    n_mcmc = 4              #   The number of concurrent MCMCs

    ary_params = []

    for i in range(0, n_mcmc):

        temp = params.copy()

        temp['jump'] = i + 1

        ary_params.append(temp)

 

    #

    #       Run the processes of MCMCs

    #

    print("The main MCMCs started.")

    results = pool.map(mcmc_func, ary_params)

    print("The main MCMCs ended.")

 

 

Run the program, the name of the input data file is asked (Figure 6).

Figure 6

 

Data in the input data file are enclosed by the lines with slash / at the heads as shown in Figure 7.

Figure 7

 

Set the name of the input data file, and press the Enter key. A name of the output file will be asked (Figure 8).

Figure 8

 

The name of the output file is an arbitrary text file name. After setting the output file name, calculation starts. When the main MCMC ends, the histogram of the samples from the posterior distribution of  will be presented (Figure 3). When the window of the histogram of  is closed by clicking the X icon at the upper right corner, the window of the histogram for  will be shown (Figure 4).

Close the window of the histogram for , then the program ends.

Output by the function print in child processes are not displayed in the window of the IDLE shell. When the program is executed in a console window (a terminal window), outputs by the child processes will be displayed (Figure 9).

Figure 9

 

Program files are archived in the zip file mcmcMultiProcPrg.zip , which can be downloaded by clicking the name mcmcMultiProcPrg.zip .

 

 

Home