I have a particular probability distribution (that is a function of a parameter $\phi$), that I would like to sample from and then do a listplot for that function.

Here is the code that I came up with:

m = 1;
\[Omega] = 1;
\[HBar] = 1;
fockelements = 8;

quadratures = 
    1/Sum[1/k!, {k, 0, fockelements}] Abs[
       Sum[(\[Alpha] E^(I \[Phi]))^n/\[Sqrt](n!) 1/
          Sqrt[2^n n!] ((m \[Omega])/(\[Pi] \[HBar]))^(1/
            4) Exp[-((m \[Omega] z^2)/(2 \[HBar]))] HermiteH[n, 
          Sqrt[(m \[Omega])/\[HBar]] z], {n, 0, 
         fockelements}]]^2, {z, -\[Infinity], \[Infinity]}] /. {\
\[Alpha] -> 1};

points = 3;
\[Phi]list = Subdivide[2 \[Pi], points];
data = RandomVariate[quadratures /. \[Phi] -> #] & /@ \[Phi]list

It seems to work but it's very slow. Origionally the code had a few functions-in-functions, and I tried putting everything in a compact form in hopes that it would speed-up, but it's still too slow.

Ideally I'd like to sample something like 1000 points, while right now it takes a very long time to do even 100 samplings.

EDIT: Not sure it helps to mention this but here are some warnings I get when I run my code:

enter image description here

Maybe this will make it more obvious what the issue is.

EDIT2: Just mentioning that I added a normalization constant to the code above as suggested by @JimB. (The code is still slow after normalization).

  • $\begingroup$ Please provide all definitions to run this. $\endgroup$ Commented Jan 12, 2021 at 9:33
  • $\begingroup$ @HenrikSchumacher, thanks for letting me know, I forgot to say what the variables are. $\endgroup$ Commented Jan 12, 2021 at 20:49
  • $\begingroup$ To make the function you give into a pdf that integrates to 1 you'll need to multiply it by Exp[-1] or use the Method->"Normalize" option. The former would probably be better. $\endgroup$
    – JimB
    Commented Jan 12, 2021 at 21:32
  • $\begingroup$ Your distribution is so very, very close to a Normal distribution with mean $\sqrt{2}\cos(\phi)$ and variance 3/4. Would you give your rationale for needing such a subtly different distribution? (I don't doubt you have a good rationale. I'm just curious.) $\endgroup$
    – JimB
    Commented Jan 12, 2021 at 21:52
  • 1
    $\begingroup$ Sorry, I was wrong. The multiplier should be 8!/109601 rather than Exp[-1]. If $n$ in the sum goes from 0 to $n_0$, then the multiplier should be $1/\sum_{k=0}^{n_0}(1/k!)$. As the number of terms goes up then the multiplier approaches Exp[-1]. $\endgroup$
    – JimB
    Commented Jan 12, 2021 at 23:13

1 Answer 1


This is more of an extended comment. Here is a slight modification of your code:

parms = {m -> 1, ω -> 1, ℏ -> 1, α -> 1, n0 -> 8};

quadratures = ProbabilityDistribution[(1/Sum[1/k!, {k, 0, n0}]) *
  Abs[Sum[(α E^(I ϕ))^n/√(n!) 1/Sqrt[2^n n!] ((m ω)/(π ℏ))^(1/4) *
  Exp[-((m ω z^2)/(2 ℏ))] HermiteH[n, Sqrt[(m ω)/ℏ] z], {n, 0, n0}]]^2,
  {z, -∞, ∞}] /. parms;

points = 10;
ϕlist = Subdivide[2 π, points];
AbsoluteTiming[Quiet[data = RandomVariate[quadratures /. ϕ -> #, 1000] & /@ ϕlist;]]
(* {5.89363, Null} *)

So that takes about 6 seconds to generate 10,000 samples (1,000 for each of 10 values of $\phi$. Getting 100,000 samples (10,000 samples for each of the 10 values of $\phi$) takes only 7 seconds. Just doing one at a time as your code is doing will be much slower. (I have to imagine that getting just a single random sample from each value of $\phi$ can't be of much use.)

The CDF when $\phi=2$ is

CDF[quadratures /. parms /. ϕ -> 2, z]

CDF of quadratures distribution with phi = 2

There likely isn't a nice closed-form inverse function to obtain random samples so I think there's not much you can do to speed things up other than to take multiple samples for each value of $\phi$.

  • $\begingroup$ Thanks JimB. Now I need to figure out how to put this result into the correct form to be plottable. $\endgroup$ Commented Jan 13, 2021 at 0:30
  • $\begingroup$ So far I tried data = Riffle[ RandomVariate[quadratures /. [Phi] -> #, 1000], #, {1, 1000, 2}] & /@ [Phi]list; and then ListPlot[ Partition[Flatten[data], 2]] - but there's some bug where it only does some of the data correctly. $\endgroup$ Commented Jan 13, 2021 at 0:34
  • $\begingroup$ Not understanding what kind of plot you want. (Each row in my example represents a different value of $\phi$.) Getting histograms for each is found with Histogram[#, "FreedmanDiaconis", "PDF"] & /@ data. $\endgroup$
    – JimB
    Commented Jan 13, 2021 at 2:04
  • $\begingroup$ I did it with: data = Riffle[ RandomVariate[quadratures /. [Phi] -> #, 10], #, {2, -1, 2}] & /@ [Phi]list;]] ListPlot[ Partition[Flatten[Reverse[data, 2]], 2] $\endgroup$ Commented Jan 13, 2021 at 2:13
  • $\begingroup$ Every 1000 elements pulls from the same value of theta, so I just need to generate a list where each element in the 1000 has its own associated theta value (ie, {{0, Rand[F[0]]}, {0, Rand[F[0]]},...(one-thousand samples with theta =0),..(more samples for other theta values),...,{pi, Rand[F[pi]]},...,{pi, Rand[F[pi]]}} $\endgroup$ Commented Jan 13, 2021 at 2:15

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.