zoukankan      html  css  js  c++  java
  • Python机器学习(五十)NumPy 线性代数

    NumPy中包含了numpy.linalg模块,提供线性代数运算功能。下表描述了该模块中的一些重要功能。

    SN函数描述
    1 dot() 两个数组的点积
    2 vdot() 两个向量的点积
    3 inner() 两个数组的内积
    4 matmul() 两个数组的矩阵乘积
    5 det() 计算矩阵的行列式
    6 solve() 解线性矩阵方程
    7 inv() 求矩阵的乘法逆矩阵

    numpy.dot()

    numpy.dot()计算两个数组的点积。

    示例

    import numpy as np  
    a = np.array([[100,200],[23,12]])  
    b = np.array([[10,20],[12,21]]) 
    
    dot = np.dot(a,b)  #[100 * 10 + 200 * 12, 100 * 20 + 200 * 21] [23*10+12*12, 23*20 + 12*21]
    print(dot)  

    输出

    [[3400 6200]
     [ 374  712]]

    numpy.vdot()

    numpy.dot()计算两个向量的点积。

    示例

    import numpy as np  
    a = np.array([[100,200],[23,12]])  
    b = np.array([[10,20],[12,21]])  
    
    vdot = np.vdot(a,b)  # 100 *10 + 200 * 20 + 23 * 12 + 12 * 21
    print(vdot)  

    输出

    5528

    numpy.inner()

    numpy.inner()计算两个数组的内积。这个函数返回一维数组内部元素乘积的和。对于n维数组,它返回元素在最后一个轴上的乘积的和。

    示例

    import numpy as np  
    a = np.array([1,2,3,4,5,6])  
    b = np.array([23,23,12,2,1,2])  
    inner = np.inner(a,b)  # 1*23 + 2*23 + 3*12 + 4*2 + 5*1 + 6*2
    print(inner)  

    输出

    130

    numpy.matmul()

    numpy.matmul()计算两个数组的矩阵乘积。

    示例

    import numpy as np  
    a = np.array([[1,2,3],[4,5,6],[7,8,9]])  
    b = np.array([[23,23,12],[2,1,2],[7,8,9]])  
    mul = np.matmul(a,b)  
    print(mul) 

    输出

    [[ 48  49  43]
     [144 145 112]
     [240 241 181]]

    numpy.linalg.det()

    numpy.linalg.det()计算矩阵的行列式。

    矩阵的行列式可以用对角元素来计算。下面是2×2矩阵的行列式

    A B
    C D

    det()计算公式:AD – BC

    示例

    import numpy as np  
    a = np.array([[1,2],[3,4]])  
    print(np.linalg.det(a)) 

    输出

    -2.0000000000000004

    numpy.linalg.solve()

    numpy.linalg.solve()解线性矩阵方程。

    这个函数用于解线性方程,学过线性代数可知,线性方程可以表示为矩阵。

    如下面的方程:

    x + y = 20
    x - y = 10
    
    可以表示为3个矩阵的方程:
    
    矩阵a:
    1   1
    1   -1
    
    矩阵b:
    x
    y
    
    矩阵c:
    20
    10
    
    方程式可以表示为:
    a x b = c

    示例

    import numpy as np  
    a = np.array([[1,1],[1,-1]])  
    c = np.array([[20], [10]])  
    b = np.linalg.solve(a, c)
    
    print(b)  

    输出

    [[15.]
     [ 5.]]

    scipy非线性方程组求解
    求解线性方程组比较简单,只需要用到一个函数(scipy.linalg.solve)就可以了。比如我们要求以下方程的解,这是一个非齐次线性方程组:

    程序代码:

    import numpy as np
    from scipy.linalg import solve
    #输出系数矩阵
    a=np.array([[3,1,-2],[1,-1,4],[2,0,3]])
    #
    b=np.array([5,-2,2.5])
    #计算
    x=solve(a,b)
    #打印结果
    print(x)

    执行结果:

    C:Anaconda3python.exe "C:Program FilesJetBrainsPyCharm 2019.1.1helperspydevpydevconsole.py" --mode=client --port=59051
    import sys; print('Python %s on %s' % (sys.version, sys.platform))
    sys.path.extend(['C:\app\PycharmProjects', 'C:/app/PycharmProjects'])
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)]
    Type 'copyright', 'credits' or 'license' for more information
    IPython 7.12.0 -- An enhanced Interactive Python. Type '?' for help.
    PyDev console: using IPython 7.12.0
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)] on win32
    runfile('C:/app/PycharmProjects/ArtificialIntelligence/test.py', wdir='C:/app/PycharmProjects/ArtificialIntelligence')
    [0.5 4.5 0.5]

    简单的方程求解

    2x-4=0的一元一次方程求解

    from sympy import *
    x= symbols('x')
    print(solve(x*2-4,x))

    执行结果:

    C:Anaconda3python.exe "C:Program FilesJetBrainsPyCharm 2019.1.1helperspydevpydevconsole.py" --mode=client --port=59147
    import sys; print('Python %s on %s' % (sys.version, sys.platform))
    sys.path.extend(['C:\app\PycharmProjects', 'C:/app/PycharmProjects'])
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)]
    Type 'copyright', 'credits' or 'license' for more information
    IPython 7.12.0 -- An enhanced Interactive Python. Type '?' for help.
    PyDev console: using IPython 7.12.0
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)] on win32
    runfile('C:/app/PycharmProjects/ArtificialIntelligence/test.py', wdir='C:/app/PycharmProjects/ArtificialIntelligence')
    [2]

    需要说明的是:solve:第一个参数为要解的方程,要求右端等于0,第二个参数为要解的未知数。还有一些 其他的参数,想了解的可以去看官方文档。

    二元一次方程方程的求解

    求解方程:

    from sympy import *
    x,y= symbols('x,y')
    print(solve([2*x-y-3,3*x+y-7],[x,y]))

    执行结果:

    C:Anaconda3python.exe "C:Program FilesJetBrainsPyCharm 2019.1.1helperspydevpydevconsole.py" --mode=client --port=59201
    import sys; print('Python %s on %s' % (sys.version, sys.platform))
    sys.path.extend(['C:\app\PycharmProjects', 'C:/app/PycharmProjects'])
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)]
    Type 'copyright', 'credits' or 'license' for more information
    IPython 7.12.0 -- An enhanced Interactive Python. Type '?' for help.
    PyDev console: using IPython 7.12.0
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)] on win32
    runfile('C:/app/PycharmProjects/ArtificialIntelligence/test.py', wdir='C:/app/PycharmProjects/ArtificialIntelligence')
    {x: 2, y: 1}

    一元二次方程的求解
    求解:

    from sympy import *
    x= symbols('x')
    print(solve(x**2+2*x+1,x))

    执行结果:

    [-1]
    [Finished in 1.0s]

    求解八元一次方程

    from sympy import *
    
    a1,b1,c1,d1,a2,b2,c2,d2=symbols('a1,b1,c1,d1,a2,b2,c2,d2')
    x0,x1,x2,y0,y1,y2=symbols('x0,x1,x2,y0,y1,y2')
    result=solve([a1+b1*x0+c1*x0*x0+d1*x0*x0*x0-y0,
        a1+b1*x1+c1*x1*x1+d1*x1*x1*x1-y1,
        a2+b2*x1+c2*x1*x1+d2*x1*x1*x1-y1,
        a2+b2*x2+c2*x2*x2+d2*x2*x2*x2-y2,
        b1+2*c1*x1+3*d1*x1*x1-b2-2*c2*x1-3*d2*x1*x1,
        c1+3*d1*x1-c2-3*d2*x1,
        c1+3*d1*x0,
        c2+3*d2*x2],[a1,b1,c1,d1,a2,b2,c2,d2])
    #设置参数
    cs=[(x0,6),(x1,9),(x2,12),(y0,0),(y1,3),(y2,0)]
    
    # print(result[a1])
    # print(result[b1])
    # print(result[c1])
    # print(result[d1])
    # print(result[a2])
    # print(result[b2])
    # print(result[c2])
    # print(result[d2])
    
    a1=result[a1].subs(cs)
    b1=result[b1].subs(cs)
    c1=result[c1].subs(cs)
    d1=result[d1].subs(cs)
    a2=result[a2].subs(cs)
    b2=result[b2].subs(cs)
    c2=result[c2].subs(cs)
    d2=result[d2].subs(cs)
    print(a1)
    print(b1)
    print(c1)
    print(d1)
    print(a2)
    print(b2)
    print(c2)
    print(d2)

    执行结果

    C:Anaconda3python.exe "C:Program FilesJetBrainsPyCharm 2019.1.1helperspydevpydevconsole.py" --mode=client --port=59277
    import sys; print('Python %s on %s' % (sys.version, sys.platform))
    sys.path.extend(['C:\app\PycharmProjects', 'C:/app/PycharmProjects'])
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)]
    Type 'copyright', 'credits' or 'license' for more information
    IPython 7.12.0 -- An enhanced Interactive Python. Type '?' for help.
    PyDev console: using IPython 7.12.0
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)] on win32
    runfile('C:/app/PycharmProjects/ArtificialIntelligence/test.py', wdir='C:/app/PycharmProjects/ArtificialIntelligence')
    3
    -9/2
    1
    -1/18
    -78
    45/2
    -2
    1/18

    Python求解各种复杂的线性/非线性方程组
    本文将要介绍几种方法去求解各种复杂的方程组,包括实数域和复数域的线性、非线性方程组,并对比这几种方法的优缺点。本文用到了numpy、scipy、sympy这三个科学计算包。

    一、线性方程组。

    线性方程组可以用numpy去求解。
    1.实数域。

    import numpy as np
    a=np.mat('1,2,3;2,4,8;9,6,3')
    b=np.mat('1;1;3')
    c=np.linalg.solve(a,b)

    输出如下:

    2.复数域。

     

    import numpy as np
    a=np.mat('1,-1j;1j,-1')
    b=np.mat('1;1')
    c=np.linalg.solve(a,b)

    输出如下:

    二、非线性方程组。

    scipy和sympy不但可以解线性方程(组),还可以求解非线性方程(组),但是也有各自的优缺点。

    1.scipy求解
    scipy.optimize里面有两个函数可以数值求解方程组,分别是root和solve,这两个函数会找到方程组的一个近似解。下面通过例子介绍这两个函数的使用。

    (1).root
    scipy.optimize.root(fun, x0, args=(), method=‘hybr’, jac=None, tol=None, callback=None, options=None)
    其中参数fun是一个函数,你需要定义这个函数,并且你定义的这个函数返回一个方程组。 参数x0是方程组解的初始猜测值,是必须要指定的参数,root函数将找到最靠近x0的一个解。参数args是函数fun的额外参数。参数method默认是‘hybr’,此参数还可以取 ‘lm’ ‘broyden1’ ‘broyden2’ ‘anderson’ ‘linearmixing’ ‘diagbroyden’ ‘excitingmixing’ ‘krylov’ ‘df-sane’,代表采用不同的算法去求解。参数tol可以认为是精度。

    先来看看实数域的例子

    from scipy.optimize import root
    
    def f1(x):
    return [x[0]+x[0]*x[1]-2,x[0]-x[1]-2]
    
    print(root(f1,[0,-1]).x)#初始猜测值[0,-1]
    print(root(f1,[0,0]).x)#初始猜测值[0,0]

    输出如下:

    求出的近似解与真解相差非常小。设置不同的初始猜测值可能得到不同的解,这是因为有的方程组的解不一定只有一个,可能有多个解,root函数会得到最接近初始猜测值的一个解。

    此外,在上面的基础上我们可以设置参数jac来提高运算速度(尤其是计算量很大时效果很明显)。
    我们需要再定义一个函数,这个函数返回值是方程组对应的雅可比矩阵,通过下面的例子说明。

    from numpy import array,mat,sin,cos,exp
    from scipy.optimize import root
    
    def f(x):
    eqs=[]
    eqs.append(x[0]*x[1]+x[1]*x[2]+sin(x[0])*exp(x[1])+x[1])
    eqs.append(x[0]*x[1]-exp(x[0])*x[1]+1)
    eqs.append(x[1]*x[2]+exp(x[1])*x[2]+1) 
    return eqs
    
    def jac1(x):#方程组对应的雅可比矩阵
    return mat([[x[1]+cos(x[0])*exp(x[1]), x[0]+x[2]+sin(x[0])*exp(x[1])+1, x[1]],
    [x[1]-exp(x[0])*x[1], x[0]-exp(x[0]), 0],
    [0 ,x[2]+exp(x[1])*x[2], x[1]+exp(x[1])]])
    
    print(root(f,[0,0,0]).x)
    print(root(f,[0,0,0],jac=jac1).x)#加上参数jac加快运算速度

    输出如下:

    再来看一个复数域的例子

    如果不设置参数method,它默认是‘hybr’,即用改进的Powell hybrid算法求解,它不适用于复数域。此时我们不能再用上面的办法,我们可以设置参数method=‘krylov’ (此算法牺牲精度换取速度,有利有弊)

    from scipy.optimize import root
    
    def f1(x):
    return [x[0]*(1j)+x[0]*x[1]+1,x[0]+x[1]-1j]
    
    print(root(f1,[1,1],method='krylov').x)
    print(root(f1,[1,1],method='krylov',tol=1e-10).x)#设置能够允许的误差为10的-10次方

    输出如下:

    第一行输出没有设置参数tol,得到的解与真解的差别还是有一点的。
    第二行输出我们设置了参数tol是10的-10次方,得到的解与真解的差别就相当小了。
    另外method参数也可以设为别的,解复数域的非线性方程组建议用method=‘krylov’,它不需要雅可比矩阵,它在处理很庞大、变量很多的方程组时比较好用。

    (2).fslove
    fsolve的用法和root很类似,但是它的功能不如root全面,fsolve其实就是用hybr算法求解, 因此它不能解复数域的方程组。下图可以看出fsolve和root的差别:

    可见,fsolve只是root的一小部分。

    scipy.optimize.fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)
    (其中参数fprime对应雅可比矩阵,其余参数和root函数基本一样。)

    from scipy.optimize import fsolve
    from numpy import array,mat
    
    def f1(x):
    return [x[0]+x[0]*x[1]-2,x[0]-x[1]-2]
    
    def jac1(x):#方程组对应的雅可比矩阵
    return mat([[1+x[1],x[0]],[1,-1]])
    
    print(fsolve(f1,[0,-1]))#初始猜测值[0,-1]
    print(fsolve(f1,[0,-1],fprime=jac1))#初始猜测值[0,-1],并设置参数prime
    print(fsolve(f1,[0,0]))#初始猜测值[0,0]
    print(fsolve(f1,[0,0],fprime=jac1))#初始猜测值[0,0],并设置参数prime

    2.sympy求解
    sympy中的solve函数可以严格求解某些方程组,nsolve可以数值求近似解。
    (1).solve
    直接看一个复数域的例子:

    from sympy import symbols,Eq,solve
    x0,x1=symbols('x0 x1')
    eqs=[Eq(x0*(1j)+x0*x1,-5),Eq(x0+x1,1j)]
    print(solve(eqs,[x0,x1]))

    输出如下:

    solve函数就会找到所有的解。但是像超越方程之类的不存在求根公式的方程,solve函数是不能求解的,只能数值求解,要用nsolve函数。
    (2).nsolve
    nsolve函数需要指定一个初始猜测解。

    from sympy import symbols,Eq,nsolve
    x0,x1=symbols('x0 x1')
    eqs=[Eq(x0*(1j)+x0*x1,-5),Eq(x0+x1,1j)]
    print(nsolve(eqs,[x0,x1],[1,1]))#初始猜测解设为[1,1]

    输出如下:

    nsolve有时候并不是很好使,初始猜测解设的不好,它有可能找不到解。

    三、scipy和sympy的优缺点分析。
    1.scipy.optimize.root求解方程组速度很快,尤其是加上参数jac或参数method='krylov’时,求解大型方程组速度会很明显的比其它办法快,即使是1000个变量的方程组,它也能很快解完。但是有的方程组有多解(比如二次方程有俩根),而scipy.optimize.root只能得到靠近初始猜测值的那个解。
    2.对一些形式比较简单的、有求根公式的方程,sympy.solve虽然能得到它所有的严格解,但是当方程组变量较多时,它求起来会很慢。而且对于不存在求根公式的复杂方程,sympy.solve是求不了的,这时要用sympy.nsolve的求数值解,速度也比scipy.optimize.root慢很多。

    四、总结
    1.线性方程组用numpy.linalg.solve足矣。(实数域复数域都可以)
    2.求解非线性方程组,如果存在求根公式的且未知量数目较少的方程组,可以用sympy求。如果是很庞大且形式较复杂的方程组,用scipy.optimize.root数值求解,最好写出雅可比矩阵以提高运算速度。(如果你懒得写雅可比矩阵,就设置参数method=‘krylov’,解起来速度也很快,只不过精度不高)。

    最后展示一段代码,测试一下scipy.optimize.root的求解速度,这里是求解含有400个未知量的非线性方程组,解20次,耗时2秒。

    from numpy import array,arange,sin,cos,sqrt,exp,pi,meshgrid,dtype,linalg,zeros
    from scipy.optimize import root,fsolve
    from matplotlib.pyplot import figure,xlabel,ylabel,show,text
    from mpl_toolkits.mplot3d import Axes3D
    from time import time
    
    def run(m,n,marker=None):
    v=1
    L,T=pi,1
    h,tao=L/m,T/n
    mlist=arange(0,m+1)
    nlist=arange(0,n+1)
    f=lambda x,t:exp(-2*t)*sin(x)*cos(x)
    phi=lambda x:sin(x)
    realsolve=lambda x,t:exp(-t)*sin(x)
    mysolves=[phi(mlist*h)]
    
    def u(x,ii,j):
    if j==1:
    return x[ii]
    elif j==1/2:
    return (x[ii]+x_former[ii])/2
    elif j==0:
    return x_former[ii]
    else:
    print("error")
    
    def jac(x):#雅可比矩阵
    a=zeros((m+1,m+1))
    a[0][0]=1
    a[-1][-1]=1
    for i in range(1,m):
    a[i][i-1]=1/6/h*((x[i+1]+x_former[i+1])/2-(x[i-1]+x_former[i-1])/2-((x[i-1]+x_former[i-1])/2+(x[i]+x_former[i])/2+(x[i+1]+x_former[i+1])/2))-v/(2*h*h)
    a[i][i]=1/tao+1/6/h*((x[i+1]+x_former[i+1])/2-(x[i-1]+x_former[i-1])/2)+v/(h*h)
    a[i][i+1]=1/6/h*((x[i+1]+x_former[i+1])/2-(x[i-1]+x_former[i-1])/2+((x[i-1]+x_former[i-1])/2+(x[i]+x_former[i])/2+(x[i+1]+x_former[i+1])/2))-v/(2*h*h)
    
    return a
    
    def eqsgroup(x):
    eqs=[x[0]]
    for i in range(1,m):
    eqs.append(1/tao*(u(x,i,1)-u(x,i,0))+1/(6*h)*(u(x,i-1,1/2)+u(x,i,1/2)+u(x,i+1,1/2))*(u(x,i+1,1/2)-u(x,i-1,1/2))-v/h/h*(u(x,i-1,1/2)-2*u(x,i,1/2)+u(x,i+1,1/2))-f(i*h,k*tao) )#
    eqs.append(x[m])
    return array(eqs)
    
    x_former=phi(mlist*h)
    start=time()
    for k in arange(0,n):
    #x_former=fsolve(eqsgroup,x_former,fprime=jac)
    x_former=root(eqsgroup,x_former,jac=jac).x
    mysolves.append(list(x_former))
    print('scipy计算时间为',time()-start,'')
    
    x,y=mlist*h,nlist*tao
    x,y=meshgrid(x, y)
    z=array(mysolves,dtype=dtype('float32'))
    z2=array(realsolve(x,y))
    dz=z-z2
    ekmax=array(linalg.norm(dz,axis=1)*sqrt(h)).max()
    
    fig = figure(figsize=(15, 10),dpi=80)
    axes3d = Axes3D(fig)
    axes3d.set_title(r'$h=pi/$'+str(m)+r' $	au=1/$'+str(n),color = 'black',size = 30) 
    axes3d.plot_surface(x,y,z)
    axes3d.plot_surface(x,y,z2)
    xlabel('x',color = 'black',size = 30)
    ylabel('t',color = 'black',size = 30)
    axes3d.set_xlim(0,pi)
    axes3d.set_xticks([0,pi/4,pi/2,0.75*pi,pi])
    axes3d.set_xticklabels([0, r'$frac{pi}{4}$', r'$frac{pi}{2}$', r'$frac{3pi}{4}$', r'$pi$'],size=20)
    axes3d.set_zlabel(r'$U(x,t)$',color = 'black',size=30)
    axes3d.set_zticks(arange(0,1,0.1))
    axes3d.set_zticklabels([0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9],size=13)
    show()
    
    
    run(400,20)

    如果用sympy去解,等很久都解不出。scipy解非线性方程组的能力不比matlab差。

    矩阵的特征值和特征向量

    import numpy as np
    A = np.array([[3,-1],[-1,3]])
    print('打印A:
    {}'.format(A))
    a, b = np.linalg.eig(A)
    print('打印特征值a:
    {}'.format(a))
    print('打印特征向量b:
    {}'.format(b))
    

    执行结果:

    C:Anaconda3python.exe "C:Program FilesJetBrainsPyCharm 2019.1.1helperspydevpydevconsole.py" --mode=client --port=54930
    import sys; print('Python %s on %s' % (sys.version, sys.platform))
    sys.path.extend(['C:\app\PycharmProjects', 'C:/app/PycharmProjects'])
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)]
    Type 'copyright', 'credits' or 'license' for more information
    IPython 7.12.0 -- An enhanced Interactive Python. Type '?' for help.
    PyDev console: using IPython 7.12.0
    Python 3.7.6 (default, Jan  8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)] on win32
    runfile('C:/app/PycharmProjects/ArtificialIntelligence/test.py', wdir='C:/app/PycharmProjects/ArtificialIntelligence')
    打印A:
    [[ 3 -1]
     [-1  3]]
    打印特征值a:
    [4. 2.]
    打印特征向量b:
    [[ 0.70710678  0.70710678]
     [-0.70710678  0.70710678]]

    一 要点

    • 假定AX=b,求解未知矩阵X 【线性代数中常遇到的运算问题】
    • 矩阵转置A^(T)
    • 矩阵的逆A^(-1)
    • 矩阵行列式的值|A|
    • 矩阵的秩 rank(A)
    • 矩阵的迹 trace(A)
    • 其它
      • 单位矩阵
      • 0向量/矩阵
      • ...

    二 示例应用

    2.1 求解AX=b中的未知参数矩阵X

    import numpy as np
    
    # Hypothsis : A*X = b
    
    A = [[2,1,2],
         [3,1,0],
         [1,1,-1]];
    b = np.transpose([-3,5,-2])# 转置
    #[or] b = np.transpose(np.array([-3,5,-2]))# 转置
    
    # 求解未知参数矩阵X
    X = np.linalg.solve(A,b) # 方式一:直接使用numpy的solve函数一键求解
    #A_inv=np.linalg.inv(A) # 方式二:先求逆运算,再点积求值
    #X=np.dot(A_inv,b) # a.dot(b) 与 np.dot(a,b) 效果相同;but np.dot(a,b)与np.dot(b,a)效果肯定是不同的(线性代数/矩阵常识)
    print("方程组的解:
    ",X);
    # [output]
    方程组的解:
     [ 4.4 -8.2 -1.8]

    2.2 利用最小二乘法拟合函数模型

    给出一组数据【5对(Xi,Yi)参数】,用最小二乘法,求形如:f(x)=a+b*x^3的经验公式。

    • 原方程(模型)
    f(x)=a+bx3f(x)=a+bx3
    • 其法方程
    ATA(ab)=ATyATA(ab)=ATy

    (ab)=(ATA)1ATy(ab)=(ATA)−1ATy
    import numpy as np
    import matplotlib.pyplot as plt  # Python 绘图工具(业界推荐)
    
    # 数据初始化
    A = [
        [1,pow(-3,3)],
        [1,pow(-2,3)],
        [1,pow(-1,3)],
        [1,pow(2,3)],
        [1,pow(4,3)]
    ];
    
    At =  np.transpose(A); # A的转置矩阵
    
    y = np.transpose([14.3,8.3,4.7,8.3,22.7]);
    
    
    # step1:求解
    ## 令 (a ,b)^T 为 未知参数X
    X = np.dot(np.dot(np.linalg.inv(np.dot(At,A)),At),y)
    print(X)
    print("a:",X[0])
    print("b:",X[1])
    
    # step2:查看拟合效果
    x = [-3,-2,-1,2,4];
    ## 1. 计算拟合数值 fitValue
    def fitValue(arg_x):
        a = X[0];
        b = X[1];
        return a + b*pow(arg_x,3);
    
    fitValues = np.zeros([5]); # 创建长为5的【一维】数组;[1,5]:创建第1行为5个元素的【二维】数组 
    for i in range(0,len(fitValues)):
        fitValues[i] = fitValue(x[i]);
        print(i,":","x(i):",x[i],"fit Y:",fitValues[i]); # just for test
        pass;
    
    ## 2. 绘图可视化
    yt = np.transpose(y); # y的转置
    plt.rcParams['figure.dpi'] = 100 #分辨率
    plt.scatter(x, yt, marker = '*',color = 'red', s = 10 ,label = 'Actual Dataset') # 真实数据集
    plt.scatter(x, fitValues, marker = 'x',color = 'green', s = 10 ,label = 'Fitting Dataset') #[拟合数据集]
    plt.legend(loc = 'best')    # 设置 图例所在的位置 使用推荐位置
    plt.show() 
    
    
    # [output]
    [ 10.67505325   0.13679816]
    a: 10.6750532504
    b: 0.136798159666
    0 : x(i): -3 fit Y: 6.98150293942
    1 : x(i): -2 fit Y: 9.58066797308
    2 : x(i): -1 fit Y: 10.5382550907
    3 : x(i): 2 fit Y: 11.7694385277
    4 : x(i): 4 fit Y: 19.430135469

    2.3 数组创建/初始化

    • numpy.linspace(start, stop[, num=50[, endpoint=True[, retstep=False[, dtype=None]]]]])
      • 一维等差数列
      • return 在指定范围内的均匀间隔的数字(组成的数组),也即返回一个等差数列
      • 参数
        • start - 起始点
        • stop - 结束点
        • num - 元素个数,默认为50,
        • endpoint - 是否包含stop数值,默认为True,包含stop值;若为False,则不包含stop值
        • retstep - 返回值形式,默认为False,返回等差数列组,若为True,则返回结果(array([samplesstep])),
        • dtype - 返回结果的数据类型,默认无,若无,则参考输入数据类型。
    import numpy as np
    
    a = np.linspace(1,10,5,endpoint= True)
    print(a) # [ 1.    3.25  5.5   7.75 10.  ]
    
    b = np.linspace(1,10,5,endpoint= False)
    print(b) #[1.  2.8 4.6 6.4 8.2]
    
    c = np.linspace(1,10,5,retstep = False)
    print(c) # [ 1.    3.25  5.5   7.75 10.  ]
    
    d = np.linspace(1,10,5,retstep = True)
    print(d) # (array([ 1.  ,  3.25,  5.5 ,  7.75, 10.  ]), 2.25)

    2.4 线性代数常用运算

    print("原矩阵A:
    ",A);
    print("原矩阵b:
    ",b);
    
    print("转置矩阵A^T:
    ",np.transpose(A)); # 转置
    print("矩阵的行列式值|A|:
    ",np.linalg.det(A)); # 方阵的行列式值:|A|
    print("矩阵的迹trace(A):
    ",np.trace(A)); 
    print("矩阵的秩rank(A):
    ",np.linalg.matrix_rank(A)); 
    print("逆矩阵A^(-1):
    ",np.linalg.inv(A)); #矩阵的逆运算(条件:矩阵A可逆(行列式值不为0)| 矩阵A为方阵)
    
    print("*"*30); # 分隔线
    
    print("N阶单位矩阵:
    ",np.eye(4));
    print(np.zeros([5])); # 创建长为5的【一维】数组;[1,5]:创建第1行为5个元素的【二维】数组 
    
    # 创建指定的初始化数组
    print(np.array([1]))  # 生成 第1行含值为1的元素的【一维】数组
    print(np.array([[56]])) # 生成 第1行含值为56的元素的【二维】数组
    np.full((3,5),3.14) # 创建一个3x5的浮点型数组,数组的值都是3.14
    # [output]
    原矩阵A:
     [[2, 1, 2], [3, 1, 0], [1, 1, -1]]
    原矩阵b:
     [-3  5 -2]
    转置矩阵A^T:
     [[ 2  3  1]
     [ 1  1  1]
     [ 2  0 -1]]
    矩阵的行列式值|A|5.0
    矩阵的迹trace(A):
     2
    矩阵的秩rank(A):
     3
    逆矩阵A^(-1):
     [[-0.2  0.6 -0.4]
     [ 0.6 -0.8  1.2]
     [ 0.4 -0.2 -0.2]]
    **********************
    N阶单位矩阵:
     [[ 1.  0.  0.  0.]
     [ 0.  1.  0.  0.]
     [ 0.  0.  1.  0.]
     [ 0.  0.  0.  1.]]
    [ 0.  0.  0.  0.  0.]
    [1]
    [[56]]
     
  • 相关阅读:
    ORACLE触发器具体解释
    秒杀多线程第四篇 一个经典的多线程同步问题
    Java中Integer类的方法
    九大排序算法再总结
    删除条目时的确认对话框
    VirtualBox安装及使用说明和虚拟机安装XP系统图文教程
    J2EE之验证码实现
    教你用笔记本破解无线路由器password
    vSphere HA状况:未知配置错误解决的方法
    HDU 2504 又见GCD
  • 原文地址:https://www.cnblogs.com/huanghanyu/p/13163037.html
Copyright © 2011-2022 走看看