In [15]:
%matplotlib inline

Curve Fitting

In [17]:
import numpy as np
from scipy.optimize import curve_fit

In [18]:
import matplotlib.pyplot as plt

In [19]:
def func(x, a, b):
    return a * x + b

In [35]:
x = np.linspace(0, 10, 100)
y = func(x, 1, 2)

In [36]:
yn = y + 2 * np.random.normal(size=len(x))

In [37]:
print curve_fit.__doc__

    Use non-linear least squares to fit a function, f, to data.

    Assumes ``ydata = f(xdata, *params) + eps``

    f : callable
        The model function, f(x, ...).  It must take the independent
        variable as the first argument and the parameters to fit as
        separate remaining arguments.
    xdata : An M-length sequence or an (k,M)-shaped array
        for functions with k predictors.
        The independent variable where the data is measured.
    ydata : M-length sequence
        The dependent data --- nominally f(xdata, ...)
    p0 : None, scalar, or N-length sequence
        Initial guess for the parameters.  If None, then the initial
        values will all be 1 (if the number of parameters for the function
        can be determined using introspection, otherwise a ValueError
        is raised).
    sigma : None or M-length sequence, optional
        If not None, these values are used as weights in the
        least-squares problem.
    absolute_sigma : bool, optional
        If False, `sigma` denotes relative weights of the data points.
        The returned covariance matrix `pcov` is based on *estimated*
        errors in the data, and is not affected by the overall
        magnitude of the values in `sigma`. Only the relative
        magnitudes of the `sigma` values matter.

        If True, `sigma` describes one standard deviation errors of
        the input data points. The estimated covariance in `pcov` is
        based on these values.

    popt : array
        Optimal values for the parameters so that the sum of the squared error
        of ``f(xdata, *popt) - ydata`` is minimized
    pcov : 2d array
        The estimated covariance of popt. The diagonals provide the variance
        of the parameter estimate. To compute one standard deviation errors
        on the parameters use ``perr = np.sqrt(np.diag(pcov))``.

        How the `sigma` parameter affects the estimated covariance
        depends on `absolute_sigma` argument, as described above.

    See Also

    The algorithm uses the Levenberg-Marquardt algorithm through `leastsq`.
    Additional keyword arguments are passed directly to that algorithm.

    >>> import numpy as np
    >>> from scipy.optimize import curve_fit
    >>> def func(x, a, b, c):
    ...     return a * np.exp(-b * x) + c

    >>> xdata = np.linspace(0, 4, 50)
    >>> y = func(xdata, 2.5, 1.3, 0.5)
    >>> ydata = y + 0.2 * np.random.normal(size=len(xdata))

    >>> popt, pcov = curve_fit(func, xdata, ydata)


In [38]:
popt, pcov = curve_fit(func, x, yn)

In [39]:

array([ 1.05473697,  1.86758014])

In [40]:

array([[ 0.00579929, -0.02899645],
       [-0.02899645,  0.19428595]])

In [50]:
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
ax.scatter(x, yn, label="Noise Data")
ax.plot(x, func(x,popt[0], popt[1]), color="red", label="Fitted Line")
ax.plot(x, y, color="green", label="True Function")

<matplotlib.legend.Legend at 0x10b4c7550>

Solve Function

In [ ]:

In [51]:
from scipy.optimize import fsolve

In [52]:
print fsolve.__doc__

    Find the roots of a function.

    Return the roots of the (non-linear) equations defined by
    ``func(x) = 0`` given a starting estimate.

    func : callable ``f(x, *args)``
        A function that takes at least one (possibly vector) argument.
    x0 : ndarray
        The starting estimate for the roots of ``func(x) = 0``.
    args : tuple, optional
        Any extra arguments to `func`.
    fprime : callable(x), optional
        A function to compute the Jacobian of `func` with derivatives
        across the rows. By default, the Jacobian will be estimated.
    full_output : bool, optional
        If True, return optional outputs.
    col_deriv : bool, optional
        Specify whether the Jacobian function computes derivatives down
        the columns (faster, because there is no transpose operation).
    xtol : float
        The calculation will terminate if the relative error between two
        consecutive iterates is at most `xtol`.
    maxfev : int, optional
        The maximum number of calls to the function. If zero, then
        ``100*(N+1)`` is the maximum where N is the number of elements
        in `x0`.
    band : tuple, optional
        If set to a two-sequence containing the number of sub- and
        super-diagonals within the band of the Jacobi matrix, the
        Jacobi matrix is considered banded (only for ``fprime=None``).
    epsfcn : float, optional
        A suitable step length for the forward-difference
        approximation of the Jacobian (for ``fprime=None``). If
        `epsfcn` is less than the machine precision, it is assumed
        that the relative errors in the functions are of the order of
        the machine precision.
    factor : float, optional
        A parameter determining the initial step bound
        (``factor * || diag * x||``).  Should be in the interval
        ``(0.1, 100)``.
    diag : sequence, optional
        N positive entries that serve as a scale factors for the

    x : ndarray
        The solution (or the result of the last iteration for
        an unsuccessful call).
    infodict : dict
        A dictionary of optional outputs with the keys:

            number of function calls
            number of Jacobian calls
            function evaluated at the output
            the orthogonal matrix, q, produced by the QR
            factorization of the final approximate Jacobian
            matrix, stored column wise
            upper triangular matrix produced by QR factorization
            of the same matrix
            the vector ``(transpose(q) * fvec)``

    ier : int
        An integer flag.  Set to 1 if a solution was found, otherwise refer
        to `mesg` for more information.
    mesg : str
        If no solution is found, `mesg` details the cause of failure.

    See Also
    root : Interface to root finding algorithms for multivariate
    functions. See the 'hybr' `method` in particular.

    ``fsolve`` is a wrapper around MINPACK's hybrd and hybrj algorithms.


In [54]:
line = lambda x: x + 3

solution = fsolve(line, -2)
print solution


In [55]:
def findIntersection(func1, func2, x0):
    return fsolve(lambda x: func1(x) - func2(x), x0)

In [56]:
funky = lambda x: np.cos(x / 5.0) * np.sin(x/2.0)
line = lambda x: 0.01 * x - 0.5

In [60]:
x = np.linspace(0, 45, 10000)
roots = findIntersection(funky, line, [15, 20, 30, 35, 40, 45])

In [79]:
fig = plt.figure(figsize=(10, 6))

ax = fig.add_subplot(111)

line1 = ax.plot(x, funky(x), label="Funky function")
line2 = ax.plot(x, line(x), label="Line")
scat = ax.scatter(roots, line(roots))

ax.legend(loc = 0)

<matplotlib.legend.Legend at 0x10f052610>

In [ ]: