[应用多元统计分析-王学民]——第一章矩阵代数Python

OLSRR

已于 2022-11-10 18:08:05 修改

阅读量1.4k

点赞数 4

分类专栏： Python 统计文章标签： python 矩阵

于 2022-10-26 16:08:05 首次发布

本文链接：https://blog.csdn.net/m0_57096988/article/details/127533684

版权

统计同时被 2 个专栏收录

9 篇文章 2 订阅

订阅专栏

Python

5 篇文章 2 订阅

订阅专栏

以下##均为输出结果

1.1定义

矩阵

import numpy as np
a = np.arange(9).reshape((3, 3))
print(a)

##[[0 1 2]
## [3 4 5]
## [6 7 8]]

列向量

x1 = np.array([[6],[9],[3]])
print(x1)

##[[6]
## [9]
## [3]]

行向量

x2 = np.array([[2,7,3,4,3]])
print(x2)
##[[2 7 3 4 3]]

向量长度

a =np.arange(10)
print(a)
b=input(a.size)
print(b)
##[0 1 2 3 4 5 6 7 8 9]
##10
arr_std = np.std(a,ddof=1)
print(arr_std)
##3.0276503540974917

范数, 欧式范数, 内积开根号(相等)

内积空间 ⊂ 赋范空间 ⊂ 距离空间 ⊂ 拓扑空间

内积诱导范数，范数诱导距离，距离诱导拓扑

向量几何

向量（1，1）

#准备工作
import matplotlib.pyplot as plt
import mpl_toolkits.axisartist as axisartist
fig=plt.figure(figsize=(4,2)) #新建画布
ax=axisartist.Subplot(fig,111)
fig.add_axes(ax)
#隐藏矩形
ax.axis[:].set_visible(False)
#坐标轴设置
ax.axis["x"]=ax.new_floating_axis(0,0,axis_direction="bottom")
ax.axis["y"]=ax.new_floating_axis(1,0,axis_direction="bottom")

ax.axis["x"].set_axisline_style("->",size=1.0)
ax.axis["y"].set_axisline_style("->",size=1.0)
ax.annotate(s='x' ,xy=(2,0) ,xytext=(2,0.1))
ax.annotate(s='y' ,xy=(0,1.0) ,xytext=(-0.5,1.0))

plt.xlim(-2,2)
plt.ylim(-2,2)
ax.set_xticks([-2,-1,0,1,2])
ax.set_yticks([-1,1])
#向量设置
plt.arrow(0,0,1,1,head_width=0.06)
plt.show()

零向量

array = np.zeros((2,3))
print(array)
#[[0. 0. 0.]
#[0. 0. 0.]]

对角线元素

a = np.arange(9).reshape((3, 3))
print(a)
##[[0 1 2]
## [3 4 5]
## [6 7 8]]
print(np.diag(a))
##[0 4 8]

非对角线元素

a = np.arange(9).reshape((3, 3))
b = a.flatten()
c=[]
for i in b:
    if i ==0:
        continue
    elif i ==4:
        continue
    elif i ==8:
        continue
    else:
        c.append(i)
print(c)
##[1, 2, 3, 5, 6, 7]

上三角阵

m = np.arange(9).reshape((3, 3))
a = np.array(m)
a1 = np.triu(a,0)
print(a1)
##[[0 1 2]
## [0 4 5]
## [0 0 8]]

下三角阵

m = np.arange(9).reshape((3, 3))
a = np.array(m)
a2 = np.tril(a,0)
print(a2)
##[[0 0 0]
## [3 4 0]
## [6 7 8]]

对角阵矩阵

m = np.arange(9).reshape((3, 3))
a = np.diag(m)
print("\n",a)
##[[0 0 0
## [0 4 0]
## [0 0 8]]

单位矩阵

e1 = np.eye(5)
print("\n",e1)
## [[1. 0. 0. 0. 0.]
## [0. 1. 0. 0. 0.]
## [0. 0. 1. 0. 0.]
## [0. 0. 0. 1. 0.]
## [0. 0. 0. 0. 1.]]

矩阵转置

m = [[1,2,3],[4,5,6]]
a = np.array(m)
a1 = a.T
print(" \n",a1)
## [[1 4]
## [2 5]
## [3 6]]

对称矩阵

