应用多元统计分析（2）第一章矩阵代数

最新推荐文章于 2022-10-17 18:55:05 发布

HarryStudyPython_ing

最新推荐文章于 2022-10-17 18:55:05 发布

阅读量737

点赞数 2

分类专栏： Python 应用多元统计分析（教材：王学明）文章标签： python 深度学习多元统计矩阵代数 numpy

本文链接：https://blog.csdn.net/qq_47805483/article/details/108621082

版权

Python 应用多元统计分析（教材：王学明）专栏收录该内容

9 篇文章 12 订阅

订阅专栏

课堂代码及注释

基础矩阵运算代码

###lecture1###
###第一章 矩阵代数

import numpy as np

#生成两个矩阵(列向量)
a=np.array([1,3,6,9])
b=np.array([20,35,56,28])

#另外的方法
np.arange(10)
np.linspace(1,10,5)

#简单矩阵
np.zeros(4)
np.zeros((4,4))
np.ones((3,4))
np.eye(4)

np.arange(9)
B=np.arange(16).reshape(4,4);B

np.triu(B)##上三角行列式
np.tril(B)##下三角行列式

#求a的长度
np.linalg.norm(a)

#矩阵乘法(a的转置乘以b)
a.T.dot(b)
np.dot(a.T, b)

#生成4*3矩阵
A=np.array([[1,3,6],[2,7,9],[4,8,5],[10,15,7]])##注意有两个[[]]

A.T.dot(A)
3*A

#计算行列式
np.linalg.det(B)

#求逆矩阵
np.linalg.inv(B)

#随机生成4*4矩阵
np.random.seed(10)##设置随机种子
B=np.random.randint(1,20,size=(4,4));B##生成随机整数（范围+size）
np.linalg.det(B)#计算行列式
np.linalg.inv(B)#求逆矩阵

np.random.rand(4,4)##随机生成矩阵（范围0到1）

###验证性质    
# 若  A : p × q ；B : q × p ，则|I p + AB| = |I q + BA|
np.random.seed(100)##设置随机种子
A=np.random.randint(1,20,size=(3,5));A## alt+shift+向下箭头 ：自动向下复制一行
B=np.random.randint(1,20,size=(5,3));B
# A : p × q ；B : q × p ，则|I p + AB| = |I q + BA|
np.linalg.det(np.eye(3)+A.dot(B))
np.linalg.det(np.eye(5)+B.dot(A))


###验证性质  
# 若 A 和 C 均为 p 阶非退化方阵，则(AC)' = C'A'
A=np.random.randint(1,20,size=(4,4));A
B=np.random.randint(1,20,size=(4,4));B

np.linalg.inv(A.dot(B))
np.linalg.inv(B).dot(np.linalg.inv(A))

###矩阵的置
np.rank(A)##求矩阵的置  有问题 结果不一样
np.linalg.matrix_rank(A)##求矩阵的置

###验证性质  
#若 A 和 C 为非退化方阵，则 rank(ABC) = rank(B)
# 准备工作 生成ABC矩阵
np.random.seed(100)##设置随机种子
A=np.random.randint(1,20,size=(4,3));A
B=np.random.randint(1,20,size=(4,4));B
np.linalg.det(B)
C=np.random.randint(1,20,size=(3,3));C
np.linalg.det(C)
#若 A 和 C 为非退化方阵，则 rank(ABC) = rank(B)
B.dot(A).dot(C)
np.linalg.matrix_rank(B.dot(A).dot(C))
np.linalg.matrix_rank(A)

### 特征根特征向量
np.random.seed(100)##设置随机种子
B=np.random.randint(1,20,size=(4,4));B
values,vector=np.linalg.eig(B);values,vector##计算特征根特征向量
values##分开赋值
vector

vector[:,1]##切片[行,列]
###切片需总结学习
sum(vector[:,0]**2)##验证是否为单位向量

B.dot(vector[:,0])
values[0]*vector[:,0]

##快速得到正交矩阵
np.random.seed(100)##设置随机种子
B=np.random.randint(1,10,size=(3,5));B##生成随机整数（范围+size）
A=B.dot(B.T);A##得到实对称矩
values,vector=np.linalg.eig(A);values,vector
vector[:,0].dot(vector[:,1])##验证是否正交

