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];

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}]]



enter image description here

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[#]] & /@ 
      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}]]



enter image description here

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?

  • $\begingroup$ Is there a known range of values of $T$? $\endgroup$
    – JimB
    Commented Apr 13 at 17:09
  • $\begingroup$ @JimB you can assume $T=1$ and then scale the results at the end. No need to worry about the value of $T$! $\endgroup$
    – Roman
    Commented Apr 13 at 18:09
  • 1
    $\begingroup$ @Roman. Good! That means HenrickShumacher 's inverse sample suggestion will work fine for the previous question at mathematica.stackexchange.com/questions/301633/…. $\endgroup$
    – JimB
    Commented Apr 13 at 18:53
  • $\begingroup$ @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. $\endgroup$ Commented Apr 23 at 15:44
  • $\begingroup$ I'm not understanding your question. First, I know nothing about using Compile. Is having a "piecewise-linear approximant" an impossibility when using Compile? $\endgroup$
    – JimB
    Commented Apr 23 at 16:04

1 Answer 1


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

Make the functions parameter dependent an real numeric


Integrate the distribution, as for as it goes

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

      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}]

inverse of Fermi-Dirac CDF

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

 ListPlot[{sample, Sort[sample]}]

enter image description here


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.