M = np.array([[2, 3, 4], [3, 45, 8], [34, 7, 0.8], [21, 31, 41]])
print(M.T)
S = np.matmul(M, M.T)
print(' \n', S)
##[[ 2.   3.  34.  21. ]
## [ 3.  45.   7.  31. ]
## [ 4.   8.   0.8 41. ]]
 
## [[  29.    173.     92.2   299.  ]
## [ 173.   2098.    423.4  1786.  ]
## [  92.2   423.4  1205.64  963.8 ]
## [ 299.   1786.    963.8  3083.  ]]

1.2矩阵的计算

例题1.2.1

m1 = [[1,4,2],[3,1,5]]
a1 = np.array(m1)
print(a1)
##[[1 4 2]
## [3 1 5]]

m2 = [[6,0,2],[3,1,4]]
a2 = np.array(m2)
print(a2)
##[[6 0 2]
## [3 1 4]]

m3 = [[1,1],[2,3]]
a3 = np.array(m3)
print(a3)
##[[1 1]
## [2 3]]

转置（A）和（B）

a11 = a1.T
print(a11)
##[[1 3]
## [4 1]
## [2 5]]

a21 = a2.T
print(a21)
##[[6 3]
## [0 1]
## [2 4]]

计算A+B;计算转置A+转置B;矩阵相乘（A*C）

#计算A+B
print("a1 + a2 = \n",a1 + a2)
##a1 + a2 = 
## [[7 4 4]
## [6 2 9]]


#计算转置A+转置B
print("a11 + a21 = \n",a11 + a21)
##a11 + a21 = 
## [[7 6]
## [4 2]
## [4 9]]

#矩阵相乘（A*C）
a4 = np.dot(a3,a1)
print("矩阵乘法a3×a1 = \n",a4)
##矩阵乘法a1×a3 = 
## [[ 4  5  7]
## [11 11 19]]

正交矩阵:

若方阵 A 满足 AA' = I =A'A, 则称 A 为正交矩阵。

m1 = [[1,0],[0,-1]]
a1 = np.array(m1)
A = ((math.sqrt(2)/2,-math.sqrt(2)/2,math.sqrt(2)/2,math.sqrt(2)/2))
A = np.array(A)
print(A)
##[[ 0.70710678 -0.70710678
##  0.70710678  0.70710678]]

m1 = [[1,0],[0,-1]]
a1 = np.array(m1)
a2 = np.array(m2)
a11 = a1.T
a3 = np.dot(a11,a1)
print(a3)
##[[1 0]
## [0 1]]

1.3行列式

例1.3.1（行列式）

m1 = [[3,1],[2,6]]
a = np.array(m1)
result = np.linalg.det(a)
print(result)
##16

例1.3.2（行列式性质）

m1 = [[3,1],[2,6]]
a = np.array(m1)
print(a)
##[[3 1]
## [2 6]]

m2 = [[1,2],[3,4]]
b = np.array(m2)
print(a)
##[[1 2]
##[3 4]]


#矩阵相乘再行列式计算

m1 = [[3,1],[2,6]]
a = np.array(m1)
m2 = [[1,2],[3,4]]
b = np.array(m2)
c = np.dot(a,b)
result = np.linalg.det(c)
print(result)
##-31.999999999999957

#直接进行行列式相乘后再相乘
m1 = [[3,1],[2,6]]
a = np.array(m1)
result1 = np.linalg.det(a)
m2 = [[1,2],[3,4]]
b = np.array(m2)
result2 = np.linalg.det(b)
print(result1*result2)
##-32.000000000000

例1.3.4（代数余子式）

m1 = [[6,5,10,1],[7,10,7,6],[9,8,12,2],[4,9,11,3]]
a = np.array(m1)
print(a)
##[[ 6  5 10  1]
## [ 7 10  7  6]
## [ 9  8 12  2]
## [ 4  9 11  3]]

n_array = np.array([[6,5,10,1],[7,10,7,6],[9,8,12,2],[4,9,11,3]])
row = 1 #选取行数
def count(n_array,i_row):
    n = n_array.shape[0]
    def count_inner(n_array,i_row):
        sum_part = 0
        if n_array.shape[0] ==1:
            return n_array[0][0]
        for i_col in range(n_array.shape[0]):
            n_array_next = np.delete(n_array,i_row-1,axis=0)
            n_array_next = np.delete(n_array_next,i_col,axis=1)
            a = n_array[i_row - 1][i_col]
            A = (count_inner(n_array_next,i_row)*((-1)**(i_row+i_col+1)))
            D =  a*A
            sum_part += D
            if n_array.shape[0] ==n:
                print ("sum_part",sum_part)
        return sum_part
    return count_inner(n_array,i_row)
