Fitting a Power Law in IDL

This Help Article tells you how to fit a power law or an exponential to a set of points. The power law has the form y = a x^b, and the exponential models y = a exp(b x). The power law or exponential increases faster than a linear function, and a simple least-squares method will fail to converge. In this case, you can use logarithms to transform the equation into a linear problem, and find the solution using IDL's curve-fitting routines.
The purpose is to fit a set of data Y[i] to a power law, y = a x^b, where x is the independent variable and y is the dependent variable. The goal is to minimize the error between Y and y. The basic idea is to use the natural logarithm to transform y = a x^b into

ln(y) = ln(a) + b ln(x)

or:

y' = a' + b x'

where

y' = ln(y)
x' = ln(x)
a' = ln(a)

Thus, the power law is now just a straight line equation, and can be solved using any curve-fitting routine.

The same method can be used when fitting an exponential law. The equation y = a exp(b x) can be transformed into

ln(y) = ln(a) + b x

or:

y' = a' + b x

where

y' = ln(y)
a' = ln(a)

Examples of both of these models are given below.

; Examples of how to use CURVEFIT and COMFIT to fit
; a power law:  y = a x^b
 
;++++++++BEGIN:  Example Routines++++++++
 
;++++++++BEGIN:  Example User-Defined Procedure++++++++
; This routine is used with CURVEFIT.
PRO curvefit_power_law, x, coefs, y, pder
 
; Computing function to fit
y = ALOG(coefs[0]) + coefs[1]*x
 
; Computing derivatives of function to fit with respect
; to its coefficients.
IF (N_PARAMS() GE 4) THEN $
    pder = [[1./coefs[0] + x*0.],[x]]
END
 
;+++++++++END:  Example User-Defined Procedure+++++++++
 
;++++++++BEGIN:  Example Fitting Procedure++++++++
PRO fit_power_law
 
; Creating data.
x = DINDGEN(50) + 1
a = 2d
b = 1.8d
y = a*(x^b)
 
; Displaying original data.
PLOT, x, y
PRINT
PRINT,'                     a               b'
PRINT,'Actual:   ', a, b
 
; Calling CURVEFIT and outputing results.
; *** use the natural log model: 
; ln(y) = ln(a) + b ln(x)
coefs = [1d, 1.D] ; Initial guesses.
weights = y*0 + 1.
y_fit = CURVEFIT(ALOG(x), ALOG(y), weights, coefs, $
   FUNCTION_NAME = 'curvefit_power_law')
 
PRINT, 'CURVEFIT: ', coefs
OPLOT, x, EXP(y_fit), PSYM=2
PLOTS, 5,2000,PSYM=2
XYOUTS, 6,2000,'CURVEFIT'
 
; Calling COMFIT and outputing results.
; *** use the LOGSQUARE model:   
; y' = a' + b log(x) + c log(x)^2
;                              where y' = log(y)
;                                    a' = log(a)
coefs = [ALOG10(2.D), 1.D, 0D] ; Initial guesses.
results = COMFIT(x, ALOG10(y), coefs, /LOGSQUARE, YFIT = y_fit)
results[0] = 10.^results[0]  ; convert first coefficient back
 
PRINT, 'COMFIT:   ', results
OPLOT, x, 10.^y_fit, PSYM=4
PLOTS, 5,1800,PSYM=4
XYOUTS, 6,1800,'COMFIT'
 
END
 
;+++++++++END:  Example Fitting Procedure+++++++++
 
;+++++++++END:  Example Routines+++++++++