Python: Plotting a power law function with exponential cutoff

Question

I have a graph between 2 functions f and g. I know it follows a power law function with exponential cutoff.

f(x) = x**(-alpha)*e**(-lambda*x)

How do I find the value of exponent alpha?

Answer 1

If you know that the points follow this law exactly, then invert the equation and put in an x and its corresponding f(x) value:

import math
alpha = -(lamda*x + math.log(f(x)))/math.log(x)

But the if the points do not exactly fit the equation you will need to do some sort of regression to determine alpha.

EDIT: Ok, so they don't fit exactly. This is getting beyond a Python question, but there may be something in numpy that can handle it. Here is a numpy linear regression recipe but your equation can't be rearranged into a linear form, so you'll have to look into non-linear regression.

Answer 2

If you have sufficiently close x points (for example one every 0.1), you can try the following:

ln(f(x)) = -alpha ln(x) - lambda x
ln(f(x))' = - alpha / x - lambda

So depending on where you have your points: If you have a lot of points near 0, you can try:

h(x) = x ln(f(x))' = -alpha - lambda x

So the limit of the function h when x goes to 0 is -alpha If you have large values of x, the function x -> ln(f(x))' tends toward lambda when x goes to infinity, so you can guess lambda and use pwdyson's expression.

If you don't have close x points, the numerical derivative will be very noisy, so I would try to guess lambda as the limit of -ln(f(x)/x for large x's...

If you don't have large values, but a large number of x's, you can try a minimization of

sum_x_i (ln(y_i) + alpha ln(x_i) + lambda x_i) ^2

on both alpha and lambda (I guess It would be more precise than the initial expression)... It is a simple least square regression (numpy.linalg.lstsq will do the job). So you have plenty of methods, the one to chose really depends on you inputs.

Answer 3

The usual and general way of doing what you want is to perform a non-linear regression (even though, as noted in another response, it is possible to linearize the problem). Python can do this quite easily with the help of the SciPy package, which is used by many scientists.

The routine you are looking for is its least-square optimization routine (scipy.optimize.leastsq). Once you wrap your head around the way this general optimization procedure works (see the example), you will probably find many other opportunities to use it. Basically, you calculate the list of differences between your measurements and their ideal value f(x), and you ask SciPy to find the parameters that make these differences as small as possible, so that your data fits the model as well as possible. This then gives you the parameter you are looking for.