2
$\begingroup$

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 = 
  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:

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

$\endgroup$
12
  • $\begingroup$ Please provide all definitions to run this. $\endgroup$ – Henrik Schumacher Jan 12 at 9:33
  • $\begingroup$ @HenrikSchumacher, thanks for letting me know, I forgot to say what the variables are. $\endgroup$ – Steven Sagona Jan 12 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 Jan 12 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 Jan 12 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 Jan 12 at 23:13
4
$\begingroup$

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

$\endgroup$
5
  • $\begingroup$ Thanks JimB. Now I need to figure out how to put this result into the correct form to be plottable. $\endgroup$ – Steven Sagona Jan 13 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$ – Steven Sagona Jan 13 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 Jan 13 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$ – Steven Sagona Jan 13 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$ – Steven Sagona Jan 13 at 2:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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