How do you generate random variables according to a given continuous probability density function?

Question

I'm trying to generate random variables according to three different probability density functions. Two of the functions are scaled normal distributions, µ = 260/3, σ = 100/3, scaled by 1.53666351546 and µ = 854/9, σ = 100/3, scaled by 1.78979465778 respectively, and the third is the sum of eleven different normal distributions and two half-normal distributions.

I tried messing with transforming the numpy.random.randn function for the first two distributions, but I couldn't get the maths to line up. I'm yet to try anything for the third distribution, because I don't even know where to start.

Using the code created by sato in this question, I was able to modify it with Peter O.'s answer to get the following code for my first distribution:

from decimal import Decimal
from scipy import stats
import numpy

Q = Decimal(0.0460999054638)/((Decimal(2)*Decimal(numpy.pi))**Decimal(0.5))

class Distribution_B(stats.rv_continuous):
    def _pdf(self, x):
        return (Q*Decimal(numpy.exp(-(((Decimal(260)-Decimal(3)*Decimal(x))**Decimal(5))/Decimal(20000)))))

distribution = Distribution_B(a = 0, b = 100)
distribution.rvs()

However, when I ran the code, I encountered the following error:

Warning (from warnings module):
  File "[numpy function_base library]", line 2411
    outputs = ufunc(*inputs)
RuntimeWarning: invalid value encountered in _cdf_single (vectorized)

Edit: I ran the code without changing anything except variable names, and it seems that the error only happens sometimes.

Answer 1

from decimal import Decimal
from scipy import stats
import numpy as np

Q = Decimal(0.0460999054638) / ((Decimal(2) * Decimal(np.pi)) ** Decimal(0.5))

class Distribution_B(stats.rv_continuous):
    def _pdf(self, x):
        x = float(x)  # Convert x to float
        return Q * Decimal(np.exp(-(((Decimal(260) / Decimal(3) - Decimal(x)) ** Decimal(2)) / (2 * (Decimal(100) / Decimal(3)) ** Decimal(2)))))

distribution = Distribution_B(a=0, b=100)
samples = distribution.rvs(size=1000)

Remember to adjust the size on the next line