print(count(n_array,row))
##280

1.4矩阵的逆

例1.4.1

a = np.array([[2,-1],[3,6]])
print(a)
##[[ 2 -1]
## [ 3  6]]

a = np.array([[2,-1],[3,6]])
result = np.linalg.det(a)
print(result)
##15.0

#伴随矩阵
def remove_line(index, arr):
    yy = []
    for i in range(len(arr)):
        hh = []
        for j in range(len(arr)):
            if (i == index[0] or j == index[1]):
                pass
            else:
                hh.append(arr[i][j])
        if (hh != []):
            yy.append(hh)
    # yy.remove([])
    return yy


dd = [[2, -1],[3,6]]

for i in range(len(dd)):
    for j in range(len(dd)):
        tm_list = remove_line(tuple([i, j]), dd)
        print((i, j), tm_list)
##(0, 0) [[6]]
##(0, 1) [[3]]
##(1, 0) [[-1]]
##(1, 1) [[2]]

例1.4.2

a = np.array([[2,-1,0],[3,6,0],[0,0,-4]])
##[[ 2 -1  0]
## [ 3  6  0]
## [ 0  0 -4]]

from fractions import Fraction
import numpy as np
np.set_printoptions(formatter={'all':lambda x: str(Fraction(x).limit_denominator())})
a = np.array([[2,-1,0],[3,6,0],[0,0,-4]])
b = np.linalg.inv(a)
print("分数形式显示逆矩阵：")
print(b)
##分数形式显示逆矩阵：
##[[2/5 1/15 0]
## [-1/5 2/15 0]
## [0 0 -1/4]]

1.5矩阵的秩

a = np.array([[1,2,3],[4,-2,5]])
print(a)
##[[ 1  2  3]
## [ 4 -2  5]]

m = np.array([[1,2,3],[4,-2,5]])
a = np.array(m)
rank = np.linalg.matrix_rank(m)
print(rank)
##2

1.6 特征值、特征向量和矩阵的迹

1.6.1

m = np.array([[3,2],[2,6]])
print(m)
##[[3 2]
## [2 6]]

a2=np.array([[3,2],[2,6]])  
b2=np.array([[3,2],[2,6]],dtype=int)  
c = np.linalg.eig(a2)     #返回矩阵a2的特征值与特征向量
print(c)
##(array([2., 7.]), array([[-0.89442719, -0.4472136 ],
##       [ 0.4472136 , -0.89442719]]))

例1.6.3

a=np.array([[2],[-4],[1]])
print(a)
##[[ 2]
## [-4]
## [ 1]]

b = np.array([[3],[5],[-1]])
print(b)
##[[ 3]
## [ 5]
## [-1]]

bt = b.T
c = np.dot(a,bt)
d = np.linalg.eig(c)
print(d)
##(array([-1.50000000e+01,  2.78104093e-15,  1.44497283e-17]), array([[ 4.36435780e-##01, -8.47157329e-01, -1.69752855e-04],
##       [-8.72871561e-01,  4.52636702e-01,  1.96214066e-01],
##       [ 2.18217890e-01, -2.78288477e-01,  9.80561070e-01]]))

at = a.T
c = np.dot(at,b)
d = np.linalg.eig(c)
print(d)
##(array([-15.]), array([[1.]]))

例1.6.5

a=np.array([[1,2,2],[2,1,2],[2,2,1]])
print(a)
##[[1 2 2]
## [2 1 2]
## [2 2 1]]

b = np.linalg.eig(a)
print(b)
##(array([-1.,  5., -1.]), array([[-0.81649658,  0.57735027,  0.        ],
##       [ 0.40824829,  0.57735027, -0.70710678],
##       [ 0.40824829,  0.57735027,  0.70710678]]))

例1.6.6

a=np.array([[1,2,2],[1,-2,2]])
print(a)
##[[ 1  2  2]
## [ 1 -2  2]]

def xsvd(H):
    U, S, V = np.linalg.svd(H)
    print('U:', U)
    print('S:', S)
    print('V:', V)
    return U, S, V


def UxSxV(U, S, V):
    S1 = np.array([[S[0], 0., 0.],
                   [0., S[1], 0.],
                   ])
    m1 = np.dot(U, S1)
    m2 = np.dot(m1, V)
    print(m2)
    return m2


def MxMtran(prom, U):
    print(prom)
    Ut = np.transpose(U)
    print(Ut)
    I = np.dot(U, Ut)
    print(I)


def main():
    H = np.array([[1., 2., 2.],
                  [1., -2.,2]]
                 )
    u, s, v = xsvd(H)
    h = UxSxV(u, s, v)
    print(h - H < 0.00001)
    MxMtran('U:', u)
    MxMtran('V:', v)


main()
##U: [[-0.70710678 -0.70710678]
## [-0.70710678  0.70710678]]
##S: [3.16227766 2.82842712]
##V: [[-4.47213595e-01 -4.15130595e-16 -8.94427191e-01]
## [ 7.32343785e-17 -1.00000000e+00  2.56459883e-16]
## [-8.94427191e-01  5.55111512e-17  4.47213595e-01]]
##[[ 1.  2.  2.]
## [ 1. -2.  2.]]
##[[ True  True  True]
## [ True  True  True]]
##U:
##[[-0.70710678 -0.70710678]
## [-0.70710678  0.70710678]]
##[[1.00000000e+00 2.32191141e-16]
## [2.32191141e-16 1.00000000e+00]]
##V:
##[[-4.47213595e-01  7.32343785e-17 -8.94427191e-01]
## [-4.15130595e-16 -1.00000000e+00  5.55111512e-17]
## [-8.94427191e-01  2.56459883e-16  4.47213595e-01]]
##[[ 1.00000000e+00  1.52994492e-16 -1.23617226e-16]
## [ 1.52994492e-16  1.00000000e+00 -6.32162446e-18]
## [-1.23617226e-16 -6.32162446e-18  1.00000000e+00]]

1.7 正定矩阵和非负定矩阵

a=np.array([[5,2],[2,1]])
print(a)
d = np.linalg.eig(a)
print(d)
##[[5 2]
## [2 1]]
##(array([5.82842712, 0.17157288]), array([[ 0.92387953, -0.38268343],
##       [ 0.38268343,  0.92387953]]))

1.8 特征值的极值问题

例 1-1.1

a = np.array([[1,2,3,4,5],[2,4,7,8,9],[3,7,10,15,20],[4,8,15,30,20],[5,9,20,20,40]])
print(a)
##[[ 1  2  3  4  5]
## [ 2  4  7  8  9]
## [ 3  7 10 15 20]
## [ 4  8 15 30 20]
## [ 5  9 20 20 40]]

a = np.array([[1,2,3,4,5],[2,4,7,8,9],[3,7,10,15,20],[4,8,15,30,20],[5,9,20,20,40]])
print(np.linalg.inv(a))
##[[ 9.78873239e+00 -2.18309859e+00 -1.85915493e+00  1.12676056e-01
##   1.40845070e-01]
## [-2.18309859e+00  7.74647887e-01  7.88732394e-01 -1.69014085e-01
##  -2.11267606e-01]
## [-1.85915493e+00  7.88732394e-01  3.94366197e-02 -8.45070423e-03
##   3.94366197e-02]
## [ 1.12676056e-01 -1.69014085e-01 -8.45070423e-03  7.32394366e-02
##  -8.45070423e-03]
## [ 1.40845070e-01 -2.11267606e-01  3.94366197e-02 -8.45070423e-03
##   3.94366197e-02]]


b = np.linalg.eig(a)
print(b)
##(array([70.33488803, 14.44024095,  1.997606  ,  0.09374538, -1.86648037]), ##array([[ 0.10513926,  0.00733125,  0.26673691,  0.95627367, -0.05730686],
##       [ 0.20596656,  0.05549834,  0.82858975, -0.22629386,  0.46554035],
##       [ 0.39707684, -0.02585507,  0.32382661, -0.18402887, -0.83840992],
##       [ 0.5462168 ,  0.78569385, -0.25851683,  0.01391472,  0.13155919],
##       [ 0.70035756, -0.61553463, -0.26569873,  0.0164776 ,  0.24443626]]))

#行列式：
a = np.array([[1,2,3,4,5],[2,4,7,8,9],[3,7,10,15,20],[4,8,15,30,20],[5,9,20,20,40]])
result1 = np.linalg.det(a)
print(result1)
##-355.00000

#矩阵的迹：
b = a.trace()
print(b)
##85