np.linalg.det(A)
np.prod(np.linalg.eig(A)[0])

np.trace(A)##计算矩阵的迹

np.abs(np.prod(np.linalg.eig(A)[0]))

### 矩阵的迹

##验证tr(AA')=tr(A'A)=所有元素的平方和
np.trace(A)
np.sum(B**2)
## 投影矩阵：正交的逆（幂）等阵


#7正定矩阵和非正定矩阵
np.random.seed(100)##设置随机种子
B=np.random.randint(1,10,size=(3,5));B
A=B.dot(B.T);A##得到实对称矩
values,vector=np.linalg.eig(A);values,vector
values,vector=np.linalg.eig(np.linalg.inv(A));values,vector

np.linalg.matrix_rank(B)

特征值的极值问题

在这里插入图片描述

###验证极值定理
## 设 A 是 p 阶对称矩阵，B 是 p 阶正定矩阵，µ1 ≥ µ2 ≥ · · · ≥ µp 是 Bö 1A 的 p 个特征值，相应的一组特征向量是 t1,t2, · · · ,tp ，则
np.random.seed(100)
L=np.random.randint(1,10,size=(3,5))
A=L.dot(L.T);A

M=np.random.randint(1,6,size=(3,5))
B=M.dot(M.T)

np.linalg.eig(A)
np.linalg.eig(B)

np.linalg.eig(np.linalg.inv(B).dot(A))

矩阵分解

lu分解/三角分解

在这里插入图片描述

##lu分解/三角分解
np.random.seed(100)
A=np.random.randint(1,6,size=(5,5));A

from scipy import *##此处错误
import scipy.linalg
from scipy.linalg import lu
lu=scipy.linalg.lu(A)
type(lu)
P=lu[0];P
L=lu[1];L
U=lu[2];U

P,L,U=scipy.linalg.lu(A);P,L,U
P.dot(L).dot(U)

谱分解（特征值分解）

在这里插入图片描述

##谱分解（特征值分解）
import numpy as np
np.random.seed(0)
A=np.random.randint(1,10,size=(4,4));A
np.linalg.det(A)

np.linalg.eig(A)##特征根分解

lamb,V=np.linalg.eig(A)
D=np.diag(lamb)

V.dot(D).dot(np.linalg.inv(V))

#补充diag()
np.diag(A)
np.diag(np.diag(A))

B=A.dot(A.T);B##生成实对称矩阵
lamb2,V2=np.linalg.eig(B)
D2=np.diag(lamb2);D2##构建对角阵diag
V2.dot(D2).dot(np.linalg.inv(V2))

V2
np.linalg.inv(V2)

奇异值分解

在这里插入图片描述

##奇异值分解
np.random.seed(10)
A=np.random.randint(1,10,size=(6,4));A

np.linalg.svd(A)
U3,d,V_T=np.linalg.svd(A)##!!!
U3,d,V_T=np.linalg.svd(A,full_matrices=False)##!!!
##对比此两行代码
D=np.diag(d)

U3.dot(np.diag(d)).dot(V_T)

##辅助full_matrices=False
M=np.zeros((6,4));M
M[:4,:4]=np.diag(d);M

Choleskey分解

在这里插入图片描述

###Choleskey分解/对于正定矩阵 A=U'U  U为上三角矩阵
np.random.seed(10)
M=np.random.randint(1,5,size=(4,5));M
B=M.dot(M.T)##构造实对称矩阵
np.linalg.eig(B)##判断是否为正定矩阵

U=np.linalg.cholesky(B)
U.dot(U.T)

QR分解

在这里插入图片描述

###QR分解/  （M*N） A=QR 正定矩阵*上三角矩阵
np.random.seed(10)
A=np.random.randint(1,10,size=(6,4));A
Q,R=np.linalg.qr(A,'complete')
np.allclose()
Q.dot(Q.T)

np.allclose(Q.dot(Q.T),np.eye(4))##判断是否为四节单位阵

补充：矩阵计算 Python 笔记

创建一个数组

在 python 中创建数组可以使用 numpy 包中的 array 函数，例如

