"Python" III, scipy--"doing scientific calculations in Python" __python

1. Least Squares fittingAssuming that there is a set of experimental data (X[i], y[i]), we know the functional relationship between them: Y = f (x), and by using these known information, you need to determine some parameter items in the function. For example, if f is a linear function f (x) = K*x+b, then parameters K and B are the values we need to determine. If these arguments are expressed in p, then we are going to find a set of P values that make the S function in the following formula the smallest:

This algorithm is called least squares fitting (least-square fitting).
--------------------------------------------------------------------------------------------------------------- ----------------
Using a known dataset (X,y) to determine a set of parameters (K,b) makes the actual Y-value minus the predicted Y-value squared and smallest.
--------------------------------------------------------------------------------------------------------------- ----------------
The function leastsq of the least squares fitting method has been provided in the Optimize library of scipy. The following is an example of data fitting with LEASTSQ:

#-*-Coding:utf-8-*-
# Note If you have Chinese in your code, you must add the above sentence to the beginning of the file, or add "#encoding: Utf-8".
import NumPy as NP
from scipy.optimize import leastsq
import pylab as pl

def func (x, p):
    "" "Data Functions used: A*sin (2*pi*k*x + theta) "" "
    A, k, theta = P return
    a*np.sin (2*np.pi*k*x+theta)   

def residuals (p , y, x): "" "The
    difference between the experimental data x, Y and the fitted function, p is the coefficients to be found for fitting
    " "" Return
    y-func (x, p)

x = np.linspace (0, -2* Np.pi
A, k, theta = ten, 0.34, NP.PI/6 # Real data function parameter
y0 = func (x, [A, K, Theta]) # real data
y1 = y0 + 2 * NP . RANDOM.RANDN (len (x)) # After adding noise to the experimental data    

p0 = [7, 0.2, 0] # First guessing function fitting parameter

# calling LEASTSQ for data fitting
# residuals a function for calculating errors
# P0 for fitting parameter's initial value
# args for experimental data to be fitted
PLSQ = LEASTSQ (residuals, P0, args= (y1, x))

print U "Real parameters:", [A, K, th ETA]
print u "fitted parameters", Plsq[0] # Experimental data fitted with parameters

Pl.plot (x, y0, label=u "real data")
Pl.plot (x, y1, label=u "experimental data with Noise" )
Pl.plot (x, func (x, plsq[0]), label=u "fitted data")
Pl.legend () pl.show ()

The function we want to fit in this example is a sinusoidal wave function, which has three parameters a, K, and Theta, corresponding to amplitude, frequency and phase angle respectively. Suppose our experimental data is a set of noise-containing data x, y1, where Y1 is added to the noise based on the real data y0.

Through the LEASTSQ function of the experimental data with the noise of x, y1 data fitting, you can find the X and real data y0 the sine relationship between the three parameters: A, K, Theta.
Finally, we can find that the fitting parameters are quite different from the real parameters, but because of the periodicity of the sine function, the function of the fitting parameters is the same as that of the real parameters.

2. function Minimum valueThe Optimize library provides several algorithms for finding the minimum value of a function: Fmin, Fmin_powell, FMIN_CG, Fmin_bfgs. The following program demonstrates fmin functionality by solving the inverse of convolution.

For a discrete linear time invariant system H, if its input is x, then its output y can be represented by the convolution of X and H:

The question now is how to compute the system's transfer function H if you know the input x and output y of the system, or how to compute the input x of the system if you know the transfer function H of the system and the output y of the system. This kind of operation is called deconvolution operation, it is very difficult, especially in the practical application, the output of the measurement system always has the error.

Using Fmin to compute deconvolution, this method can only be used on a small number of series, so there is no great practical value, but to evaluate the performance of the Fmin function is good.

#-*-Coding:utf-8-*-# This procedure uses various Fmin functions to find the convolution inverse import scipy.optimize as opt import numpy as NP def Test_fmin_convolve (FMI Nfunc, X, h, y, yn, x0): "" "X (*) H = y, (*) indicates convolution yn the result of adding noise to the base of Y x0 to solve X's initial value" "" Def Convolve _func (h): "" "" "The Power Fmin of Yn-x (*) h will be calculated to make this power the smallest" "" Return Np.sum (yn-n P.convolve (x, h)) # **2 calls the Fmin function to x0 to the initial value H0 = Fminfunc (Convolve_func, x0) print fminfunc.__name__ print "---------------------" # outputs the relative error between x (*) H0 and y print "error of Y:", Np.sum (Np.convolve (x, H0)-y) **2)/np.sum (y*
    *2) # Output the relative error between H0 and H print "error of H:", Np.sum ((h0-h) **2)/np.sum (h**2) print Def test_n (M, N, Nscale):  "" "" "randomly produces x, h, y, yn, x0 and other series, call the various fmin functions to solve the length of b m x, n is the length of H, nscale the intensity of interference" "" x = Np.random.rand (m) h =

    Np.random.rand (n) y = np.convolve (x, h) yn = y + np.random.rand (len (y)) * Nscale x0 = Np.random.rand (n) Test_fmin_convolve (OPT.FMin, X, h, y, yn, x0) test_fmin_convolve (Opt.fmin_powell, X, h, y, yn, x0) test_fmin_convolve (OPT.FMIN_CG, X, H, y , yn, x0) test_fmin_convolve (OPT.FMIN_BFGS, X, h, y, yn, x0) if __name__ = "__main__": Test_n (200, 20, 0.1)

