1
$\begingroup$

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.

$\endgroup$
6
  • 4
    $\begingroup$ 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) $\endgroup$
    – Bob Hanlon
    Apr 9 at 0:29
  • 3
    $\begingroup$ 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). $\endgroup$ Apr 9 at 1:46
  • $\begingroup$ @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? $\endgroup$ Apr 13 at 12:47
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 13 at 14:07
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 13 at 14:14

1 Answer 1

3
$\begingroup$
$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] &

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

For a -> 0 this is a GammaDistribution

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

(* True *)

Add limit definition

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

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

enter image description here

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 *)
$\endgroup$
1
  • $\begingroup$ Thanks! The only problem is that for small numbers of sampled points, the timing remains large. Say, 1 point still requires 0.01 s. $\endgroup$ Apr 13 at 12:54

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.