Generate copula-correlated samples with specified marginals in Python

Question

I have N random variables (X1,...,XN) each of which is distributed over a specific marginal (normal, log-normal, Poisson...) and I want to generate a sample of p joint realizations of these variables Xi, given that the variables are correlated with a given Copula, using Python 3. I know that R is a better option but i want to do it in Python.

Following this method I managed to do so with a Gaussian Copula. Now I want to adapt the method to use a Archimedean Copula (Gumbel, Frank...) or a Student Copula. Ath the beginning of the Gaussian copula method, you draw a sample of p realizations from a multivariate normal distribution. To adapt this to another copula, for instance a bivariate Gumbel, my idea is to draw a sample from the joint distribution of a bivariate Gumbel, but I am not sure on how to implement this.

I have tried using several Python 3 packages : copulae, copula and copulas all provide the noption to fit a particular copula to a dataset but do not allow to draw a random sample from a given copula.

Can you provide some algorithmic insight on how to draw multivariate random samples from a given Copula with uniform marginals?

Answer 1

Please check this page "Create a composed distribution". I think this is what you're looking for

For example, if you have 2 distributions x1: Uniform on [1, 3] and x2: Normal(0,2) and if you know the dependence structure copula, you can build the multidimensional distribution X = (x1, x2) very easily.

import openturns as ot

x1 = ot.Uniform(1, 3)
x2 = ot.Normal(0, 2)
copula = ot.IndependentCopula()

X = ot.ComposedDistribution([x1, x2], copula)

X.getSample(5) will return a sample of size = 5:

>>>   [ X0         X1         ]
0 : [  1.87016    0.802719  ]
1 : [  1.72333    2.73565   ]
2 : [  1.00422    2.00869   ]
3 : [  1.47887    1.4831    ]
4 : [  1.51031   -0.0872247 ]

You can view the cloud in 2D

import matplotlib.pyplot as plt

sample = dist.getSample(1000)
plt.scatter(sample[:, 0], sample[:, 1], s=2)

If you choose copula = ot.ClaytonCopula(2) the result will be:

with GumbelCopula(2):