import numpy as np
x = np.array([1, 2, 3, 4])
x

创建一个矩阵

在 python 中可以用函数 np.matrix() 函数来创建矩阵。在 numpy 中，matrix是 array 的一个小分支，包含于 array，所以 matrix 拥有 array 的所有特性。创建矩阵的方法有很多种，下面举例分析：

np.matrix(x) #返回x的copy
np.mat(x) #不返回copy,等价于np.asmatrix(x)

使用 reshape 函数可以控制矩阵的行数和列数，默认按行排列，参数 order 设置为’F’ 则可以按列排列。例如：

x = np.arange(12)
np.mat(x).reshape((3, 4))
np.mat(x).reshape((3, 4), order='F')

矩阵转置

在这里插入图片描述

A = np.mat(x).reshape((3, 4))
A.T
np.transpose(A)

矩阵相加减

在 python 中对同行同列矩阵相加减，可用符号：“+”、“-”，例如：

A = B = np.mat(range(12)).reshape((3, 4))
A + B
A - B

数与矩阵相乘

A 为 m × n 矩阵,c 为实数，在 Python 中求 cA 可用符号：”*”, 例如：

c = 2
c * A 1

矩阵相乘

A 为 m×n 矩阵,B 为 n×k 矩阵, 在 Python 中求 AB，可以使用函数 np.dot(A,B), 或者直接 A * B，例如：

A = np.mat(range(12)).reshape((3, 4))
A = np.mat(range(12)).reshape((4, 3))
A * B 
np.dot(A, B)

矩阵对角元素相关运算

例如要取一个方阵的对角元素

A = np.mat(range(16)).reshape((4, 4))
np.diag(A)
A.diagonal()

对一个列表或一维数组应用 np.diag() 函数，将产生以该列表或数组为对角元素的数组，例如：

np.diag(np.diag(A))
np.diag([1, 2, 3])

np.eye() 和 np.identity() 函数，输入一个正整数 n，可以产生 n × n 的对角矩阵。例如：

np.eye(3)
np.identity(3)

矩阵求逆

在这里插入图片描述

x = np.random.randn(4, 4)
a = np.mat(x)
np.linalg.inv(a)
np.linalg.solve(a, np.eye(4))
a.I
a.I * a

矩阵的特征值和特征向量

矩阵 A 的谱分解为 A = UΛUT，其中 ΛU 是由 A 的特征值组成的对角矩阵，U 的列为 A 的特征值对于的特征向量，在 Python 中可以使用 np.linalg. eig() 得到 U 和 ΛU。例如：

A = np.eye(4) + 1
values, vectors = np.linalg.eig(A)
vectors.dot(np.diag(values)).dot(vectors.T)
vectors.T.dot(vectors)

矩阵的 Choleskey 分解

对于正定矩阵 A，可对其进行 Cholesky 分解，即 A = PTP，其中 P 为下三角矩阵它要求矩阵的所有特征值必须大于零，故分解的下三角的对角元也是大于零的，在 python 中可以用函数 np.linalg.cholesky() 进行 Cholesky 分解，例如：

P = np.linalg.cholesky(A)
P.dot(P.T)

若矩阵为对称正定矩阵，可以利用 Cholesky 分解求行列式的值，如：

np.prod(P.diagonal() ** 2)
np.linalg.det(A)

矩阵奇异值分解

在这里插入图片描述

A = np.mat(range(18)).reshape((3, 6))
u, s, v_T = np.linalg.svd(A, full_matrices=False)
u * np.diag(s) * v_T 
u * u.T
v_T * v_T.T

矩阵 QR 分解

在这里插入图片描述

A = np.mat(range(16)).reshape((4, 4))
q, r = np.linalg.qr(A)
q * r 
q * q.T

矩阵广义逆 (Moore-Penrose)

在这里插入图片描述

A_ = np.linalg.pinv(A)
#验证性质1 
A * A_ * A 
#验证性质2 
A_ * A * A_ 
#验证性质3 
(A * A_).H
A * A_ 
#验证性质4 
(A_ * A).H
A_ * A

矩阵 Kronecker 积

在这里插入图片描述

