Fast sampling from Fermi-Dirac distribution under specific conditions

This question is a continuation of the following one.

Let us assume a Fermi-Dirac (FD) distribution (times $$p^{2}$$):

AvalFD = Integrate[x^2/(Exp[x/T] + 1), {x, 0, Infinity}, Assumptions -> T > 0];
DistrFD[T_,x_]=1/AvalFD*x^2/(Exp[x/T]+1);


I need to make a fast point sampler from it. There are a few complications:

• I need to sample a very small number of points (in practice, 1 or 2) but many times (so there should be no lower time threshold for calling the sampler).
• The temperature is a dynamic parameter - it changes from time to time (so it may be complicated to precompute the inverse CDF).
• I need to sample inside a compiled code (so the sampling code must be compilable).

This is my approach:

FastFDsampler =
Compile[{{npoints, _Integer}, {T, _Real}},
Module[{pts, weights},
pts = RandomReal[{0, 10 T}, 50 npoints];
weights = pts^2/(Exp[pts/T] + 1);
RandomSample[weights -> pts, npoints]
]
, CompilationTarget -> "C", RuntimeOptions -> "Speed",
RuntimeAttributes -> {Listable}];
FastFDsampler[1, 3]; // AbsoluteTiming
sss = FastFDsampler[10000, 3]; // AbsoluteTiming
Show[Histogram[sss, 30, "ProbabilityDensity"],
Plot[DistrFD[3, x], {x, 0, 30}]]


{0.0000379,Null}

{0.0112687,Null}

It first samples points uniformly within some range, then computes the weights from the distribution DistrFD, and then resamples. The method is not very good, as one must first generate a much larger amount of points than the desired number of final points (the difference between the uniform distribution and FD is huge); also, there is an upper limit on the maximal value of the generated value (10T in the example) imposed by hand.

As an alternative, I came up with the following approach:

FastFDsampler2 =
Compile[{{npoints, _Integer}, {T, _Real}},
Module[{ptsBoltzmann, weightsRatio},
(*Sampling points according to the Boltzmann distribution \
E^2Exp[-E/T] - to account for the correct high-energy tail*)
ptsBoltzmann = -T*Total[Log[#]] & /@
RandomVariate[
UniformDistribution[{0, 1}], {5 npoints, 3}];
(*Ratio of the weigths - Fermi-Dirac over Boltzmann*)
weightsRatio = Exp[ptsBoltzmann/T]/(Exp[ptsBoltzmann/T] + 1);
(*Resampling the points*)
RandomSample[weightsRatio -> ptsBoltzmann, npoints]
], CompilationTarget -> "C", RuntimeOptions -> "Speed"]
FastFDsampler2[1, 3]; // AbsoluteTiming
sss2 = FastFDsampler2[10000, 3]; // AbsoluteTiming
Show[Histogram[{sss2, -3*Total[Log[#]] & /@
RandomVariate[UniformDistribution[{0, 1}], {10^4, 2 + 1}]}, 100,
"ProbabilityDensity"], Plot[DistrFD[3, x], {x, 0, 30}]]
Show[Histogram[{sss, sss2}, 100, "ProbabilityDensity"],
Plot[DistrFD[3, x], {x, 0, 30}]]


{0.0000323,Null}

{0.0107477,Null}

It first samples points according to the Boltzmann distribution, then computes the weights as the ratio of the values of FD and Boltzmann distributions, and then resamples. I assume that it would require a much smaller number of points - simply because the Boltzmann and FD distributions differ only in the domain of small x.

Is there any faster/more reliable approach for the specified conditions?

• Is there a known range of values of $T$?
– JimB
Commented Apr 13 at 17:09
• @JimB you can assume $T=1$ and then scale the results at the end. No need to worry about the value of $T$! Commented Apr 13 at 18:09
• @Roman. Good! That means HenrickShumacher 's inverse sample suggestion will work fine for the previous question at mathematica.stackexchange.com/questions/301633/….
– JimB
Commented Apr 13 at 18:53
• @JimB : would it be really an option given that I woild need to use interpolations? The built-in interpolation is very slow, and a custom one would still require some way of handling. Commented Apr 23 at 15:44
• I'm not understanding your question. First, I know nothing about using Compile. Is having a "piecewise-linear approximant" an impossibility when using Compile?
– JimB
Commented Apr 23 at 16:04

Perhaps its a good idea to use the InverseFunction of the CDF.

Make the functions parameter dependent an real numeric

    AvalFD[T_,x_]=Integrate[x^2/(Exp[x/T]+1),{x,0,Infinity},
Assumptions->T>0]//N;
DistrFD[T_,x_]:=1/AvalFD[T,x]*x^2/(Exp[x/T]+1);


Integrate the distribution, as for as it goes

     ICDFa[p_] :=
Evaluate[
InverseFunction[(
Function @@ {x,
Integrate[DistrFD[1, xi], {xi, 0, x},
Assumptions :> x > 0]})]
[p]]

ICDFa[y] // FullForm


$$\text{InverseFunction}\left[x \mapsto 1.10921 x \text{Li}_2\left(-e^{-x}\right)+1.10921 \text{Li}_3\left(-e^{-x}\right)-0.554605 x^2 \log \left(e^{-x}+1\right)+1.\right][y]$$

Making the function Listable makes it possible to use samples from uniform distribution in $$(0,1)$$ as argument

  SetAttributes[ICDFa, {NumericFunction, Listable}]


  Timing[sample = ICDFa[RandomReal[{0, 1}, 1000]];]
{1.25, Null}

ListPlot[{sample, Sort[sample]}]