What you need to do is numerically evaluate the integral. There are a few ways of doing this in matlab. There is the "quad" family of functions. I have had problems with them when you want infinite bnds and it won't necessarily create a pdf that is monotonically increasing (since it approximates the function, and then integrates).
So probably the best way (?) is just numerically sample the function over reasonable bnds and a small spacing. Then just do a trapazoidal numeric integration.
So first figure out some reasonable bnds. From your formulation, you seem to know
a priori what the standard deviation is
.
So let
By Chebchev's inequality, you are guaranteed to miss at most 1% of the total area under the curve by using this as a bound. For most distributions, it is substantially better than that. You may want to crank that sucker down to 4 or 5, say, rather than 10.
Code:
x = linspace(-B,B,10000);
pdf = a*exp(-(b*x).^c);
% perform trapazoidal cumulative integration
cdf = cumtrapz(x,pdf);
You will be able to tell how well you did by looking at 1-cdf(end). If that is very small, then chances are you have a good sampling of the pdf. If you don't want your cdf to be quite that big you still need to calculate the cdf over a big range and small spacing (as I have done) and
then you can downsample.
For example:
Code:
x_small = -B:0.05:B;
cdf_small = interp1(x,cdf,x_small,'linear','extrap');