Example Multilevel MCMC
Defines a hierarchy of simple Gaussian models and applies Multilevel MCMC to it.
# tell python where the MUQ libraries are installed
import sys
sys.path.insert(0, "/my/muq/dir/build/install/lib/")
# import ploting tools
import matplotlib as mpl
#mpl.use('TKAgg')
import matplotlib.pyplot as plt
# import numpy, which we use for linear algebra in python
import numpy as np
import h5py
import random
# import the MUQ libraries
import muq.Utilities as mu # import MUQ utilities module
import muq.Modeling as mm # import MUQ modeling module
import muq.Approximation as ma # import MUQ approximation module
import muq.SamplingAlgorithms as msa
class IdentityModpiece(mm.PyModPiece):
def __init__(self):
# initialize the muqModeling::ModPiece with zero inputs and one output, there is one output that has n components
#mm.PyModPiece.__init__(self, [], [n])
super(IdentityModpiece, self).__init__([2], [2])
def EvaluateImpl(self, inputs):
# this function does do anything (we have already computed the discretization)
# but ModPieces need an EvaluateImpl function
self.outputs = inputs
return
identityModpiece = IdentityModpiece()
cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
mu = np.array([.8, 2.3]) # mean
targetDensity = mm.Gaussian(mu, cov * 2.0).AsDensity()
cacheModPiece = mm.OneStepCachePiece(targetDensity)
problem_coarse = msa.SamplingProblem(cacheModPiece,cacheModPiece)
mu = np.array([1.0, 2.0]) # mean
targetDensity = mm.Gaussian(mu, cov).AsDensity()
cacheModPiece = mm.OneStepCachePiece(targetDensity)
problem_fine = msa.SamplingProblem(cacheModPiece,cacheModPiece)
nmcmc = 1000
trueCoeff = mu
# MCMC
options = dict()
options['NumSamples'] = nmcmc
options['PrintLevel'] = 3
options['KernelList'] = 'Kernel1'
options['Kernel1.Method'] = 'MHKernel'
options['Kernel1.Proposal'] = 'MyProposal'
options['Kernel1.MyProposal.Method'] = 'AMProposal'
options['Kernel1.MyProposal.InitialVariance'] = 0.1
options['Kernel1.MyProposal.AdaptSteps'] = 100
options['Kernel1.MyProposal.AdaptStart'] = 500
# create the MCMC sampler
mcmc = msa.SingleChainMCMC(options, problem_fine)
samps = mcmc.Run([trueCoeff])
print(mcmc.GetSamples().Mean())
print(mcmc.GetQOIs().Mean())
print(samps.Mean())
# MLMCMC
mioptions = dict()
mioptions['NumSamples_0'] = nmcmc * 10
mioptions['NumSamples_1'] = nmcmc
mioptions['Subsampling'] = 10
mioptions['Proposal.Method'] = 'AMProposal'
mioptions['Proposal.InitialVariance'] = 0.1
mioptions['Proposal.AdaptSteps'] = 100
mioptions['Proposal.AdaptStart'] = 500
mimcmc = msa.MIMCMC(mioptions, trueCoeff, [problem_coarse, problem_fine])
mimcmc.Run([trueCoeff])
print(mimcmc.MeanParam())
print(mimcmc.MeanQOI())
print(cacheModPiece.HitRatio())
# Let's dive into the MIMCMC data structures to do some post processing
# Iterate over all multiindices we have
mlindices = mimcmc.GetIndices()
for i in range(0,mlindices.Size()):
# Get the MIMCMCBox representing the telescoping sum component associated with this index
box = mimcmc.GetMIMCMCBox(mlindices.at(i))
boxIndices = box.GetBoxIndices()
print(mlindices.at(i).GetVector())
# Print the contribution to the multilevel telescoping sum of this box
print(box.MeanParam())
fig=plt.figure()
ax=fig.add_axes([0,0,1,1])
plt.xlim(-5, 8)
plt.ylim(-5, 8)
# Plot all samples of every MCMC chain of this box
for i in range(0,boxIndices.Size()):
boxIndex = boxIndices.at(i)
boxChain = box.GetChain(boxIndex)
samplesMat = boxChain.GetSamples().AsMatrix()
if boxIndex.GetValue(0) == 0:
ax.scatter(samplesMat[0,:], samplesMat[1,:], color='b', alpha=0.3, label="Coarse Samples")
else:
ax.scatter(samplesMat[0,:], samplesMat[1,:], color='r', alpha=0.3, label="Fine Samples")
# Let's access an individual sample
samples = boxChain.GetSamples()
print("Sample 10")
print(samples[10].state)
# and if it had a coarse proposal, get that as well:
if samples[10].HasMeta("coarseSample"):
print ("from coarse sample:")
print(samples[10].GetMetaSamplingState("coarseSample").state)
ax.legend()
# Plot all mean values of every MCMC chain of this box
for i in range(0,boxIndices.Size()):
boxIndex = boxIndices.at(i)
boxChain = box.GetChain(boxIndex)
mean = boxChain.GetSamples().Mean()
plt.axvline(x=mean[0], color='b' if boxIndex.GetValue(0) == 0 else 'r');
plt.axhline(y=mean[1], color='b' if boxIndex.GetValue(0) == 0 else 'r');
Completed code:
# tell python where the MUQ libraries are installed
import sys
sys.path.insert(0, "/my/muq/dir/build/install/lib/")
# import ploting tools
import matplotlib as mpl
#mpl.use('TKAgg')
import matplotlib.pyplot as plt
# import numpy, which we use for linear algebra in python
import numpy as np
import h5py
import random
# import the MUQ libraries
import muq.Utilities as mu # import MUQ utilities module
import muq.Modeling as mm # import MUQ modeling module
import muq.Approximation as ma # import MUQ approximation module
import muq.SamplingAlgorithms as msa
class IdentityModpiece(mm.PyModPiece):
def __init__(self):
# initialize the muqModeling::ModPiece with zero inputs and one output, there is one output that has n components
#mm.PyModPiece.__init__(self, [], [n])
super(IdentityModpiece, self).__init__([2], [2])
def EvaluateImpl(self, inputs):
# this function does do anything (we have already computed the discretization)
# but ModPieces need an EvaluateImpl function
self.outputs = inputs
return
identityModpiece = IdentityModpiece()
cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
mu = np.array([.8, 2.3]) # mean
targetDensity = mm.Gaussian(mu, cov * 2.0).AsDensity()
cacheModPiece = mm.OneStepCachePiece(targetDensity)
problem_coarse = msa.SamplingProblem(cacheModPiece,cacheModPiece)
mu = np.array([1.0, 2.0]) # mean
targetDensity = mm.Gaussian(mu, cov).AsDensity()
cacheModPiece = mm.OneStepCachePiece(targetDensity)
problem_fine = msa.SamplingProblem(cacheModPiece,cacheModPiece)
nmcmc = 1000
trueCoeff = mu
# MCMC
options = dict()
options['NumSamples'] = nmcmc
options['PrintLevel'] = 3
options['KernelList'] = 'Kernel1'
options['Kernel1.Method'] = 'MHKernel'
options['Kernel1.Proposal'] = 'MyProposal'
options['Kernel1.MyProposal.Method'] = 'AMProposal'
options['Kernel1.MyProposal.InitialVariance'] = 0.1
options['Kernel1.MyProposal.AdaptSteps'] = 100
options['Kernel1.MyProposal.AdaptStart'] = 500
# create the MCMC sampler
mcmc = msa.SingleChainMCMC(options, problem_fine)
samps = mcmc.Run([trueCoeff])
print(mcmc.GetSamples().Mean())
print(mcmc.GetQOIs().Mean())
print(samps.Mean())
# MLMCMC
mioptions = dict()
mioptions['NumSamples_0'] = nmcmc * 10
mioptions['NumSamples_1'] = nmcmc
mioptions['Subsampling'] = 10
mioptions['Proposal.Method'] = 'AMProposal'
mioptions['Proposal.InitialVariance'] = 0.1
mioptions['Proposal.AdaptSteps'] = 100
mioptions['Proposal.AdaptStart'] = 500
mimcmc = msa.MIMCMC(mioptions, trueCoeff, [problem_coarse, problem_fine])
mimcmc.Run([trueCoeff])
print(mimcmc.MeanParam())
print(mimcmc.MeanQOI())
print(cacheModPiece.HitRatio())
# Let's dive into the MIMCMC data structures to do some post processing
# Iterate over all multiindices we have
mlindices = mimcmc.GetIndices()
for i in range(0,mlindices.Size()):
# Get the MIMCMCBox representing the telescoping sum component associated with this index
box = mimcmc.GetMIMCMCBox(mlindices.at(i))
boxIndices = box.GetBoxIndices()
print(mlindices.at(i).GetVector())
# Print the contribution to the multilevel telescoping sum of this box
print(box.MeanParam())
fig=plt.figure()
ax=fig.add_axes([0,0,1,1])
plt.xlim(-5, 8)
plt.ylim(-5, 8)
# Plot all samples of every MCMC chain of this box
for i in range(0,boxIndices.Size()):
boxIndex = boxIndices.at(i)
boxChain = box.GetChain(boxIndex)
samplesMat = boxChain.GetSamples().AsMatrix()
if boxIndex.GetValue(0) == 0:
ax.scatter(samplesMat[0,:], samplesMat[1,:], color='b', alpha=0.3, label="Coarse Samples")
else:
ax.scatter(samplesMat[0,:], samplesMat[1,:], color='r', alpha=0.3, label="Fine Samples")
# Let's access an individual sample
samples = boxChain.GetSamples()
print("Sample 10")
print(samples[10].state)
# and if it had a coarse proposal, get that as well:
if samples[10].HasMeta("coarseSample"):
print ("from coarse sample:")
print(samples[10].GetMetaSamplingState("coarseSample").state)
ax.legend()
# Plot all mean values of every MCMC chain of this box
for i in range(0,boxIndices.Size()):
boxIndex = boxIndices.at(i)
boxChain = box.GetChain(boxIndex)
mean = boxChain.GetSamples().Mean()
plt.axvline(x=mean[0], color='b' if boxIndex.GetValue(0) == 0 else 'r');
plt.axhline(y=mean[1], color='b' if boxIndex.GetValue(0) == 0 else 'r');
Defines a hierarchy of simple Gaussian models and applies Multilevel MCMC to it.
# tell python where the MUQ libraries are installed
import sys
sys.path.insert(0, "/my/muq/dir/build/install/lib/")
# import ploting tools
import matplotlib as mpl
#mpl.use('TKAgg')
import matplotlib.pyplot as plt
# import numpy, which we use for linear algebra in python
import numpy as np
import h5py
import random
# import the MUQ libraries
import muq.Utilities as mu # import MUQ utilities module
import muq.Modeling as mm # import MUQ modeling module
import muq.Approximation as ma # import MUQ approximation module
import muq.SamplingAlgorithms as msa
class IdentityModpiece(mm.PyModPiece):
def __init__(self):
# initialize the muqModeling::ModPiece with zero inputs and one output, there is one output that has n components
#mm.PyModPiece.__init__(self, [], [n])
super(IdentityModpiece, self).__init__([2], [2])
def EvaluateImpl(self, inputs):
# this function does do anything (we have already computed the discretization)
# but ModPieces need an EvaluateImpl function
self.outputs = inputs
return
identityModpiece = IdentityModpiece()
cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
mu = np.array([.8, 2.3]) # mean
targetDensity = mm.Gaussian(mu, cov * 2.0).AsDensity()
cacheModPiece = mm.OneStepCachePiece(targetDensity)
problem_coarse = msa.SamplingProblem(cacheModPiece,cacheModPiece)
mu = np.array([1.0, 2.0]) # mean
targetDensity = mm.Gaussian(mu, cov).AsDensity()
cacheModPiece = mm.OneStepCachePiece(targetDensity)
problem_fine = msa.SamplingProblem(cacheModPiece,cacheModPiece)
nmcmc = 1000 trueCoeff = mu # MCMC options = dict() options['NumSamples'] = nmcmc options['PrintLevel'] = 3 options['KernelList'] = 'Kernel1' options['Kernel1.Method'] = 'MHKernel' options['Kernel1.Proposal'] = 'MyProposal' options['Kernel1.MyProposal.Method'] = 'AMProposal' options['Kernel1.MyProposal.InitialVariance'] = 0.1 options['Kernel1.MyProposal.AdaptSteps'] = 100 options['Kernel1.MyProposal.AdaptStart'] = 500 # create the MCMC sampler mcmc = msa.SingleChainMCMC(options, problem_fine) samps = mcmc.Run([trueCoeff]) print(mcmc.GetSamples().Mean()) print(mcmc.GetQOIs().Mean()) print(samps.Mean())
# MLMCMC mioptions = dict() mioptions['NumSamples_0'] = nmcmc * 10 mioptions['NumSamples_1'] = nmcmc mioptions['Subsampling'] = 10 mioptions['Proposal.Method'] = 'AMProposal' mioptions['Proposal.InitialVariance'] = 0.1 mioptions['Proposal.AdaptSteps'] = 100 mioptions['Proposal.AdaptStart'] = 500 mimcmc = msa.MIMCMC(mioptions, trueCoeff, [problem_coarse, problem_fine]) mimcmc.Run([trueCoeff]) print(mimcmc.MeanParam()) print(mimcmc.MeanQOI()) print(cacheModPiece.HitRatio())
# Let's dive into the MIMCMC data structures to do some post processing
# Iterate over all multiindices we have
mlindices = mimcmc.GetIndices()
for i in range(0,mlindices.Size()):
# Get the MIMCMCBox representing the telescoping sum component associated with this index
box = mimcmc.GetMIMCMCBox(mlindices.at(i))
boxIndices = box.GetBoxIndices()
print(mlindices.at(i).GetVector())
# Print the contribution to the multilevel telescoping sum of this box
print(box.MeanParam())
fig=plt.figure()
ax=fig.add_axes([0,0,1,1])
plt.xlim(-5, 8)
plt.ylim(-5, 8)
# Plot all samples of every MCMC chain of this box
for i in range(0,boxIndices.Size()):
boxIndex = boxIndices.at(i)
boxChain = box.GetChain(boxIndex)
samplesMat = boxChain.GetSamples().AsMatrix()
if boxIndex.GetValue(0) == 0:
ax.scatter(samplesMat[0,:], samplesMat[1,:], color='b', alpha=0.3, label="Coarse Samples")
else:
ax.scatter(samplesMat[0,:], samplesMat[1,:], color='r', alpha=0.3, label="Fine Samples")
# Let's access an individual sample
samples = boxChain.GetSamples()
print("Sample 10")
print(samples[10].state)
# and if it had a coarse proposal, get that as well:
if samples[10].HasMeta("coarseSample"):
print ("from coarse sample:")
print(samples[10].GetMetaSamplingState("coarseSample").state)
ax.legend()
# Plot all mean values of every MCMC chain of this box
for i in range(0,boxIndices.Size()):
boxIndex = boxIndices.at(i)
boxChain = box.GetChain(boxIndex)
mean = boxChain.GetSamples().Mean()
plt.axvline(x=mean[0], color='b' if boxIndex.GetValue(0) == 0 else 'r');
plt.axhline(y=mean[1], color='b' if boxIndex.GetValue(0) == 0 else 'r');
Completed code:
# tell python where the MUQ libraries are installed
import sys
sys.path.insert(0, "/my/muq/dir/build/install/lib/")
# import ploting tools
import matplotlib as mpl
#mpl.use('TKAgg')
import matplotlib.pyplot as plt
# import numpy, which we use for linear algebra in python
import numpy as np
import h5py
import random
# import the MUQ libraries
import muq.Utilities as mu # import MUQ utilities module
import muq.Modeling as mm # import MUQ modeling module
import muq.Approximation as ma # import MUQ approximation module
import muq.SamplingAlgorithms as msa
class IdentityModpiece(mm.PyModPiece):
def __init__(self):
# initialize the muqModeling::ModPiece with zero inputs and one output, there is one output that has n components
#mm.PyModPiece.__init__(self, [], [n])
super(IdentityModpiece, self).__init__([2], [2])
def EvaluateImpl(self, inputs):
# this function does do anything (we have already computed the discretization)
# but ModPieces need an EvaluateImpl function
self.outputs = inputs
return
identityModpiece = IdentityModpiece()
cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
mu = np.array([.8, 2.3]) # mean
targetDensity = mm.Gaussian(mu, cov * 2.0).AsDensity()
cacheModPiece = mm.OneStepCachePiece(targetDensity)
problem_coarse = msa.SamplingProblem(cacheModPiece,cacheModPiece)
mu = np.array([1.0, 2.0]) # mean
targetDensity = mm.Gaussian(mu, cov).AsDensity()
cacheModPiece = mm.OneStepCachePiece(targetDensity)
problem_fine = msa.SamplingProblem(cacheModPiece,cacheModPiece)
nmcmc = 1000
trueCoeff = mu
# MCMC
options = dict()
options['NumSamples'] = nmcmc
options['PrintLevel'] = 3
options['KernelList'] = 'Kernel1'
options['Kernel1.Method'] = 'MHKernel'
options['Kernel1.Proposal'] = 'MyProposal'
options['Kernel1.MyProposal.Method'] = 'AMProposal'
options['Kernel1.MyProposal.InitialVariance'] = 0.1
options['Kernel1.MyProposal.AdaptSteps'] = 100
options['Kernel1.MyProposal.AdaptStart'] = 500
# create the MCMC sampler
mcmc = msa.SingleChainMCMC(options, problem_fine)
samps = mcmc.Run([trueCoeff])
print(mcmc.GetSamples().Mean())
print(mcmc.GetQOIs().Mean())
print(samps.Mean())
# MLMCMC
mioptions = dict()
mioptions['NumSamples_0'] = nmcmc * 10
mioptions['NumSamples_1'] = nmcmc
mioptions['Subsampling'] = 10
mioptions['Proposal.Method'] = 'AMProposal'
mioptions['Proposal.InitialVariance'] = 0.1
mioptions['Proposal.AdaptSteps'] = 100
mioptions['Proposal.AdaptStart'] = 500
mimcmc = msa.MIMCMC(mioptions, trueCoeff, [problem_coarse, problem_fine])
mimcmc.Run([trueCoeff])
print(mimcmc.MeanParam())
print(mimcmc.MeanQOI())
print(cacheModPiece.HitRatio())
# Let's dive into the MIMCMC data structures to do some post processing
# Iterate over all multiindices we have
mlindices = mimcmc.GetIndices()
for i in range(0,mlindices.Size()):
# Get the MIMCMCBox representing the telescoping sum component associated with this index
box = mimcmc.GetMIMCMCBox(mlindices.at(i))
boxIndices = box.GetBoxIndices()
print(mlindices.at(i).GetVector())
# Print the contribution to the multilevel telescoping sum of this box
print(box.MeanParam())
fig=plt.figure()
ax=fig.add_axes([0,0,1,1])
plt.xlim(-5, 8)
plt.ylim(-5, 8)
# Plot all samples of every MCMC chain of this box
for i in range(0,boxIndices.Size()):
boxIndex = boxIndices.at(i)
boxChain = box.GetChain(boxIndex)
samplesMat = boxChain.GetSamples().AsMatrix()
if boxIndex.GetValue(0) == 0:
ax.scatter(samplesMat[0,:], samplesMat[1,:], color='b', alpha=0.3, label="Coarse Samples")
else:
ax.scatter(samplesMat[0,:], samplesMat[1,:], color='r', alpha=0.3, label="Fine Samples")
# Let's access an individual sample
samples = boxChain.GetSamples()
print("Sample 10")
print(samples[10].state)
# and if it had a coarse proposal, get that as well:
if samples[10].HasMeta("coarseSample"):
print ("from coarse sample:")
print(samples[10].GetMetaSamplingState("coarseSample").state)
ax.legend()
# Plot all mean values of every MCMC chain of this box
for i in range(0,boxIndices.Size()):
boxIndex = boxIndices.at(i)
boxChain = box.GetChain(boxIndex)
mean = boxChain.GetSamples().Mean()
plt.axvline(x=mean[0], color='b' if boxIndex.GetValue(0) == 0 else 'r');
plt.axhline(y=mean[1], color='b' if boxIndex.GetValue(0) == 0 else 'r');
MUQ is on Slack! - Join our slack workspace to connect with other users, get help, and discuss new features.
Test Status
Acknowledgments
This material is based upon work supported by the National Science Foundation under Grant No. 1550487.
Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.
This material is based upon work supported by the US Department of Energy, Office of Advanced Scientific Computing Research, SciDAC (Scientific Discovery through Advanced Computing) program under awards DE-SC0007099 and DE-SC0021226, for the QUEST and FASTMath SciDAC Institutes.