Generate random points inside a Epanechnikov distribution

I am trying to generate points inside an Epanechnikov distribution given its parameters. If anyone has a sample code or can help me with this, I will be glad.

Answer 1

[Roger Stafford]

(See my second answer for a much faster method)

The following code generates random x values in accordance with the Epanechnikov distribution with parameters mu and r, as I understand it. How this relates to the 2D diagram you show, I cannot determine, but perhaps you can make use of this method in achieving what you have in mind.

mu = 7; r = 5; % <-- Choose parameters
n = 8192; % <-- Select desired number of points.
x = rand(1,n);
for k = 1:n
  t = roots([-1/4,0,3/4,2/4-x(k)]); % Get three roots of cubic
  x(k) = t(t>=-1&t<=1)*r + mu; % Discard two of them
end

% The following plot shows that x undoubtedly has 
% the correct Epanechnikov cumulative distribution:
xs = sort(x);
p = linspace(0,1,n);
x1 = linspace(mu-r,mu+r,n);
p1 = (x1-mu)/r;
p1 = (-p1.^3+3*p1+2)/4; % Epanechnikov cubic
plot(x1,p1,'bo',xs,p,'y-')

This code uses the 'roots' function to solve the above cubic polynomial, and is therefore not very efficient. I am sure there is a faster method of solving this cubic, but I am not aware of it at the moment.

[Walter Roberson]

It is possible to create exact formulas for solutions to cubics, but because the formulas are fairly long, using those formulas would likely be slower than using roots()

[Roger Stafford]

@Walter: I found a very much faster method on the internet which uses 'sin' and 'asin'. I have entered it as a second "answer". I have checked that it actually does solve that cubic.

Answer 2

[Roger Stafford]

Courtesy of the internet, I found a much, much faster method of solving the cubic below using the 'sin' and 'asin' functions.

The following code generates random x values in accordance with the Epanechnikov distribution with parameters mu and r, as I understand it. How this relates to the 2D diagram you show, I cannot determine, but perhaps you can make use of this method in achieving what you have in mind.

   mu = 7; r = 5; % <-- Choose parameters
   n = 8192; % <-- Select desired number of points
   x = 2*sin(asin(2*rand(1,n)-1)/3)*r+mu; % <-- Much faster

   % The following plot shows that x undoubtedly has 
   % the correct Epanechnikov cumulative distribution:
   x = sort(x);
   p = linspace(0,1,n);
   x1 = linspace(mu-r,mu+r,n);
   p1 = (x1-mu)/r;
   p1 = (-p1.^3+3*p1+2)/4; % Epanechnikov cubic
   plot(x1,p1,'bo',x,p,'y-')