I have a list of npts elements energies obeying some arbitrary distribution. By performing some operation many times, it would be converted into the sample obeying a distribution with the shape e^2*Exp[-e/T], where I call the elements by e and where T is determined from the value of npts and the total Total[energies]. To avoid performing this operation, I simply resample the npts points according to the final distribution (the corresponding array is called resampled). However, as this procedure violates the conservation law (Total[resampled] is not equal to Total[energies]) due to the probabilistic nature of sampling, I then manually rescaled all the values by the ratio Total[enegies]/Total[resampled].

Below, there is an example of my approach:

(*Preliminary definitions*)
norm[power_, Temp_] = 
  Integrate[p^power Exp[-p/Temp], {p, 0, Infinity}, 
   Assumptions -> power > 0 && Temp > 0];
DistrNormalized[p_, power_, Temp_] = (p^power Exp[-p/Temp])/
  norm[power, Temp];
Emean[power_, Temp_] = 
  Integrate[DistrNormalized[p, power, Temp]*p, {p, 0, Infinity}, 
   Assumptions -> power > 0 && Temp > 0];
Teff[power_, nparticles_, etot_] = 
 temp /. Solve[nparticles*Emean[power, temp] == etot, temp][[1]]
DistrSampler[nparticles_, T_, power_] := Module[{pts, weights, sample},
  pts = RandomReal[{0, 15 T}, Max[10*nparticles, 10^5]];
  weights = #^power Exp[-#/T] & /@ pts;
  sample = RandomChoice[weights -> pts, nparticles];
(*The actual example*)
(*Number of points*)
npts = 5*10^3;
(*Random energies*)
energies = RandomReal[{0, 5}, npts];
(*The total value of the energy*)
tot = Total[energies];
(*The effective temperature of the equilibrium shape*)
Tval = Teff[2, npts, tot];
(*Sampling the same number of points according to the distribution e^2Exp[-e/T] *)
resample = DistrSampler[npts, Tval, 2];
(*Rescaling their values to account for the conservation law*)
resample = resample*tot/Total[resample];

However, I feel that this naive way might be wrong: the operation changes not just the mean but also the variance of the sample in a way that does not correspond to the original distribution e^2 Exp[-e/T]. In my real simulation, which is related to the physics of binary interacting particles, it leads to the violation of the kinetic equilibrium.

What can be a more adequate but still fast way of producing resample?

  • $\begingroup$ I'm not a physicist so I don't understand the objective (or even if the objective makes any sense) but one can simplify things by noting that Teff can be more directly programmed with Teff[power_, nparticles_, etot_] := 2^(1/(-1 + power)) (etot/(nparticles Gamma[2 + power]))^(1/(-1 + power)). $\endgroup$
    – JimB
    Mar 18 at 1:14
  • $\begingroup$ Related: stats.stackexchange.com/questions/243260/…. $\endgroup$
    – JimB
    Mar 18 at 3:49

1 Answer 1


This is an extended comment. Is the question about the following: You want to have $n$ random samples from a distribution with density

$$f(e)=\frac{e^p \exp \left(-\frac{e}{T}\right)}{T^{p+1} \Gamma (p+1)}$$

where $p$ and $T$ are parameters and $0 \leq e < \infty$ conditioned on the sum of the samples is $s$ ?

$$\sum _{i=1}^n e(i)=s$$

It's not clear how $s$ is to be set. Is that an arbitrary value or the value of a particular random sample?

If you just resample the points from a particular sample and scale them to add up to the original sum, one can certainly do that (as one can do a lot of things to numbers and the numbers don't complain) but I'm not seeing what objective that achieves.


Not the answer you're looking for? Browse other questions tagged or ask your own question.