0
$\begingroup$

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];
  sample
]
(*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?

$\endgroup$
2
  • $\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

2
$\begingroup$

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.

$\endgroup$

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