Monte Carlo sampling with UQpy

We’ll be using UQpy’s Monte Carlo sampling functionalities. We also use Matplotlib to display results graphically.

Additionally, this demonstration opts to use Numpy’s random state management to ensure that results are reproducible between notebook runs.

from UQpy.sampling import MonteCarloSampling
import matplotlib.pyplot as plt

from numpy.random import RandomState

Step-by-step: continuous univariate distribution

First, we import UQpy’s normal distribution class.

from UQpy.distributions import Normal

We’ll start by constructing two identical standard normal distributions normal1 and normal2

normal1 = normal2 = Normal()

Next, we’ll construct a MonteCarloSampling object mc to generate random samples following those distributions. Here, we specify an optional initial number of samples, nsamples to be generated at the object’s construction. For teh purposes of this demonstration, we also supply a random seed random_state.

We access the generated samples via the samples attribute.

mc = MonteCarloSampling(distributions=[normal1, normal2],
                        nsamples=5,
                        random_state=RandomState(123))

mc.samples

To generate more samples on mc after construction, we call mc.run and once again specify nsamples.

mc.run(nsamples=2, random_state=RandomState(23))

mc.samples

We can transform the samples onto the unit hypercube via applying the probability integral transformation on the samples to yield similar samples from the uniform distribution. We call mc.transform_u01, from which results are stored in the samplesU01 attribute.

mc.transform_u01()

mc.samplesU01

We can visualize the (untransformed) samples by plotting them on axes of each distribution’s range.

fig, ax = plt.subplots()
plt.title('Samples')

plt.scatter(x=mc.samples[:, 0],
            y=mc.samples[:, 1],
            marker='o')

plt.setp(ax, xlim=(-1.7, 1.7), ylim=(-2.6, 2.6))
ax.yaxis.grid(True)
ax.xaxis.grid(True)

As well, we can visualize each distribution’s sample densities via histograms.

fig, ax = plt.subplots(1, 2)

for i in (0, 1):
    ax[i].set_title('Distribution ' + str(i + 1))
    ax[i].hist(mc.samples[:, i])
    ax[i].yaxis.grid(True)
    ax[i].xaxis.grid(True)

plt.setp(ax, xlim=(-3, 3), ylim=(0, 2));

Additional Examples

Continuous multivariate distribution

We’ll use UQpy’s multivariate normal class.

from UQpy.distributions import MultivariateNormal

And we construct a multivariate normal distribution mvnormal specifying parameters mean with a vector of mean values and cov with a covariance matrix.

mvnormal = MultivariateNormal(mean=[1, 2],
                              cov=[[4, -0.9],
                                   [-0.9, 1]])

With this distribution, we construct a MonteCarloSampling object mvmc and generate five samples on construction.

mvmc = MonteCarloSampling(distributions=mvnormal,
                          nsamples=5,
                          random_state=RandomState(456))

mvmc.samples

Mixing a multivariate and a univariate continuous distribution

Here, we use one of our normal distributions and our multivariate normal distribution. Notice how each distribution has its own bundle (array) of samples per run of sampling—even when that bundle contains a single value.

mixedmc = MonteCarloSampling(distributions=[normal1, mvnormal],
                             nsamples=5,
                             random_state=RandomState(789))

mixedmc.samples

Mixing a continuous and a discrete distribution

We’ll use UQpy’s binomial distribution class for our discrete distribution.

from UQpy.distributions import Binomial

With that, we’ll construct a binomial distribution binomial with five trials n and a 40% probability p of success per trial.

binomial = Binomial(n=5, p=0.4)

And we construct a MonteCarloSampling object cdmv with five initial samples using our binomial distribution and one of our normal distributions.

cdmv = MonteCarloSampling(distributions=[binomial, normal1],
                          nsamples=5,
                          random_state=RandomState(333))

cdmv.samples