1

I would like to generate random numbers from a CDF. I have looked at the inverse sampling method but I am struggling to find the inverse of my CDF given that it is a 6th order polynomial.

$F_X \left ({x} \right)= Ax^2-Bx^4 + Cx^6$

Where $A$, $B$ and $C$ are real positive and non-zero. Is there a different method to generate a random number based on this CDF? If so, I would really appreciate any guidance (especially for closed-form solutions).

My "tries": I have learnt about the Box-Mueller but it cannot be applied for this case (correct me if I am wrong). Secondly, I have had a look at rejection sampling but I don't quite understand how to apply it (working on MATLab by the way).

Edit: The region of interest is $0 < x < 80$

Aviram Serra
  • 95
  • You might consider working with the density, its derivative, instead and using an acceptance-rejection method. Otherwise simply use the CDF and the inverse method, with Newton’s method to solve for what you need. – A rural reader Mar 13 '21 at 22:59
  • @Aruralreader thanks bud! I was having a look at the inverse sampling method, I have found the roots and am now trying to figure out the inverse. I have also found an interesting method which is the Lagrange inversion (to solve a polynomial) but it doesnt seem useful in my scenario. – Aviram Serra Mar 14 '21 at 01:15
  • Unless you have $x$ bounded, you don't have a CDF as $F_X(x)$ will be larger than 1 when $x$ gets large in absolute value. Do you know what the bounds are or are those also to be estimated from data? – JimB Mar 14 '21 at 21:10
  • @JimB thanks for the response. The region is $0 < x < 80$. You can see the specific plot here link – Aviram Serra Mar 14 '21 at 21:13
  • You need to add the information about the region for $x$ to your question. – JimB Mar 14 '21 at 21:54
  • You can use WolframAlpha as you did with the plot of the cdf to find the inverse function: solve (3.4675*10^-4)*x^2 -(3.5116*10^-8)*x^4+(8.361*10^-13)*x^6 ==p for x. – JimB Mar 15 '21 at 00:49
  • 1
    If you want a quick way to knock this out, start by figuring out the maximum of $f(x) = 2Ax - 4Bx^3 + 6Cx^5$; call it $f_{\text{max}}$. Then generate $X$ uniformly between $0$ and $80$ and $U$ between $0$ and $f_{\text{max}}$. If you get $U \leq f(X)$ accept $X$, otherwise repeat until you do. Not optimal for sure but you do avoid having to solve an equation for each $X$ you generate. Doesn't answer your question, which is about solving $F(x) = u$, but such as it is maybe it's helpful. – A rural reader Mar 15 '21 at 01:05

1 Answers1

1

You can use Wolfram Alpha as you did with the plot of the cdf to find the inverse function: solve (3.4675*10^-4)*x^2 -(3.5116*10^-8)*x^4+(8.361*10^-13)*x^6 ==p for x.

If you have Mathematica, you can generate samples in the following manner:

(* CDF and PDF *)
cdf = 3.4675 10^(-4) x^2 - 3.5116 10^(-8) x^4 + 8.361 10^(-13) x^6;
pdf = D[cdf, x];

(* Generate random samples *)
n = 1000;
data = x /. Solve[cdf == #][[4]] & /@ RandomReal[{0, 1}, n];


(* Fit a nonparametric density estimate and plot results *)
skd = SmoothKernelDistribution[data, 
   Automatic, {"Bounded", {0, 80}, "Gaussian"}];
Plot[{pdf, PDF[skd, x]}, {x, 0, 80}, 
  PlotLegends -> {"True density", "Estimated density"}]

True and estimated pdf

The appropriate solution for $x$ given $p$ (a random sample from a Uniform[0,1] distribution) in the following equation

$$3.4675*10^{-4} x^2 - 3.5116*10^{-8} x^4 + 8.361*10^{-13} x^6 = p$$

is given by

x /. Solve[cdf == p, x][[4]]

General result for obtaining a random sample

The imaginary parts cancel out (eventually). I don't know if this completely qualifies as a case of casus irreducibilis.

JimB
  • 1,861
  • Yes, it's a casus irreducibilis for any probability $\ p\ $, unless $\ p\in F_X(\mathbb{Q})\ $, which can be true for only a countable number of $\ p\ $. The polynomial $\ F_X(x)\ $ is a cubic in $\ x^2\ $ whose derivative has zeroes at $\ x^2\approx 6400.17\ $ (where $\ F_X(x)\approx1.00003\ $), and at $\ x^2\approx21599.34\ $ (where $\ F_X(x)\approx-0.468\ $). So for any $\ p\in(0,1]\ $ the equation $\ F_X(x)=p\ $ has six real roots, exactly one of which lies in the interval $\ (0,80]\ $. – lonza leggiera Mar 15 '21 at 11:29