A = np.mat(range(4)).reshape((2, 2))
B = np.mat([1]*4).reshape((2, 2))
np.kron(A, B)

矩阵的形状

与数组一样，使用 shape 属性可以得到一个矩阵的形状，例如：

A = np.mat(range(12)).reshape((3, 4))
A.shape
A.shape[0] #返回行数 
A.shape[1] #返回列数

矩阵的行和、列和、行平均、列平均

np.sum(A, axis=1) #行和 
np.sum(A, axis=0) #列和 
np.mean(A, axis=1) #行平均 
np.mean(A, axis=0) #列平均

取矩阵的上、下三角部分

在这里插入图片描述

A = np.mat(range(1, 17)).reshape((4, 4))
np.triu(A)
np.tril(A)
np.triu(A, k=1)
np.triu(A, k=-1)

行列式的值

x = np.mat(np.random.rand(4, 4))
np.linalg.det(x)

向量化算子

在这里插入图片描述

def vec(x):
	return np.mat(x).ravel().T
def vech(x):
	x = np.mat(x)
	t = x[np.tril(x) == x]
	return t 
x = np.arange(12).reshape((3, 4))
vec(x)
vech(x)

时间序列的滞后值

def tslag(x, k=0, trim=False):
	"""
	在时间序列分析中，我们常常要用到一个序列的滞后序列,x为一个数组，k为正整 数，k+1代表滞	后阶数， 若trim为假，则返回序列与原序列长度相同，但含有NA值； 若trim项为真，则返回序列中不含有NA值。
	"""
	c = np.array([])
	x = x.astype(float) 
	for i in range(k+1):
		y = x.copy()
		y[(i+1):] = x[:-(i+1)]
		y[:(i+1)] = np.nan
		c = np.append(c, y)
	c = c.reshape((x.shape[0], (k+1)), order='F')
	if trim == True:
		c = c[~np.isnan(c).any(axis=1)]
	return c

x = np.arange(20)
tslag(x, k=3)
tslag(x, k=3, trim=True)

第二次作业

##设置随机种子为学号最后两位，随机生成 5*5 矩阵 A ，编写 python 程序：
import numpy as np
np.random.seed(47)
A=np.random.randint(0,10,size=(5,5))
A


##      1.求矩阵 A 的行列式、逆、秩、特征根和迹;
np.linalg.det(A)##行列式
np.linalg.inv(A)##逆
np.linalg.matrix_rank(A)##秩
values , vectors = np.linalg.eig(A)
values ##特征根
vectors ##特征向量
np.trace(A) ##迹


##      2.对矩阵 A 进行谱分解并验证;
np.linalg.det(A)
np.linalg.eig(A)##特征根分解

lamb,V=np.linalg.eig(A)
D=np.diag(lamb)

V.dot(D).dot(np.linalg.inv(V)) == A ##FALSE原因 精度问题


##      3.利用 A 矩阵生成一个对称正定阵，进行 Choleskey 分解并验证
M=A.dot(A.T) ; M##构造实对称矩阵
np.linalg.eig(M)##判断是否为正定矩阵

U=np.linalg.cholesky(M) ; U
U.dot(U.T) == M  ##FALSE原因 精度问题
np.array_equal(U.dot(U.T),M) ##FALSE原因 精度问题


##      4.随机生成 4*6 矩阵 B 进行 svd 分解并验证;
np.random.seed(47)
B=np.random.randint(1,10,size=(4,6)) ;  B
np.linalg.svd(B)
U,d,V_T=np.linalg.svd(B,full_matrices=False)
D=np.diag(d)
U.dot(np.diag(d)).dot(V_T)

U.dot(np.diag(d)).dot(V_T) == B ##FALSE原因 精度问题

##      5.验证行列式性质 13
np.random.seed(47)
A1=np.random.randint(1,10,size=(3,4)) ; A1
B1=np.random.randint(10,size=(4,3)) ; B1

np.linalg.det(np.eye(3)+A1.dot(B1)) 
np.linalg.det(np.eye(4)+B1.dot(A1))
np.linalg.det(np.eye(3)+A1.dot(B1)) == np.linalg.det(np.eye(4)+B1.dot(A1)) ##FALSE原因 精度问题

