Generate bivariate random numbers from a joint probability distribution in Python

Question

I have two distributions over two parameters h and t. h is Weibull distributed while t is conditioned on h and it is log-normal distributed:

f_H = np.exp(-(h / alpha) ** beta) * (beta / alpha) * (h / alpha) ** (beta - 1)
f_TIH = np.exp(-(np.log(t) - mu_h) ** 2.0 / (2.0 * sigma_h ** 2)) / (t * sigma_h * np.sqrt(2.0 * np.pi))

where:

mu_h = a0 + a1 * h ** a2
sigma_h = b0 + b1 * np.exp(b2 * h)

and:

a0 = 0.7
a1 = 0.282
a2 = 0.167
b0 = 0.07
b1 = 0.3449
b2 = -0.2073
alpha = 1.76
beta = 1.59

The joint PDF for h and t is then given as:

f_joint = f_H * f_TIH

My question is how can I sample random values for h and t from the joint PDF?

Answer 1

Since you have clear conditioning, you can generate $h$ with a numpy.random.weibull: $$ h = \alpha \ \mathtt{weilbull}(\beta) $$

and then generate $t$ from log-normal distribution numpy.random.lognormal: $$ t = \mathtt{lognormal}(\mu_h, \sigma_h) $$

As a quick check, you can generate several millions of points and check that density plot looks like a joint PDF.