# How to generate random value under the given distribution fast?

Consider a distribution

distr[EN_, Enu_, mN_] =
Exp[-(Sqrt[EN^2 + mN^2]/Enu)]*Sqrt[(Sqrt[Enu^2-mN^2]/2)^2 - EN^2] Sqrt[
EN^4 - (Enu/4)^4]


where mN is a free positive parameter, Enu>Sqrt[5]mN, and EN ranges within

ENmin[Enu_, mN_] = Max[mN, Enu/4];
ENmax[Enu_, mN_] = Sqrt[Enu^2-mN^2]/2;


distr is a probability distribution before the normalization, where Enu,mN are parameters, and EN is a variable.

Next, consider a dataset with Enu:

EnuVals = RandomReal[{0,60},10^6];


and a grid for mN:

mNvals=Table[x,{x,0.1,1.,0.1}];


My goal is to generate random values of EN (following the distribution distr) for the given values of Enu and mN. I.e. the desired final data would have rows

tabfinal={Enuvals[[i]],mNvals[[j]],ENrand[Enuvals[[i]],mNvals[[j]]]}


I do not have an idea how to quickly compute ENrand given the distribution distr. For particular values of Enu,mN, say Enuval,mNval, I would do the following:

ENvalsRandom = RandomReal[{ENmin[Enuval, mNval],ENmax[Enuval, mNval]}, 10000];
weights1 =
Hold@Compile[{{ENvalsRandom, _Real, 1}},
Table[distr[ENvalsRandom[[i]],Enuval,mNval], {i, 1, 10000, 1}],
CompilationTarget -> "C", RuntimeOptions -> "Speed"] /.
DownValues@distr//
ReleaseHold;
weights = weights1[ENvalsRandom];
weightedData = WeightedData[tabRPPS, weights];
edist = EmpiricalDistribution[weightedData];
ENvalRandomTrue = RandomVariate[edist, 1]


But this way is really slow since it requires building the empirical distribution Length[mNvals]*Length[EnuVals] times.

Could you please suggest how to optimize it?

• I can't understand what random variable is of interest. And what is distr? Is that a probability density function? cumulative probability density function? Something else? Are all of the arguments for distr parameters? Is one of those the random variable of interest?
– JimB
Sep 6, 2022 at 20:57
• @JimB : Thanks! I have updated the question. distr is a probability distribution before the normalization, where Enu,mN are parameters, and EN is a variable. Sep 6, 2022 at 21:04
• The random selection of Enu values will need to be changed from RandomReal[{0, 60},10^6] to RandomReal[{Sqrt[5] mN, 60},10^6] because if $0<mN\leq 1$, then $Enu>\sqrt{5}mN$.
– JimB
Sep 9, 2022 at 2:47

I suspect that you'll have to give up on "fast" for your un-normalized density function, as you'll need to normalize it numerically with specific values for mN and Enu and use the "invert the cdf" approach where both the "invert" and "cdf" parts both need to be found numerically. If there was a closed-form for the cdf (which seems doubtful), that would help.

Here is a brute-force approach when specifying values for mN and Enu:

(* Set parameters *)
mN = 1;
Enu = 5;

(* Constant of integration *)
const = NIntegrate[Exp[-(Sqrt[en^2 + mN^2]/Enu)]*
Sqrt[(Sqrt[Enu^2 - mN^2]/2)^2 - en^2] Sqrt[en^4 - (Enu/4)^4],
{en, Max[mN, Enu/4], Sqrt[Enu^2 - mN^2]/2}]
(* 3.03052 *)

(* cdf of random variable EN *)
cdf[EN_] := Module[{},
If[EN < Max[mN, Enu/4], 0,
If[EN > Sqrt[Enu^2 - mN^2]/2, 1,
NIntegrate[Exp[-(Sqrt[en^2 + mN^2]/Enu)]*
Sqrt[(Sqrt[Enu^2 - mN^2]/2)^2 - en^2] Sqrt[en^4 - (Enu/4)^4]/const,
{en, Max[mN, Enu/4], EN}] ]]]

(* The cdf has a uniform(0,1) so we take a random sample from a uniform distribution
and then solve for the corresponding EN *)
(* Take n random samples *)
SeedRandom[12345];
n = 1000;
(* Starting value for EN *)
en0 = (Max[mN, Enu/4] + Sqrt[Enu^2 - mN^2]/2)/2;
AbsoluteTiming[
sample = Table[EN /. FindRoot[cdf[EN] == RandomReal[{0, 1}], {{EN, en0}}], {i, n}];
]
(* {20.7947, Null} *)

(* Show results *)
Show [Histogram[sample, "FreedmanDiaconis", "PDF", ImageSize -> Medium,
PlotRange -> {{1, 2.5}, Automatic}],
Plot[Exp[-(Sqrt[en^2 + mN^2]/Enu)]*
Sqrt[(Sqrt[Enu^2 - mN^2]/2)^2 - en^2] Sqrt[en^4 - (Enu/4)^4]/const,
{en, Max[mN, Enu/4], Sqrt[Enu^2 - mN^2]/2}]]