I am currently dealing with exponential decay (of muons). The PDF is the following: https://wikimedia.org/api/rest_v1/media/math/render/svg/fdacadb747bac9932ffd02ee35cf3cbd0bc0b4c8

And I am told that mouns have a theoretical mean lifetime $\tau=2$. Also lifetime below $tmin=1$ or above $tmax=10$ is not possible, because of device issues.

Initially I was asked to write the pdf for the exp. decay. That I did, by integrating within the interval, normalizing and finding the normalization constant. This was my code (python):

def pdf_exp_d(t,tau,tmin,tnax):
    C = 1/(tau*(np.exp(-t_min/tau)-np.exp(-t_max/tau)))
    return C*np.exp(-t/tau)

Then I had to generate random numbers (50 of them) within the $[t_{min};t_{max}]$. As a help I was given the following line:

## def exponentialDecayPDF(t, tau= tau, tmin=tmin, tmax=tmax):

Now, I know how to produce random numbers with exp distribution. I can do that by directly using the np.random.exponential, or I can use inverse sampling methode. But the problem is, that I cannot produce random numbers with this distribution within an interval. Is there some sort of hint that I can get,as to how I should proceed?

Thank you


1 Answer 1


The inverse transform sampling method is generally used to allow sampling uniform random variables and then transforming them such that the result has the suitable distribution. To do this you need an inverse CDF for your function:

import sympy as sp
Pr, s, t, tau, t_min, t_max = sp.symbols('Pr, s, t, tau, t_min, t_max', real=True, positive=True)
C = 1 / (tau * (sp.exp(-t_min / tau) - sp.exp(-t_max / tau)))
pdf = C*sp.exp(-t/tau)
cdf = sp.integrate(pdf, (t, t_min, s))
cdf = cdf.args[0][0] # Select the relevant solution
icdf = sp.solve(cdf-Pr, s) # Solve for the inverse

$F^{-1}(\mathrm{Pr}) = \displaystyle t_{max} + t_{min} + \tau \log{\left(\frac{1}{- \mathrm{Pr} \; e^{\frac{t_{max}}{\tau}} + \mathrm{Pr} \; e^{\frac{t_{min}}{\tau}} + e^{\frac{t_{max}}{\tau}}} \right)}$

Then you can generate a sample as

import matplotlib.pyplot as plt
import numpy as np
from numpy import exp, log
def inv_cdf(Pr, t_min, t_max, tau):
    return t_max + t_min + tau*log(1/(-Pr*exp(t_max/tau) + Pr*exp(t_min/tau) + exp(t_max/tau)))                            
U_vals = np.random.rand(1000)
T_vals = inv_cdf(U_vals, 1,10,2)


Histogram of resulting samples


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.