0
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – JimB
    Sep 6, 2022 at 20:57
  • $\begingroup$ @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. $\endgroup$ Sep 6, 2022 at 21:04
  • 1
    $\begingroup$ 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$. $\endgroup$
    – JimB
    Sep 9, 2022 at 2:47

1 Answer 1

2
$\begingroup$

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

Histogram and true density

$\endgroup$

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.