其他笔记

import numpy as np
np.random.seed(202)
A=np.random.randint(0,10,size=(5,5))
A

##求矩阵 A 的行列式、逆、秩、特征根和迹;
np.linalg.det(A)
np.linalg.inv(A)
np.linalg.matrix_rank(A)
values , vectors = np.linalg.eig(A)
vectors
np.trace(A)

##验证公式
np.random.seed(202)
A=np.random.randint(10,size=(3,4))
B=np.random.randint(10,size=(4,3))

np.linalg.det(np.eye(3)+A.dot(B)) == np.linalg.det(np.eye(4)+B.dot(A))

##验证公式(还有疑问)
np.random.seed(202)
A=np.random.randint(10,size=(3,3))
B=np.random.randint(10,size=(4,4))
np.linalg.matrix_rank(A)
np.linalg.matrix_rank(B)
AI=np.linalg.inv(A) ; AI
BI=np.linalg.inv(B) ; BI
Z_1 = np.zeros((3,4)) 
Z_2 = np.zeros((4,3))

C=np.hstack((np.vstack((A,Z_2)),np.vstack((Z_1,B)))) ; C
D=np.hstack((np.vstack((AI,Z_2)),np.vstack((Z_1,BI)))) ; D

####补充 np.stack(),np.hstack(),np.vstack

np.array_equal(np.linalg.inv(C),D)

np.linalg.inv(C) == D

np.linalg.inv(C)
D

##验证公式
np.random.seed(202)
A=np.random.randint(10,size=(3,3))
C=np.random.randint(10,size=(4,4))
B=np.random.randint(10,size=(3,4))
np.linalg.matrix_rank(A)
np.linalg.matrix_rank(C)

np.linalg.matrix_rank(A.dot(B).dot(C)) == np.linalg.matrix_rank(B)

##验证公式
np.random.seed(202)
A=np.random.randint(10,size=(3,4))
B=np.random.randint(10,size=(4,3))
AB=A.dot(B)
BA=B.dot(A)
values1, vectors1 = np.linalg.eig(AB)
values2, vectors2 = np.linalg.eig(BA)
values1
values2
values1 == values2

##验证公式
np.random.seed(202)
A=np.random.randint(10,size=(5,5))
values, vectors = np.linalg.eig(A)

np.linalg.det(A) == np.prod(values)

np.linalg.det(A)
np.prod(values)
##说明：np.prod()函数用来计算所有元素的乘积，对于有多个维度的数组可以指定轴，如axis=1指定计算每一行的乘积。
##补充：np.prod()函数

##验证公式
np.random.seed(202)
A=np.random.randint(10,size=(4,5))

np.trace(A.dot(A.T))
np.trace(A.T.dot(A))
np.sum(A**2)

np.trace(A.dot(A.T)) == np.trace(A.T.dot(A))
np.trace(A.dot(A.T)) == np.sum(A**2)

##验证公式
np.random.seed(202)
A=np.random.randint(10,size=(5,5))
values, vectors = np.linalg.eig(A)

np.trace(A) == np.sum(values)

np.trace(A)
np.sum(values)

##验证公式
A=np.eye(4) * 3 + 1 ; A
B=A**(1/2) ; B 

np.array_equal(A,B**2)

##验证公式(存在疑问)
np.random.seed(202)

C = np.random.randint(10,size=(4,4)) ; C
A = C.dot(C.T) ; A

D = np.random.randint(10,size=(4,5))
B = D.dot(D.T)

np.linalg.eig(B)

L = np.linalg.inv(B).dot(A) ; L
values, vectors = np.linalg.eig(L)
values

t1=vectors[0]
t2=vectors[1]
t3=vectors[2]
t4=vectors[3]

t1.reshape(1,4).dot(A).dot(t1.reshape(4,1))/t1.reshape(1,4).dot(B).dot(t1.reshape(4,1))

t3.reshape(1,4).dot(A).dot(t3.reshape(4,1))/t3.reshape(1,4).dot(B).dot(t3.reshape(4,1))

import pandas as pd
dat=pd.read_csv("C:\python\hw01/trees.csv")
B=np.mat(dat)
A=B.T.dot(B)