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 $\displaystyle \sigma_x$.

So let

$\displaystyle B = 10\, \sigma_x $

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');