.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/inference/mle/plot_complex_probability_model.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_inference_mle_plot_complex_probability_model.py: Complex probability distribution model ============================================== Here we define a bivariate probability model, with a dependence structure defined using a gumbel copula. The goal of inference is to learn the paremeters of the Gaussian marginals and the copula parameter, i.e., the model has 5 unknown parameters. .. GENERATED FROM PYTHON SOURCE LINES 14-15 Initially we have to import the necessary modules. .. GENERATED FROM PYTHON SOURCE LINES 18-26 .. code-block:: default import matplotlib.pyplot as plt from UQpy.inference import DistributionModel, MLE from UQpy.distributions import Normal from UQpy.inference import MinimizeOptimizer from UQpy.distributions import JointIndependent, JointCopula, Gumbel from UQpy.sampling import ImportanceSampling .. GENERATED FROM PYTHON SOURCE LINES 27-29 First data is generated from a true model. A distribution with copulas does not possess a fit method, thus sampling is performed using importance sampling/resampling. .. GENERATED FROM PYTHON SOURCE LINES 32-52 .. code-block:: default # dist_true exhibits dependence between the two dimensions, defined using a gumbel copula dist_true = JointCopula(marginals=[Normal(), Normal()], copula=Gumbel(theta=2.)) # generate data using importance sampling: sample from a bivariate gaussian without copula, then weight samples u = ImportanceSampling(proposal = JointIndependent(marginals=[Normal(), Normal()]), log_pdf_target = dist_true.log_pdf, nsamples=500) print(u.samples.shape) print(u.weights.shape) # Resample to obtain 5,000 data points u.resample(nsamples=5000) data_2 = u.unweighted_samples print('Shape of data: {}'.format(data_2.shape)) fig, ax = plt.subplots() ax.scatter(data_2[:, 0], data_2[:, 1], alpha=0.2) ax.set_title('Data points from true bivariate normal with gumbel dependency structure') plt.show() .. image-sg:: /auto_examples/inference/mle/images/sphx_glr_plot_complex_probability_model_001.png :alt: Data points from true bivariate normal with gumbel dependency structure :srcset: /auto_examples/inference/mle/images/sphx_glr_plot_complex_probability_model_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (500, 2) (500,) Shape of data: (5000, 2) .. GENERATED FROM PYTHON SOURCE LINES 53-57 To define a model for inference, the user must create a custom file, here bivariate_normal_gumbel.py, to compute the log_pdf of the distribution, given a bivariate data matrix and a parameter vector of length 5. Note that for any probability model that is not one of the simple univariate pdfs supported by UQpy, such a custom file will be necessary. .. GENERATED FROM PYTHON SOURCE LINES 60-67 .. code-block:: default d_guess = JointCopula(marginals=[Normal(loc=None, scale=None), Normal(loc=None, scale=None)], copula=Gumbel(theta=None)) print(d_guess.get_parameters()) candidate_model = DistributionModel(n_parameters=5, distributions=d_guess) print(candidate_model.list_params) .. rst-class:: sphx-glr-script-out .. code-block:: none {'loc_0': None, 'scale_0': None, 'loc_1': None, 'scale_1': None, 'theta_c': None} ['loc_0', 'scale_0', 'loc_1', 'scale_1', 'theta_c'] .. GENERATED FROM PYTHON SOURCE LINES 68-71 When calling MLEstimation, the function minimize from the scipy.optimize package is used by default. The user can define bounds for the optimization, a seed, the algorithm to be used, and set the algorithm to perform several optimization iterations, starting at a different random seed every time. .. GENERATED FROM PYTHON SOURCE LINES 74-87 .. code-block:: default optimizer = MinimizeOptimizer(bounds=[[-5, 5], [0, 10], [-5, 5], [0, 10], [1.1, 4]], method="SLSQP") ml_estimator = MLE(inference_model=candidate_model, data=data_2, optimizer=optimizer) ml_estimator = MLE(inference_model=candidate_model, data=data_2, optimizer=optimizer, initial_parameters=[1., 1., 1., 1., 4.]) print('ML estimates of the mean={0:.3f} and std. dev={1:.3f} of 1st marginal (true: 0.0, 1.0)'. format(ml_estimator.mle[0], ml_estimator.mle[1])) print('ML estimates of the mean={0:.3f} and std. dev={1:.3f} of 2nd marginal (true: 0.0, 1.0)'. format(ml_estimator.mle[2], ml_estimator.mle[3])) print('ML estimates of the copula parameter={0:.3f} (true: 2.0)'.format(ml_estimator.mle[4])) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/utilities/MinimizeOptimizer.py:29: OptimizeWarning: Unknown solver options: catol return minimize(function, initial_guess, args=args, /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:63: RuntimeWarning: divide by zero encountered in log c = exp(-(((-log(u)) ** theta + (-log(v)) ** theta) ** (1 / theta))) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:65: RuntimeWarning: invalid value encountered in divide pdf_val = (c * 1 / u * 1 / v /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:66: RuntimeWarning: divide by zero encountered in log * ((-log(u)) ** theta + (-log(v)) ** theta) ** (-2 + 2 / theta) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:67: RuntimeWarning: divide by zero encountered in log * (log(u) * log(v)) ** (theta - 1) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:67: RuntimeWarning: invalid value encountered in multiply * (log(u) * log(v)) ** (theta - 1) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:68: RuntimeWarning: divide by zero encountered in log * (1 + (theta - 1) * ((-log(u)) ** theta + (-log(v)) ** theta) ** (-1 / theta))) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:66: RuntimeWarning: divide by zero encountered in power * ((-log(u)) ** theta + (-log(v)) ** theta) ** (-2 + 2 / theta) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:65: RuntimeWarning: invalid value encountered in multiply pdf_val = (c * 1 / u * 1 / v /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/copulas/Gumbel.py:68: RuntimeWarning: divide by zero encountered in power * (1 + (theta - 1) * ((-log(u)) ** theta + (-log(v)) ** theta) ** (-1 / theta))) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/scipy/stats/_distn_infrastructure.py:2029: RuntimeWarning: divide by zero encountered in divide x = np.asarray((x - loc)/scale, dtype=dtyp) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/scipy/stats/_distn_infrastructure.py:2071: RuntimeWarning: divide by zero encountered in divide x = np.asarray((x - loc)/scale, dtype=dtyp) /home/docs/checkouts/readthedocs.org/user_builds/uqpyproject/envs/stable/lib/python3.9/site-packages/UQpy/distributions/collection/JointCopula.py:121: RuntimeWarning: divide by zero encountered in log return np.log(c_) + logpdf_val ML estimates of the mean=0.152 and std. dev=1.038 of 1st marginal (true: 0.0, 1.0) ML estimates of the mean=0.165 and std. dev=1.066 of 2nd marginal (true: 0.0, 1.0) ML estimates of the copula parameter=2.087 (true: 2.0) .. GENERATED FROM PYTHON SOURCE LINES 88-89 Again, some known parameters can be fixed during learning. .. GENERATED FROM PYTHON SOURCE LINES 92-104 .. code-block:: default d_guess = JointCopula(marginals=[Normal(loc=None, scale=None), Normal(loc=0., scale=1.)], copula=Gumbel(theta=None)) candidate_model = DistributionModel(n_parameters=3, distributions=d_guess) optimizer = MinimizeOptimizer(bounds=[[-5, 5], [0, 10], [1.1, 4]], method="SLSQP") ml_estimator = MLE(inference_model=candidate_model, data=data_2, optimizer=optimizer, initial_parameters=[1., 1., 4.]) print('ML estimates of the mean={0:.3f} and std. dev={1:.3f} of 1st marginal (true: 0.0, 1.0)'. format(ml_estimator.mle[0], ml_estimator.mle[1])) print('ML estimates of the copula parameter={0:.3f} (true: 2.0)'.format(ml_estimator.mle[2])) .. rst-class:: sphx-glr-script-out .. code-block:: none ML estimates of the mean=0.025 and std. dev=0.983 of 1st marginal (true: 0.0, 1.0) ML estimates of the copula parameter=1.923 (true: 2.0) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.359 seconds) .. _sphx_glr_download_auto_examples_inference_mle_plot_complex_probability_model.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/SURGroup/UQpy/master?urlpath=lab/tree/notebooks/auto_examples/inference/mle/plot_complex_probability_model.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_complex_probability_model.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_complex_probability_model.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_