How to avoid infinities in the following example?

Question

The task

Consider a decaying particle traveling along some axis $x$, with the differential probability to decay given by the Poisson distribution: $$ \frac{dP_{\text{decay}}}{dx} = \frac{\exp[-x/l_{\text{decay}}]}{l_{\text{decay}}} $$ I would like to generate its decay points inside the given interval $x_{\text{min}},x_{\text{max}}$.

The problem

My implementation of this task encounters numeric instabilities if $l_{\text{decay}}\ll x_{\text{min}}$.

My code

xmin = 14;
xmax = 34;
(*CDF, inverse CDF, and the values of CDF corresponding to xmin,xmax*)
CDFpoisson[xdecay_, ldecay_] = 
 Integrate[Exp[-x/ldecay]/ldecay, {x, 0, xdecay}]
invCDFpoisson[u_, ldecay_] = 
  x /. Solve[CDFpoisson[x, ldecay] == u, x][[1]] /. {C[1] -> 0};
uxmin[ldecay_] = 
  u /. Solve[invCDFpoisson[u, ldecay] == xmin, u, Reals][[1]];
uxmax[ldecay_] = 
  u /. Solve[invCDFpoisson[u, ldecay] == xmax, u, Reals][[1]];
(*Random values of u corresponding to the interval, and random decay points*)
DecayPointsData[ldecay_] := Block[{},
TestPoints = RandomReal[{uxmin[ldecay], uxmax[ldecay]}, 10^6];
Table[invCDFpoisson[TestPoints[[i]], ldecaytest], {i, 1, 
Length[TestPoints], 1}]
]

The code works properly for large values of ldecay:

ldecaytest=0.7;
TabDecayPoints = DecayPointsData[ldecaytest];
integralVal = Integrate[Exp[-x/ldecaytest]/ldecaytest, {x, 14, 34}];
Show[LogLogPlot[Exp[-x/ldecaytest]/(
  ldecaytest*integralVal), {x, xmin, xmax}, 
  PlotRange -> {{xmin, xmax}, {3.5*10^-6, 3}}, Frame -> True, 
  ImageSize -> "Large"], 
 Histogram[TabDecayPoints, 100, "ProbabilityDensity", 
  ScalingFunctions -> {"Log", "Log"}]]

However, it breaks down at smaller values:

ldecaytest = 0.4;
TabDecayPoints = DecayPointsData[ldecaytest];

Infinite expression 1/0. encountered.

The reason is that uxmin, uxmax are approximated by 1 for large ratio xmin/ldecaytest, which results in infinite values of invCDFpoisson.

Question

Could you please tell me how to avoid these infinities without a significant reduction in speed?

Answer 1

First of all: your code can be simplified significantly

Clear[dist]
dist[ldecay_] = ExponentialDistribution[1/ldecay] (* What you call a Poisson dist is actually an exponential distribution *)

CDF[dist[ldecay], x] (* check the CDF *) 

This should compute the values you want. Note that you need to use arbitrary precision to calculate these values:

{xmin, xmax} = {14, 34};
DecayPointsData[ldecay_, n_] := Block[{TestPoints},
  TestPoints = RandomReal[CDF[dist[ldecay], {xmin, xmax}], n, WorkingPrecision -> Precision[ldecay]];
  InverseCDF[dist[ldecay], TestPoints]
];

DecayPointsData[0.4`20, 10^3]

Edit

Note that CDF[dist[0.4], {xmin, xmax}] gives numbers very close to 1. These are difficult to represent because you can't use powers of 10 to represent them. Instead, you should consider using the survival function:

SurvivalFunction[dist[0.4], {xmin, xmax}]

As you can see, this can easily be computed with machine precision.