--------------------------------------------------------------------------------------------------------------- ----------------

if __name__ = = "__main__"

The function of the above code is to execute the statement below if the current file is compiled and run. If this file is called within another file through import, the IF statement is not executed.

--------------------------------------------------------------------------------------------------------------- ----------------
3. Solving nonlinear equations optimize the Fsolve function in the library can be used to solve the nonlinear equations. Its basic invocation form is as follows:

Fsolve (func, x0)

Func (x) is a function of calculating the error of the equations, and its parameter x is a vector representing a set of possible solutions for each unknown of the system of equations, Func returns the error resulting from the generation of x into the equations; x0 is the initial value of the unknown vector. If you want to solve the following equation group:

F1 (U1,U2,U3) = 0
F2 (u1,u2,u3) = 0
f3 (u1,u2,u3) = 0

Then func can be defined as follows:

def func (x):
    u1,u2,u3 = x return
    [F1 (U1,U2,U3), F2 (U1,U2,U3), F3 (U1,U2,U3)]

The following is a practical example that solves the solution of the following equation group:

5*x1 + 3 = 0
4*x0*x0-2*sin (x1*x2) = 0
x1*x2-1.5 = 0

The procedure is as follows:

From scipy.optimize import fsolve from
math import Sin,cos

def f (x):
    x0 = float (x[0))
    x1 = float (x[1)) 
  x2 = float (x[2]) return
    [
        5*x1+3,
        4*x0*x0-2*sin (x1*x2),
        x1*x2-1.5
    ] Result

= Fsolve ( F, [1,1,1])

print result
print f (Result)

The output is:

[ -0.70622057-0.6        -2.5       ]
[0.0, -9.1260332624187868e-14, 5.3290705182007514e-15]

Because the Fsolve function is an array of arguments passed when the function f is called, the calculation speed will be reduced if the elements in the array are evaluated directly, so here we first use the FLOAT function to convert the elements in the array to the standard floating-point number in Python. It then invokes the function in the standard math library to perform the operation.

When the offset group is resolved, the fsolve will automatically calculate the Jacobian matrix of the equations, if there are many unknowns in the equations, and when the Jacobian matrix is sparse, the transfer of a function that calculates the Jacobian matrix will greatly increase the speed of the computation.
Jacobian matrix
The Jacobian matrix is a matrix in which the first derivative is arranged in a certain way, it gives the optimal linear approximation of the differentiable equation and the given point, so it is similar to the derivative of the multivariate function. For example, the preceding function f1,f2,f3 and the u1,u2,u3 of unknown unknowns are as follows:

Using the fsolve example of the Jacobian matrix, the function J of the Jacobian matrix is passed to the Fsolve by the Fprime parameter, the function J and function f have an unknown solution vector parameter x, and function J computes the Jacobian matrix of the nonlinear equations at the vector x point. Because of the few unknowns in this example, the calculation of Jacobian matrices by programs does not lead to a rise in computational speed.

#-*-Coding:utf-8-*-from
scipy.optimize import fsolve from
math import Sin,cos
def f (x):
    x0 = float (x[ 0])
    x1 = float (x[1])
    x2 = float (x[2)) return
    [
        5*x1+3,
        4*x0*x0-2*sin (x1*x2),
        x1*x2-1.5
    ]
    
def j (x):
    x0 = float (x[0])
    x1 = float (x[1))
    x2 = float (x[2)) return    
    [
        [0, 5, 0],
        [8*x0, -2* X2*cos (X1*X2), -2*x1*cos (X1*X2)],
        [0, x2, X1]
    ] Result
 
= Fsolve (f, [1,1,1], fprime=j)
print result
print f (Result)

4. Numerical IntegrationNumerical integration is a numerical solution to a definite integral, for example, a numerical integral is used to calculate the area of a shape. Let's consider how to calculate the area of a semicircle radius of 1, according to the area of the circle formula, the area should be equal to PI/2. The unit semicircle (semicircle above the x-axis) curve can be represented by the following function:

def half_circle (x): Return
    (1-x**2) **0.5

The following procedure calculates the area of the unit semicircle by using the classic small rectangle to calculate the sum of the areas:

>>> N = 10000
>>> x = Np.linspace ( -1, 1, N)
>>> dx = 2.0/n # This value is actually the distance between each of the two adjacent points above X, that is, lifting Line width
>>> y = half_circle (x) # Each point of the ordinate, each point to the next extreme to see the high consistent. So the Y-value of each point times the DX is the area of the rectangle between the current point and the next.
>>> DX * np.sum (y) * 2 # twice times the area to see the PI value
3.1412751679988937

The coordinates of a series of points in a circle calculated by the above method can also be numerically integrated using Numpy.trapz:

>>> import NumPy as NP
>>> Np.trapz (y, x) * 2 # twice times
3.1415893269316042

This function calculates the area of a polyline and an X axis with x,y as the vertex coordinates. The same number of split points, the result of the Trapz function is closer to the exact value of some.

If we call the Quad function in the Scipy.integrate library, we will get very precise results:

>>> from scipy import integrate
>>> pi_half, err = Integrate.quad (half_circle,-1, 1)
>> > pi_half*2
3.1415926535897984

This section I moved over the content of a lot less, want to learn more content of the reader suggested to go to the original address to learn.
This article references from: HTTP://OLD.SEBUG.NET/PAPER/BOOKS/SCIPYDOC/SCIPY_INTRO.HTML#ID1