Logarithmic and Exponential Curve Fit in Python - Numpy

With numpy function "polyfit":

X,y : data to be fitted

import numpy as np

1. Exponential fit

cf = np.polyfit(X, np.log(y), 1)

will return two coefficients, who will compose the equation:

exp(cf[1])*exp(cf[0]*X)


2. Logarithm fit:

cf = np.polyfit(np.log(X), y, 1)

will return two coefficients, who will compose the equation:

cf[0]*log(X)+cf[1]

Polynomial Curve Fitting

The code below shows how easily you can do a Polynomial Curve Fitting with Python and Numpy.

import numpy as np

# sample x and y data - example
x = [7.76,10.11,11.89,14.81,15.49]
y = [1.851,1.971,1.953,1.842,1.805]

# the polyfit functions does the nth degree polynomial best fit on the data, 
# returning the polynomial coefficients

n = 4   # 4th degree polynomial, you can change for whatever degree you want
coefs = np.polyfit(x,y,n)

# The poly1d function applies the polynomial function to our calculated coefficients
polyf = np.poly1d(coefs)

#if we want to apply our polynomial function to a range of x values
xf = np.linspace(0,20)
yf = polyf(xf)

Alternating Block Hyetograph Method with Python

To generate hypothetic storm events, we can use some methods as Alternating block method, Chicago method, Balanced method, SCS Storms among others.

In this post I show a way to generate a hypothetic storm event using the Alternating block method.

This python code uses the numpy library. The altblocks functions uses as input the idf parameters as list, and total duration, delta time and return period as floats.

import numpy as np

def altblocks(idf,dur,dt,RP):
    aDur = np.arange(dt,dur+dt,dt)    # in minutes
    aInt = (idf[0]*RP**idf[1])/((aDur+idf[2])**idf[3])  # idf equation - in mm/h
    aDeltaPmm = np.diff(np.append(0,np.multiply(aInt,aDur/60.0)))
    aOrd=np.append(np.arange(1,len(aDur)+1,2)[::-1],np.arange(2,len(aDur)+1,2))
    prec = np.asarray([aDeltaPmm[x-1] for x in aOrd])
    aAltBl = np.vstack((aDur,prec))
    return aAltBl

Make numpy array of 'datetime' between two dates

A simple way to create an array of dates (time series), between two dates:

We can use the numpy arange - https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html , function which is most used to create arrays using start / stop / step arguments.

Syntax:
numpy.arange([start, ]stop, [step, ]dtype=None)

In case of datetime values, we need to specify the step value, and the correct type and unit of the timestep in the dtype argument

. dtype='datetime64[m]' will set the timestep unit to minutes;
. dtype='datetime64[h]' will set the timestep unit to hours;
. dtype='datetime64[D]' will set the timestep unit to days;
. dtype='datetime64[M]' will set the timestep unit to months;
. dtype='datetime64[Y]' will set the timestep unit to months;

For example:

import numpy as np

dates = np.arange('2017-06-01', '2017-06-02', 15, dtype='datetime64[m]')
# 15 is the timestep value, dtype='datetime64[m] means that the step is datetime minutes

This example will create an array of 96 values, between 01jun2017 and 02jun2017, with a time step of 15 minutes.

Exponential curve fit in numpy

With numpy function "polyfit" we can easily fit diferent kind of curves, not only polynomial curves.

According to the users manual, the numpy.polyfit does:

"
Least squares polynomial fit.

Fit a polynomial p(x) = p[0] * x**deg + ... + p[deg] of degree deg to points (x, y). Returns a vector of coefficients p that minimises the squared error.
"

If we use X and y as arrays with our data, the code:

coef = np.polyfit(X, np.log(y), 1)

will return two coefficients, who will compose the equation:

exp(coef[1])*exp(coef[0]*X)

Giving you the exponential curve that better fits our data - X and y.
The polyfit function can receive weight values, which we can use in case of giving less importance to very small values, for example. We can use a weight function as following:

coef = np.polyfit(X, np.log(y), 1, w=np.sqrt(y))

Giving more weight to higher values.

To retrieve the R-squared index of our exponenctial curve, we can use de scikit r2_score, as following:
y_pred = np.exp(coefs[1])*np.exp(coefs[0]*X)

from sklearn.metrics import r2_score

r2s = r2_score(y, y_pred, sample_weight=None, multioutput=None)

Numpy - Accumulated and Incremental series

In Hydrology, it is always needed to deal with time-series of variables, as flow series or precipitation series, with the variable being incremental or accumulated.

Numpy has a great way to transform between accumulated and incremental series.

To accumulate a incremental series use the method

   numpy.cumsum(incrementalSeries)

And to transform a accumulated array to a incremental one, use:

    numpy.diff(accumulatedSeries)