In [15]:

    

%matplotlib inline


    

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``

    Parameters
    ----------
    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.

    Returns
    -------
    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
    --------
    leastsq

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

    Examples
    --------
    >>> 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]:

    

popt


    

    
Out[39]:





array([ 1.05473697,  1.86758014])

In [40]:

    

pcov


    

    
Out[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")
ax.legend(loc=2)


    

    
Out[50]:





<matplotlib.legend.Legend at 0x10b4c7550>






    

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.

    Parameters
    ----------
    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
        variables.

    Returns
    -------
    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:

        ``nfev``
            number of function calls
        ``njev``
            number of Jacobian calls
        ``fvec``
            function evaluated at the output
        ``fjac``
            the orthogonal matrix, q, produced by the QR
            factorization of the final approximate Jacobian
            matrix, stored column wise
        ``r``
            upper triangular matrix produced by QR factorization
            of the same matrix
        ``qtf``
            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.

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

    

In [54]:

    

line = lambda x: x + 3

solution = fsolve(line, -2)
print solution


    

    




[-3.]

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)


    

    
Out[79]:





<matplotlib.legend.Legend at 0x10f052610>






    

In [ ]: