Example Multiindex MCMC
Defines a two dimensional hierarchy of simple Gaussian models and applies Multiindex 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
cov = np.array([[1.0,0.3],[0.2,1.5]]) # covariance
m = np.array([.8, 2.4]) # mean
targetDensity = mm.Gaussian(m, cov * 2.0).AsDensity()
problem_coarse12 = msa.SamplingProblem(targetDensity)
m = np.array([.8, 2.1]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov * 1.3).AsDensity()
problem_coarse1 = msa.SamplingProblem(targetDensity)
m = np.array([0.95, 2.4]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov * 1.5).AsDensity()
problem_coarse2 = msa.SamplingProblem(targetDensity)
m = np.array([1.0, 2.0]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov).AsDensity()
problem_fine = msa.SamplingProblem(targetDensity)
nmcmc = 10000
trueCoeff = m
# 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.01
options['Kernel1.MyProposal.AdaptSteps'] = 50
options['Kernel1.MyProposal.AdaptStart'] = 100
# create the MCMC sampler
mcmc = msa.SingleChainMCMC(options, problem_fine)
samps = mcmc.Run([trueCoeff])
print(samps.Mean())
nmimcmc = 1000
#MIMCMC
mioptions = dict()
mioptions['NumSamples_0_0'] = nmimcmc * 100 # Number of samples per level
mioptions['NumSamples_1_0'] = nmimcmc * 10
mioptions['NumSamples_0_1'] = nmimcmc * 10
mioptions['NumSamples_1_1'] = nmimcmc
mioptions['Subsampling'] = 10
mioptions['Proposal.Method'] = 'AMProposal'
mioptions['Proposal.InitialVariance'] = 0.01
mioptions['Proposal.AdaptSteps'] = 50
mioptions['Proposal.AdaptStart'] = 100
# Let's define a set of multiindices for our models
multiindexset = mu.MultiIndexFactory.CreateFullTensor(orders=[1,1])
# And put our models in a list, ordered like the multiindices
models = []
for i in range(0,multiindexset.Size()):
print(multiindexset.at(i).GetVector())
if (multiindexset.at(i).GetVector() == [0, 0]).all():
models.append(problem_coarse12)
elif (multiindexset.at(i).GetVector() == [1, 0]).all():
models.append(problem_coarse1)
elif (multiindexset.at(i).GetVector() == [0, 1]).all():
models.append(problem_coarse2)
elif (multiindexset.at(i).GetVector() == [1, 1]).all():
models.append(problem_fine)
else:
print ("Model not defined for index!")
# Now, plug models into MIMCMC
mimcmc = msa.MIMCMC(mioptions, trueCoeff, multiindexset, models)
mimcmc.Run([trueCoeff])
print(mimcmc.MeanParam())
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
cov = np.array([[1.0,0.3],[0.2,1.5]]) # covariance
m = np.array([.8, 2.4]) # mean
targetDensity = mm.Gaussian(m, cov * 2.0).AsDensity()
problem_coarse12 = msa.SamplingProblem(targetDensity)
m = np.array([.8, 2.1]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov * 1.3).AsDensity()
problem_coarse1 = msa.SamplingProblem(targetDensity)
m = np.array([0.95, 2.4]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov * 1.5).AsDensity()
problem_coarse2 = msa.SamplingProblem(targetDensity)
m = np.array([1.0, 2.0]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov).AsDensity()
problem_fine = msa.SamplingProblem(targetDensity)
nmcmc = 10000
trueCoeff = m
# 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.01
options['Kernel1.MyProposal.AdaptSteps'] = 50
options['Kernel1.MyProposal.AdaptStart'] = 100
# create the MCMC sampler
mcmc = msa.SingleChainMCMC(options, problem_fine)
samps = mcmc.Run([trueCoeff])
print(samps.Mean())
nmimcmc = 1000
#MIMCMC
mioptions = dict()
mioptions['NumSamples_0_0'] = nmimcmc * 100 # Number of samples per level
mioptions['NumSamples_1_0'] = nmimcmc * 10
mioptions['NumSamples_0_1'] = nmimcmc * 10
mioptions['NumSamples_1_1'] = nmimcmc
mioptions['Subsampling'] = 10
mioptions['Proposal.Method'] = 'AMProposal'
mioptions['Proposal.InitialVariance'] = 0.01
mioptions['Proposal.AdaptSteps'] = 50
mioptions['Proposal.AdaptStart'] = 100
# Let's define a set of multiindices for our models
multiindexset = mu.MultiIndexFactory.CreateFullTensor(orders=[1,1])
# And put our models in a list, ordered like the multiindices
models = []
for i in range(0,multiindexset.Size()):
print(multiindexset.at(i).GetVector())
if (multiindexset.at(i).GetVector() == [0, 0]).all():
models.append(problem_coarse12)
elif (multiindexset.at(i).GetVector() == [1, 0]).all():
models.append(problem_coarse1)
elif (multiindexset.at(i).GetVector() == [0, 1]).all():
models.append(problem_coarse2)
elif (multiindexset.at(i).GetVector() == [1, 1]).all():
models.append(problem_fine)
else:
print ("Model not defined for index!")
# Now, plug models into MIMCMC
mimcmc = msa.MIMCMC(mioptions, trueCoeff, multiindexset, models)
mimcmc.Run([trueCoeff])
print(mimcmc.MeanParam())
Defines a two dimensional hierarchy of simple Gaussian models and applies Multiindex 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
cov = np.array([[1.0,0.3],[0.2,1.5]]) # covariance m = np.array([.8, 2.4]) # mean targetDensity = mm.Gaussian(m, cov * 2.0).AsDensity() problem_coarse12 = msa.SamplingProblem(targetDensity) m = np.array([.8, 2.1]) # mean #cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance targetDensity = mm.Gaussian(m, cov * 1.3).AsDensity() problem_coarse1 = msa.SamplingProblem(targetDensity) m = np.array([0.95, 2.4]) # mean #cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance targetDensity = mm.Gaussian(m, cov * 1.5).AsDensity() problem_coarse2 = msa.SamplingProblem(targetDensity) m = np.array([1.0, 2.0]) # mean #cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance targetDensity = mm.Gaussian(m, cov).AsDensity() problem_fine = msa.SamplingProblem(targetDensity)
nmcmc = 10000 trueCoeff = m # 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.01 options['Kernel1.MyProposal.AdaptSteps'] = 50 options['Kernel1.MyProposal.AdaptStart'] = 100 # create the MCMC sampler mcmc = msa.SingleChainMCMC(options, problem_fine) samps = mcmc.Run([trueCoeff]) print(samps.Mean())
nmimcmc = 1000
#MIMCMC
mioptions = dict()
mioptions['NumSamples_0_0'] = nmimcmc * 100 # Number of samples per level
mioptions['NumSamples_1_0'] = nmimcmc * 10
mioptions['NumSamples_0_1'] = nmimcmc * 10
mioptions['NumSamples_1_1'] = nmimcmc
mioptions['Subsampling'] = 10
mioptions['Proposal.Method'] = 'AMProposal'
mioptions['Proposal.InitialVariance'] = 0.01
mioptions['Proposal.AdaptSteps'] = 50
mioptions['Proposal.AdaptStart'] = 100
# Let's define a set of multiindices for our models
multiindexset = mu.MultiIndexFactory.CreateFullTensor(orders=[1,1])
# And put our models in a list, ordered like the multiindices
models = []
for i in range(0,multiindexset.Size()):
print(multiindexset.at(i).GetVector())
if (multiindexset.at(i).GetVector() == [0, 0]).all():
models.append(problem_coarse12)
elif (multiindexset.at(i).GetVector() == [1, 0]).all():
models.append(problem_coarse1)
elif (multiindexset.at(i).GetVector() == [0, 1]).all():
models.append(problem_coarse2)
elif (multiindexset.at(i).GetVector() == [1, 1]).all():
models.append(problem_fine)
else:
print ("Model not defined for index!")
# Now, plug models into MIMCMC
mimcmc = msa.MIMCMC(mioptions, trueCoeff, multiindexset, models)
mimcmc.Run([trueCoeff])
print(mimcmc.MeanParam())
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
cov = np.array([[1.0,0.3],[0.2,1.5]]) # covariance
m = np.array([.8, 2.4]) # mean
targetDensity = mm.Gaussian(m, cov * 2.0).AsDensity()
problem_coarse12 = msa.SamplingProblem(targetDensity)
m = np.array([.8, 2.1]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov * 1.3).AsDensity()
problem_coarse1 = msa.SamplingProblem(targetDensity)
m = np.array([0.95, 2.4]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov * 1.5).AsDensity()
problem_coarse2 = msa.SamplingProblem(targetDensity)
m = np.array([1.0, 2.0]) # mean
#cov = np.array([[1.0,0.8],[0.8,1.5]]) # covariance
targetDensity = mm.Gaussian(m, cov).AsDensity()
problem_fine = msa.SamplingProblem(targetDensity)
nmcmc = 10000
trueCoeff = m
# 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.01
options['Kernel1.MyProposal.AdaptSteps'] = 50
options['Kernel1.MyProposal.AdaptStart'] = 100
# create the MCMC sampler
mcmc = msa.SingleChainMCMC(options, problem_fine)
samps = mcmc.Run([trueCoeff])
print(samps.Mean())
nmimcmc = 1000
#MIMCMC
mioptions = dict()
mioptions['NumSamples_0_0'] = nmimcmc * 100 # Number of samples per level
mioptions['NumSamples_1_0'] = nmimcmc * 10
mioptions['NumSamples_0_1'] = nmimcmc * 10
mioptions['NumSamples_1_1'] = nmimcmc
mioptions['Subsampling'] = 10
mioptions['Proposal.Method'] = 'AMProposal'
mioptions['Proposal.InitialVariance'] = 0.01
mioptions['Proposal.AdaptSteps'] = 50
mioptions['Proposal.AdaptStart'] = 100
# Let's define a set of multiindices for our models
multiindexset = mu.MultiIndexFactory.CreateFullTensor(orders=[1,1])
# And put our models in a list, ordered like the multiindices
models = []
for i in range(0,multiindexset.Size()):
print(multiindexset.at(i).GetVector())
if (multiindexset.at(i).GetVector() == [0, 0]).all():
models.append(problem_coarse12)
elif (multiindexset.at(i).GetVector() == [1, 0]).all():
models.append(problem_coarse1)
elif (multiindexset.at(i).GetVector() == [0, 1]).all():
models.append(problem_coarse2)
elif (multiindexset.at(i).GetVector() == [1, 1]).all():
models.append(problem_fine)
else:
print ("Model not defined for index!")
# Now, plug models into MIMCMC
mimcmc = msa.MIMCMC(mioptions, trueCoeff, multiindexset, models)
mimcmc.Run([trueCoeff])
print(mimcmc.MeanParam())
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.