Calculating the inverse of a CDF which is a polynomial

446 Views Asked by At

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$

1

There are 1 best solutions below

1
On BEST ANSWER

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.