Up

Random Variation of Posterior Distribution

Simulation with PyStan using pickle

Posterior distribution depends on data, which is considered to be sampled from a population. Since data varies depending on sampling, posterior distribution also varies corresponding to the data. To investigate variation of posterior distribution, Python scripts with Stan were developed for simulation in Ubuntu.

As a population, two-dimensional multivariate normal distribution with means 0, variances 1, and correlation coefficient 0.5 was used. Each set of data of size n=100 was sampled and analyzed by Bayesian method. Each of Bayesian analysis produced a posterior distribution. This experiment of sampling and analysis was done 20 times, so produced 20 posterior distributions. Posterior distributions of correlation coefficient were displayed in Figure 1.

Figure 1. Posterior distributions of correlation coefficient for sample size 100

The Bayesian analysis was done by Stan script in Listing 1, which is used by the Python script in Listing2.

Listing 1. Stan script for two-dimensional multivariate normal distribution (file name stanBiNormal.stan).

data{

int<lower=2> N;

vector[2] XY[N];

}

parameters{

real mux;

real muy;

real<lower=0> sigmax;

real<lower=0> sigmay;

real<lower = -1, upper = 1> rho;

}

transformed parameters {

vector[2] mu;

cov_matrix[2] Sgm;

mu[1] = mux;

mu[2] = muy;

Sgm[1][1] = sigmax * sigmax;

Sgm[1][2] = rho * sigmax * sigmay;

Sgm[2][1] = rho * sigmax * sigmay;

Sgm[2][2] = sigmay * sigmay;

}

model{

for (i in 1:N)

XY[i] ~ multi_normal(mu, Sgm);

}

The Python script, which used the Stan script in Listing1, is shown in Listing 2. Run the script in Listing 2, the graph in Figure 1 will be shown. In sampling by MCMC, the Stan script is compiled only once, and the constructed binary code is saved in a file by pickle, which is loaded by pickle again to use for sampling by MCMC. That is, compiling is one time, sampling is repeated using the compiled code saved in a file. The Python script was run in Ubuntu.

Listing 2. A Python script, which uses the Stan script in Listing 1.

import numpy as np

import numpy.random as nprn

import pystan

from pystan import StanModel

import pickle

import matplotlib.pyplot as plt

import seaborn as sb

import time

rho = 0.5

n = 100

sm = StanModel(file = 'stanBiNormal.stan')

with open('model.pkl', 'wb') as f:

pickle.dump(sm, f)

for i in range(20):

sm = pickle.load(open('model.pkl', 'rb'))

xy = nprn.multivariate_normal([0, 0], [[1, rho], [rho, 1]], n)

Data = {'N': n, 'XY': xy}

print('\n\nSampling started...i = {}\n'.format(i))

time.sleep(1)

def f_init():

return dict(mux = 0.0, muy = 0.0, rho = 0.0,

sigmax = 1.0, sigmay = 1)

fit = sm.sampling(data = Data, init = f_init)

DistRho = fit['rho']

sb.kdeplot(DistRho, clip = [-1, 1], alpha = 0.5, color = 'b')

plt.hlines([0], -1, 1, 'k')

plt.plot([rho], [0], 'ro', markersize = 10)

plt.title('Distribution of \nPosterior distributions of $\\rho$' +

', n = {}'.format(n))

plt.show()

The posterior distributions in Figure 1 are satisfactory in that they cover the parameter true value 0.5 shown by the small red disk. However, when we set the data sample size to be 3, we get the posterior distributions shown in Figure 2.

Figure 2. Posterior distributions for sample size 3

Compared with the distributions in Figure 1, the posterior distributions in Figure 2 are unsatisfactory. Sample size as small as 3 provides little information about correlation.

Up