"vinh" wrote in message <iqfq3s$gqa$1@newscl01ah.mathworks.com>...
> Hello,
>
> I've done my fair amount of searching and the answer seems within reach, however my attempts at solving this problem haven't worked out too well.
>
> I wish to use the rand function to generate random numbers w/ uniform distribution. My domain is from 1 to 1 so I know how to extend my result from rand to those boundaries. However, my distribution function is given as
>
> g = 1  x^2
>
> My thought process has been to integrate my distribution function to get my density function x  x^3/3 ... the inverse of this should allow me to use my result from rand to generate random numbers in accordance with my distribution function g. However, running a histogram on the result is completely incorrect...I've also tried using the inverse of my distribution function g but that doesn't work either. Could someone help me with the correct way to think about this?
           
I imagine you have already found that your g(x) needs to be multiplied by 3/4 to make it a valid probability density for x between 1 and +1, and you have undoubtedly found that the cumulative distribution function would then be
G(x) = 1/4*x^3+3/4*x+1/2
The hard part of your task is to then find the inverse to this function so that you can use 'rand'. In other words given y you are trying to find x such that y = 1/4*x^3+3/4*x+1/2 for any y in 0 <= y <= 1. Cubic equations like this do have explicit solutions but they are much messier than quadratic equations. The solution to this one for x in 1 <= x <= 1 and y in 0 <= y <= 1, courtesy of Cardan's formulas (1545), is as follows. Given a vector of y values, do this:
a = 12*y;
b = 2*sqrt(y.*(1y));
x = real((1/21/2*sqrt(3)*i)*(a+b*i).^(1/3)+...
(1/2+1/2*sqrt(3)*i)*(ab*i).^(1/3));
Then x will be the required corresponding roots. You can test this by computing 1/4*x.^3+3/4*x+1/2 to see if you come back to the original y values again (within round off error of course.)
Hence if you do
y = rand(1,n);
the above procedure will produce corresponding x values that possess the density distribution you desire, namely g(x) = 3/4*(1x^2).
Note: In theory the above values of x should already be real without applying the 'real' operator, but due to round off error they will usually have very small imaginary parts, so the 'real' operator is an easy way to remove these errors.
Roger Stafford
