# How to quickly sample from the following distributions

Consider the following distributions:

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


Here, T>0 is a parameter, and 0<x<Infinity.

How can one sample n random points according to the distribution?

weightscomp =
Compile[{{data, _Real},{sign,_Real}, {T, _Real}}, Module[{vals},
vals = data;
vals^2/(Exp[vals/T]+stat)
]
, CompilationTarget -> "C", RuntimeOptions -> "Speed",
RuntimeAttributes -> {Listable}]
DistrSampler[nparticles_, T_, stat_] := Module[{(*pts,weights,sample*)},
sign=If[stat=="FD",1.,-1.];
pts = RandomReal[{0, 15 T}, Max[10*nparticles, 10^5]];
weights = weightscomp[pts, T,sign];
sample = RandomSample[weights -> pts, nparticles];
sample
]


But it is relatively slow. In principle, to produce a small number of points adequately, I need to sample a huge number of random points and weights for them. I am especially interested in a code that may be compiled.

• Your normalization of the distributions do not appear to be correct. Presumably you want DistrBE[T_, x_] = 1/AvalBE*x^2/(Exp[x/T] - 1) and DistrFD[T_, x_] = 1/AvalFD*x^2/(Exp[x/T] + 1) Apr 9 at 0:29
• Use inverse sampling: en.wikipedia.org/wiki/Inverse_transform_sampling. The CDF seems to be easy enough to obtain symbolically. To invert it, you could form a piecewise-linear approximant and invert that (just reflect the graph about the angle bisector of the coordinate system). Apr 9 at 1:46
• @HenrikSchumacher : for my purposes, I will need to sample 1-2 points from the mentioned distributions, assuming a dynamically updating temperature inside Do many times. How do you think, would the method you proposed (using the piecewise-linear approximant) work efficiently? Apr 13 at 12:47
• That means you discard the precious distribution after two samples? Oh, then forget about the method with the piecewise-linear that I mentioned. The cost of building this approximation will ammortized if the number of samples is substantially greater than the number of pieces in the piecewise-linear approximation. Apr 13 at 14:07
• Maybe you are better off with some simple rejection sampling. PDF[GammaDistribution[2, 2 T], x] seems to be a decent approximation of DistrBE[T,x] to make it work. Apr 13 at 14:14

\$Version

(* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *)

Clear["Global*"]


Define your distributions with ProbabilityDistribution so that they have behave like the built-in distributions.

DistBE[T_] = ProbabilityDistribution[x^2/(Exp[x/T] - 1), {x, 0, Infinity},
Assumptions -> T > 0, Method -> "Normalize"];

({#, #[DistBE[T]]} & /@ {PDF[#, x] &, CDF[#, x] &, Mean, Variance,
DistributionParameterAssumptions}) /.
{Function[PDF[Slot[1], x]] :> PDF,
Function[CDF[Slot[1], x]] :> CDF} //
Grid[#, Frame -> All] &


Plot[Evaluate@Table[PDF[DistBE[T], x], {T, 1, 2.5, 0.5}],
{x, 0, 15},
PlotRange -> All,
PlotLegends -> Placed[LineLegend[Range[1, 2.5, 0.5],
LegendLabel -> T], {.7, .7}]]


Generating 10^5 samples:

SeedRandom[1234];
AbsoluteTiming[(dataBE = RandomVariate[DistBE[1], 10^5])[[1 ;; 4]]]

(* {0.075116, {0.729519, 1.877, 0.246637, 5.51877}} *)


Similarly with DistFD

DistFD[T_] = ProbabilityDistribution[x^2/(Exp[x/T] + 1), {x, 0, Infinity},
Assumptions -> T > 0, Method -> "Normalize"];

({#, #[DistFD[T]]} & /@ {PDF[#, x] &, CDF[#, x] &, Mean, Variance,
DistributionParameterAssumptions}) /.
{Function[PDF[Slot[1], x]] :> PDF,
Function[CDF[Slot[1], x]] :> CDF} //
Grid[#, Frame -> All] &


Plot[Evaluate@Table[PDF[DistFD[T], x], {T, 1, 2.5, 0.5}],
{x, 0, 15},
PlotRange -> All,
PlotLegends -> Placed[LineLegend[Range[1, 2.5, 0.5],
LegendLabel -> T], {.7, .7}]]


Generating 10^5 samples:

SeedRandom[1234];
AbsoluteTiming[(dataFD = RandomVariate[DistFD[1], 10^5])[[1 ;; 4]]]

(* {0.075829, {1.16232, 2.36629, 0.539816, 5.92021}} *)


EDIT: The two distributions can be reduced to a single two-parameter distribution

dist[T_, a_] = ProbabilityDistribution[x^2/(Exp[x/T] + a), {x, 0, Infinity},
Assumptions -> T > 0 && a >= -1, Method -> "Normalize"];

({#, #[dist[T, a]]} & /@ {PDF[#, x] &, CDF[#, x] &, Mean, Variance,
DistributionParameterAssumptions}) /.
{Function[PDF[Slot[1], x]] :> PDF, Function[CDF[Slot[1], x]] :> CDF} //
Grid[#, Frame -> All] &


For a -> 0 this is a GammaDistribution

Limit[PDF[dist[T, a], x], a -> 0] == PDF[GammaDistribution[3, T], x]

(* True *)


dist[T_, 0] = GammaDistribution[3, T];

Plot[Evaluate@Table[Tooltip[PDF[dist[1.5, a], x], a], {a, -1, 2}], {x, 0, 10}]


SeedRandom[1234];
AbsoluteTiming[(dataBE2 = RandomVariate[dist[1, -1], 10^5])[[1 ;; 4]]]

(* {0.075883, {0.729519, 1.877, 0.246637, 5.51877}} *)

dataBE === dataBE2

(* True *)

SeedRandom[1234];
AbsoluteTiming[(dataFD2 = RandomVariate[dist[1, 1], 10^5])[[1 ;; 4]]]

(* {0.075216, {1.16232, 2.36629, 0.539816, 5.92021}} *)

dataFD === dataFD2

(* True *)
`
• Thanks! The only problem is that for small numbers of sampled points, the timing remains large. Say, 1 point still requires 0.01 s. Apr 13 at 12:54