Create a function that will return the cdf at any given point to implement the bisection code.

Trying to implement a random sampling algorithm via the transformation S = Fs-1(U) where U~U(0,1).

My given pdf is c*e^-x^3 for x>0, o for x<0. Solved c to be 1.1198. I used numerical integration to find the cdf, ~0.999958. But now I do not know how to invert the cdf to create random sampling. I tried to set up the bisection method below, to find a value x where h(x) = u, maintain values xmin and xmax where h(xmin) < u and h(xmax) > u. A point x is chosen at the midpoint between xmin and xmax. Based on the value of h(x), either xmin or xmax is updated to halve the difference between them, but I do not know how to properly do it.

I need to create a function that will return the cdf at any given point and ultimately plot an estimated pdf from these samples.

syms y
expr = 1.1198*exp(-y.^3); 
F = int(expr,0,inf);

h= F;
u = 1;
xmin = 1.5;
xmax = 3;
tol = 0.000000000001;

while xmax - xmin > tol
    x  = (xmax+xmin)/2 ;
    
    if h > u
        xmax = x  ;
    else
        xmin = x ;
    end
    XX=x
end

end
fprintf('It is %g\n', x)

[Image Analyst]

I formatted your post as code so people with the Symbolic Toolbox can copy it easily.

Did you try cumsum or cumsum?

Answer 1

[Torsten]

syms c x real
f = c*exp(-x.^3);
cnum = solve(int(f,x,0,Inf)==1,c);
f = subs(f,c,cnum);
F = int(f,0,x);
F = matlabFunction(F);
y = rand(500,1);
Finv = arrayfun(@(y)fsolve(@(x)F(x)-y,1,optimset('Display','off')),y);
Finv = Finv(abs(imag(Finv))<1e-8);
hold on
x = 0:0.1:5;
fnum = matlabFunction(f);
plot(x,fnum(x))
histogram(Finv,'Normalization','pdf')
hold off