# Speeding up code doing random sampling

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;

ProbabilityDistribution[
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
ListPlot[data]


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: 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).

• Please provide all definitions to run this. Jan 12 at 9:33
• @HenrikSchumacher, thanks for letting me know, I forgot to say what the variables are. Jan 12 at 20:49
• 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.
– JimB
Jan 12 at 21:32
• 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.)
– JimB
Jan 12 at 21:52
• 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].
– JimB
Jan 12 at 23:13

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] 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$$.

• Thanks JimB. Now I need to figure out how to put this result into the correct form to be plottable. Jan 13 at 0:30
• 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. Jan 13 at 0:34
• 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.
– JimB
Jan 13 at 2:04
• I did it with: data = Riffle[ RandomVariate[quadratures /. [Phi] -> #, 10], #, {2, -1, 2}] & /@ [Phi]list;]] ListPlot[ Partition[Flatten[Reverse[data, 2]], 2] Jan 13 at 2:13
• 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, Rand[F]},...(one-thousand samples with theta =0),..(more samples for other theta values),...,{pi, Rand[F[pi]]},...,{pi, Rand[F[pi]]}} Jan 13 at 2:15