数据挖掘_实验二 数据预处理常用方法与应用
实验二 数据预处理常用方法与应用
【实验目的】
1 掌握随机抽样的常用方法。
2 掌握变量规范化的常用方法。
3 掌握常见距离的定义及计算方法。
4 掌握二元变量相似度、连续变量余弦相似度计算方法。
5 掌握Mahalanobis距离的计算方法。
【实验类型】
设计型
【实验学时】
2学时
【实验环境】
Windows 7以上操作系统
Python3.0以上版本
Pycharm开发环境
Spyder开发环境
【实验要求】
一、请随机生成1000个二维数据点，从中分别随机抽取约75%，50%，25%的数据点，并绘图展示抽取数据点的情况。
参考输出结果:
code:
import matplotlib.pyplot as plt
import random as rd
x = [rd.gauss(0, 1) for i in range(1000)]
y = [rd.gauss(0, 1) for i in range(1000)]
def f(x, y, p):
xx = []
yy = []
for j in range(1000):
if rd.random() >= p:
xx.append(x[j])
yy.append(y[j])
return xx, yy
def paint(xx, yy, str_1, str_2):
plt.subplot(str_1)
plt.title(str_2+str(len(xx)))
plt.scatter(xx, yy, s=[0.5])
x1, y1 = f(x, y, 0.25)
x2, y2 = f(x, y, 0.5)
x3, y3 = f(x, y, 0.75)
paint(x, y, 221, "all--")
paint(x1, y1, 222, "75%--")
paint(x2, y2, 223, "50%--")
paint(x3, y3, 224, "25%--")
plt.show()
二、请计算生成正弦函数序列和余弦函数序列，并绘图展示。将正弦序列和余弦序列添加随机偏差（噪音）后再次绘图，观察噪音对数据的影响。
参考输出结果:
code:
from math import *
import matplotlib.pyplot as plt
import numpy as np
x = np.arange(-6*pi, 6*pi, 0.1)
y = np.sin(x)
z = np.cos(x)
y1 = y+np.random.randn(len(y))
z1 = y+np.random.randn(len(z))
def paint(xx, yy, zz, str_1):
plt.subplot(str_1)
plt.plot(xx, yy)
plt.plot(xx, zz)
paint(x, y, z, 211)
paint(x, y1, z1, 212)
plt.show()
三、验证几类变量变换的效果，包括：指数变换、> 1幂变换、对数变换、（0，1）幂变换、< 0幂变换、正弦变换，绘图演示各类变换的效果。
参考输出结果：
code:
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(1, 50, 100)
y1 = np.exp(x)
y2 = x**3
y3 = np.log10(x)
y4 = np.sqrt(x)
y5 = 1/x
y6 = np.sin(x)
def paint(x, y, str_1, str_2):
plt.subplot(str_2)
plt.title(r"{}".format(str_1))
plt.plot(x, y)
paint(x, y1, "$y=e^x$", 231)
paint(x, y2, "$y=x^3$", 232)
paint(x, y3, "$y=log_{10}x$", 233)
paint(x, y4, "$y=\sqrt{x}$", 234)
paint(x, y5, '$y=\\frac{1}{x}$', 235)
paint(x, y6, "$y=sin(x)$", 236)
plt.show()
"""
输出函数式参考： https://stackoverflow.com/questions/4028267/how-to-render-latex-markup-using-python/58496145#58496145
"""
四、为了提出变量量纲对数据分析的影响，经常需要对变量进行一定的变换，其中常用的一种是标准化（normalization）变换（即将数据序列转换为期望值为0和方差为1的序列），转换公式为：
其中
为序列的均值，
为序列的样本标准差，即
%5E2%7D)
- (这儿借用了谷歌的api需要科学网络才能显示这四个公式)
请随机生成一个序列，并生成它的标准化序列，将两个序列通过绘图展示出来。
参考输出结果（其中随机生成的序列为期望值1000，标准差200，长度1000的正态分布序列）：
code:
from math import *
import matplotlib.pyplot as plt
import numpy as np
y1 = np.random.normal(1000, 200, 1000)
sd = sqrt(sum((y1-np.average(y1))**2))/(len(y1)-1)
y2 = (y1-np.average(y1))/sd
print(np.mean(y2), np.std(y2))
plt.subplot(211)
plt.plot(y1)
plt.subplot(212)
plt.plot(y2)
plt.show()
五、距离计算
向量间的闵可夫斯基距离的定义如下：
其中当r=1时成为街区距离，为不同分量差绝对值的和
当r=2时即为欧几里得距离
当r=+infinity时，等于所有分量差绝对值中最大的一个
编写函数实现三种距离的计算，并对下面4个二维向量分别计算其三种距离定义的距离矩阵。
参考输出结果如下：
code:
import math
import re
infinity = float("inf")
def d(point1, point2, r):
global infinity
sum = 0.0
if r != infinity:
for i in range(len(point1)):
sum += (abs(point1[i]-point2[i])**r)
sum **= (1/r)
else:
for i in range(len(point1)):
xyabs_ = abs(point1[i]-point2[i])
sum = sum if sum > xyabs_ else xyabs_
return sum
def ma(list_1, r):
list_2 = []
for i in list_1:
list_3 = []
for j in list_1:
list_3.append(d(i, j, r))
list_2.append(list_3)
return list_2
def out(list_):
for i in list_:
for j in i:
print("{}\t".format(j), end="")
print()
print()
list1 = [[0, 2], [2, 0], [3, 1], [5, 1]]
print("Data Matrix:")
out(list1)
print("City Block Distance:")
out(ma(list1, 1))
print("Euclid Distance:")
out(ma(list1, 2))
print("Distance with Element Max Difference:")
out(ma(list1, infinity))
"""
题目中的距离公式不是太容易理解，下面的文档有一些不错的解释
闵可夫斯基距离详解：https://wenku.baidu.com/view/35296aa3ec630b1c59eef8c75fbfc77da3699700.html
"""
"""
尝试了一下参考图中的那种输出，没有成功；
此代码的输出：
Data Matrix:
0 2
2 0
3 1
5 1
City Block Distance:
0.0 4.0 4.0 6.0
4.0 0.0 2.0 4.0
4.0 2.0 0.0 2.0
6.0 4.0 2.0 0.0
Euclid Distance:
0.0 2.8284271247461903 3.1622776601683795 5.0990195135927845
2.8284271247461903 0.0 1.4142135623730951 3.1622776601683795
3.1622776601683795 1.4142135623730951 0.0 2.0
5.0990195135927845 3.1622776601683795 2.0 0.0
Distance with Element Max Difference:
0 2 3 5
2 0 1 3
3 1 0 2
5 3 2 0
"""
六、二元变量的相似度衡量可以使用简单匹配系数（SMC）或Jaccard系数
设x和y是两个对象，都由n个二元属性组成，这样两个对象可生成如下四个量（频率）：
简单匹配系数（Simple Matching Coefficient, SMC）定义为：
Jaccard系数（Jaccard Coefficient）的定义为：
编写两个函数分别实现连个系数的计算，并生成两个随机序列调用函数计算。
参考输出结果:
code
import numpy as np
def fun_count(a, b):
count00 = 0
count11 = 0
count10 = 0
count01 = 0
for i in range(0, len(a)):
if a[i] == 0:
if a[i] == b[i]:
count00 += 1
else:
count01 += 1
else:
""" a[i]==1 """
if a[i] == b[i]:
count11 += 1
else:
count10 += 1
return count00, count01, count10, count11
def fun_SCM(a, b):
f00, f01, f10, f11 = fun_count(a, b)
SCM = ((f11+f00)*1.0)/(f00+f01+f10+f11)
return SCM
def fun_J(a, b):
f00, f01, f10, f11 = fun_count(a, b)
J = (f11*1.0)/(f01+f10+f11)
return J
a = np.random.randint(0, 2, 100)
b = np.random.randint(0, 2, 100)
print("Vector a is:")
print(a)
print("Vector b is:")
print(b)
print("SCM coefficient:")
print(fun_SCM(a, b))
print("Jaccard coefficient")
print(fun_J(a, b))
七、余弦相似度是衡量数值向量相似度的一种方法，设x和y是两个n元向量，则余弦相似度定义为
编写函数实现余弦相似度的计算，并调用函数计算下面两个向量的余弦相似度。
参考输出结果：
code
def fun_xy(a, b):
sum = 0.0
for i in range(len(a)):
sum += (a[i]*b[i])
return sum
def fun_len(x):
sum = 0.0
for i in x:
sum += i**2
return sum**(1/2)
def fun_cosxy(x, y):
return (fun_xy(x, y)*1.0)/(fun_len(x)*fun_len(y))
x = (3, 2, 0, 5, 0, 0, 0, 2, 0, 0)
y = (1, 0, 0, 0, 0, 0, 0, 1, 0, 2)
print("Vetor Angle Consine:{}".format(fun_cosxy(x, y)))
八、从统计学的观点看，计算多维向量的距离时，应该使用mahalanobis距离，其定义为：
属性间协方差的定义如下所示：
实现欧几里得距离、使用方差逆矩阵调整的距离、使用协方差逆矩阵调整的距离三类距离的计算函数，并构造两个随机向量验证三种距离的不同。
参考输出结果：
code
from math import *
import re
import numpy as np
# 计算两个矩阵的乘积
# numpy的对应函数 np.linalg.dot(X,Y)
def i_matrix_mul(X, Y):
rrow = np.size(X, axis=0)
rcol = np.size(Y, axis=1)
rmid = np.size(X, axis=1)
rel = np.linspace(0.0, 0.0, rrow*rcol).reshape(rrow, rcol)
for r in range(rrow):
for c in range(rcol):
rel[r, c] = 0.0
for mid in range(rmid):
rel[r, c] += X[r, mid] * Y[mid, c]
return rel
# 计算期望值（均值），对应np.average
def i_mean(x):
return sum(x)/len(x)
# 计算样本方差
# 注意，numpy.var计算的是总体的，使用时需要调整
def i_var(x):
x_mean = i_mean(x)
return sum((x-x_mean)**2)/(len(x)-1)
# 计算样本标准差
# 注意，numpy.std计算的是总体的，使用时需要调整
def i_std(x):
x_mean = i_mean(x)
return sqrt(i_var(x))
# 序列标准化
def i_norm_vec(x):
x_mean = i_mean(x)
x_std = i_std(x)
return (x-x_mean)/x_std
# 计算两个向量的协方差
# numpy的对应方法 np.cov(x,y),但是返回协方差矩阵
def i_cov(x, y):
x_av = np.average(x)
y_av = np.average(y)
return sum((x-x_av)*(y-y_av))/(len(x)-1)
# 计算一个numpy矩阵的协方差（每行是一个变量）
# numpy的对应方法 np.cov(X)
def i_cov_matrix(X):
vc = np.size(X, axis=0)
cm = np.linspace(0.0, 0.0, vc**2).reshape(vc, vc)
for i in range(vc):
for j in range(i, vc):
cm[j, i] = cm[i, j] = i_cov(X[i, ...], X[j, ...])
return cm
# 计算方差矩阵，（一般不用这个函数，这里为实验演示添加）
def i_var_matrix(X):
vc = np.size(X, axis=0)
cm = np.linspace(0.0, 0.0, vc**2).reshape(vc, vc)
for i in range(vc):
cm[i, i] = i_cov(X[i, ...], X[i, ...])
return cm
# 计算两个向量的相关系数（先标准化再计算协方差）
# numpy的对应方法 np.corrcoef(x,y)
def i_cor(x, y):
x_s = i_norm_vec(x)
y_s = i_norm_vec(y)
return i_cov(x_s, y_s)
# 计算一个numpy矩阵的相关系数阵（每行是一个变量）
# numpy的对应方法 np.corrcoef(X)
def i_cor_matrix(X):
vc = np.size(X, axis=0)
cm = np.linspace(0.0, 0.0, vc**2).reshape(vc, vc)
for i in range(vc):
for j in range(i, vc):
cm[i, j] = i_cov(X[i, ...], X[i, ...])
return cm
# 计算矩阵的逆，注意，未检查奇异矩阵
# numpy的对应函数np.lialq.inv(Y)
# 其实当不能从下面找到交换行时，就可以判断奇异矩阵了
def i_matrx_inv(Y):
X = np.copy(Y) # 不破换原来的矩阵
n = np.size(X, axis=0)
rel = np.identity(n)
""" 先变为右上三角举证 """
for i in range(n):
""" 如果（i，i)位置为0则从下面行向上交换 """
if abs(X[i, i]) < 1e-5:
for r in range(i+1, n):
if X[r, i] > 1e-5:
X[i, ...], X[r, ...] = X[r, ...], X[i, ...]
rel[i, ...], rel[r, ...] = rel[r, ...], rel[i, ...] # 伴随矩阵要一块变化
break
""" 如果找不到，可以触发异常，是奇异矩阵 """
""" 当前行（i，j）位置变为1 """
t = float(X[i, i])
for j in range(i, n):
""" 只需处理i后面的列就可以，因为前面的已经都变成0了 """
X[j, j] /= t
for j in range(n):
"""伴随矩阵要一块变化且要变化所有的列 """
rel[i, j] /= t
""" 将所有下面行的相应位置变为0 """
for row in range(i+1, n):
rh = float(X[row, i])
for col in range(i, n):
""" 这里不必处理所有的列，因为前面已经都变为0了 """
X[row, col] -= rh*X[i, col]
for col in range(n):
""" 伴随矩阵要一块变化并且要处理所有的列 """
rel[row, col] -= rh*rel[i, col]
for i in range(n-1, 0, -1):
for row in range(i-1, -1, -1):
rh = float(X[row, i])
for col in range(n):
""" 伴随矩阵变化，要处理所有的列 """
rel[row, col] -= rh*rel[i, col]
return rel
# 计算欧几里得距离
# 也可以用numpy的计算模的函数完成计算np.linalg.norm(x-y)
def dist(x, y):
return sqrt(sum([(x[i]-y[i]) ** 2 for i in range(len(x))]))
# 距离计算的例子,生成两个相关序列，观察其不同的距离
np.random.seed(856479)
a = np.random.randn(100) # 标准正态分布的序列
b = 50*a+10000 # 期望值为10000,方差为2500的正态分布，且与a完全线性相关
m = np.vstack((a, b))
# 首尾两个点的距离
fp = np.array([a[0], b[0]])
lp = np.array([a[99], b[99]])
print("Two data points:")
print(fp)
print(lp)
print()
# 完全未调整的距离
print("Unadjusted diatance :{}".format(dist(fp, lp)))
print()
# 计算相关矩阵
dv = (fp-lp).reshape(1, 2)
dv_t = dv.T
print("Difference vector")
print(dv)
print("Transposed difference vecyor:")
print(dv_t)
print()
# 经过方差调整的距离
print("Adjusted only by variance:")
m_var = i_var_matrix(m)
m_var_inv = i_matrx_inv(m_var)
print("Variance Matrix:")
print(m_var)
print("Inver variance matrix:")
print(m_var_inv)
print("Distance by variance:")
print(np.sqrt(i_matrix_mul(i_matrix_mul(dv, m_var_inv), dv_t)))
print()
# 经协方差矩阵调整的距离
print("Adjusted by covariance：")
m_cov = i_cov_matrix(m)
m_cov_inv = i_matrx_inv(m_cov)
print("Covariance matrix：")
print(m_cov)
print("Inverse covariance matrix:")
print(m_cov_inv)
print("Distance by covariance:")
print(np.sqrt(i_matrix_mul(i_matrix_mul(dv, m_cov_inv), dv_t)))
【实验目的】
1 掌握随机抽样的常用方法。
2 掌握变量规范化的常用方法。
3 掌握常见距离的定义及计算方法。
4 掌握二元变量相似度、连续变量余弦相似度计算方法。
5 掌握Mahalanobis距离的计算方法。
【实验类型】
设计型
【实验学时】
2学时
【实验环境】
Windows 7以上操作系统
Python3.0以上版本
Pycharm开发环境
Spyder开发环境
【实验要求】
一、请随机生成1000个二维数据点，从中分别随机抽取约75%，50%，25%的数据点，并绘图展示抽取数据点的情况。
参考输出结果:
code:
import matplotlib.pyplot as plt
import random as rd
x = [rd.gauss(0, 1) for i in range(1000)]
y = [rd.gauss(0, 1) for i in range(1000)]
def f(x, y, p):
xx = []
yy = []
for j in range(1000):
if rd.random() >= p:
xx.append(x[j])
yy.append(y[j])
return xx, yy
def paint(xx, yy, str_1, str_2):
plt.subplot(str_1)
plt.title(str_2+str(len(xx)))
plt.scatter(xx, yy, s=[0.5])
x1, y1 = f(x, y, 0.25)
x2, y2 = f(x, y, 0.5)
x3, y3 = f(x, y, 0.75)
paint(x, y, 221, "all--")
paint(x1, y1, 222, "75%--")
paint(x2, y2, 223, "50%--")
paint(x3, y3, 224, "25%--")
plt.show()
二、请计算生成正弦函数序列和余弦函数序列，并绘图展示。将正弦序列和余弦序列添加随机偏差（噪音）后再次绘图，观察噪音对数据的影响。
参考输出结果:
code:
from math import *
import matplotlib.pyplot as plt
import numpy as np
x = np.arange(-6*pi, 6*pi, 0.1)
y = np.sin(x)
z = np.cos(x)
y1 = y+np.random.randn(len(y))
z1 = y+np.random.randn(len(z))
def paint(xx, yy, zz, str_1):
plt.subplot(str_1)
plt.plot(xx, yy)
plt.plot(xx, zz)
paint(x, y, z, 211)
paint(x, y1, z1, 212)
plt.show()
三、验证几类变量变换的效果，包括：指数变换、> 1幂变换、对数变换、（0，1）幂变换、< 0幂变换、正弦变换，绘图演示各类变换的效果。
参考输出结果：
code:
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(1, 50, 100)
y1 = np.exp(x)
y2 = x**3
y3 = np.log10(x)
y4 = np.sqrt(x)
y5 = 1/x
y6 = np.sin(x)
def paint(x, y, str_1, str_2):
plt.subplot(str_2)
plt.title(r"{}".format(str_1))
plt.plot(x, y)
paint(x, y1, "$y=e^x$", 231)
paint(x, y2, "$y=x^3$", 232)
paint(x, y3, "$y=log_{10}x$", 233)
paint(x, y4, "$y=\sqrt{x}$", 234)
paint(x, y5, '$y=\\frac{1}{x}$', 235)
paint(x, y6, "$y=sin(x)$", 236)
plt.show()
"""
输出函数式参考： https://stackoverflow.com/questions/4028267/how-to-render-latex-markup-using-python/58496145#58496145
"""
四、为了提出变量量纲对数据分析的影响，经常需要对变量进行一定的变换，其中常用的一种是标准化（normalization）变换（即将数据序列转换为期望值为0和方差为1的序列），转换公式为：
其中
为序列的均值，
为序列的样本标准差，即
- (这儿借用了谷歌的api需要科学网络才能显示这四个公式)
请随机生成一个序列，并生成它的标准化序列，将两个序列通过绘图展示出来。
参考输出结果（其中随机生成的序列为期望值1000，标准差200，长度1000的正态分布序列）：
code:
from math import *
import matplotlib.pyplot as plt
import numpy as np
y1 = np.random.normal(1000, 200, 1000)
sd = sqrt(sum((y1-np.average(y1))**2))/(len(y1)-1)
y2 = (y1-np.average(y1))/sd
print(np.mean(y2), np.std(y2))
plt.subplot(211)
plt.plot(y1)
plt.subplot(212)
plt.plot(y2)
plt.show()
五、距离计算
向量间的闵可夫斯基距离的定义如下：
其中当r=1时成为街区距离，为不同分量差绝对值的和
当r=2时即为欧几里得距离
当r=+infinity时，等于所有分量差绝对值中最大的一个
编写函数实现三种距离的计算，并对下面4个二维向量分别计算其三种距离定义的距离矩阵。
参考输出结果如下：
code:
import math
import re
infinity = float("inf")
def d(point1, point2, r):
global infinity
sum = 0.0
if r != infinity:
for i in range(len(point1)):
sum += (abs(point1[i]-point2[i])**r)
sum **= (1/r)
else:
for i in range(len(point1)):
xyabs_ = abs(point1[i]-point2[i])
sum = sum if sum > xyabs_ else xyabs_
return sum
def ma(list_1, r):
list_2 = []
for i in list_1:
list_3 = []
for j in list_1:
list_3.append(d(i, j, r))
list_2.append(list_3)
return list_2
def out(list_):
for i in list_:
for j in i:
print("{}\t".format(j), end="")
print()
print()
list1 = [[0, 2], [2, 0], [3, 1], [5, 1]]
print("Data Matrix:")
out(list1)
print("City Block Distance:")
out(ma(list1, 1))
print("Euclid Distance:")
out(ma(list1, 2))
print("Distance with Element Max Difference:")
out(ma(list1, infinity))
"""
题目中的距离公式不是太容易理解，下面的文档有一些不错的解释
闵可夫斯基距离详解：https://wenku.baidu.com/view/35296aa3ec630b1c59eef8c75fbfc77da3699700.html
"""
"""
尝试了一下参考图中的那种输出，没有成功；
此代码的输出：
Data Matrix:
0 2
2 0
3 1
5 1
City Block Distance:
0.0 4.0 4.0 6.0
4.0 0.0 2.0 4.0
4.0 2.0 0.0 2.0
6.0 4.0 2.0 0.0
Euclid Distance:
0.0 2.8284271247461903 3.1622776601683795 5.0990195135927845
2.8284271247461903 0.0 1.4142135623730951 3.1622776601683795
3.1622776601683795 1.4142135623730951 0.0 2.0
5.0990195135927845 3.1622776601683795 2.0 0.0
Distance with Element Max Difference:
0 2 3 5
2 0 1 3
3 1 0 2
5 3 2 0
"""
六、二元变量的相似度衡量可以使用简单匹配系数（SMC）或Jaccard系数
设x和y是两个对象，都由n个二元属性组成，这样两个对象可生成如下四个量（频率）：
简单匹配系数（Simple Matching Coefficient, SMC）定义为：
Jaccard系数（Jaccard Coefficient）的定义为：
编写两个函数分别实现连个系数的计算，并生成两个随机序列调用函数计算。
参考输出结果:
code
import numpy as np
def fun_count(a, b):
count00 = 0
count11 = 0
count10 = 0
count01 = 0
for i in range(0, len(a)):
if a[i] == 0:
if a[i] == b[i]:
count00 += 1
else:
count01 += 1
else:
""" a[i]==1 """
if a[i] == b[i]:
count11 += 1
else:
count10 += 1
return count00, count01, count10, count11
def fun_SCM(a, b):
f00, f01, f10, f11 = fun_count(a, b)
SCM = ((f11+f00)*1.0)/(f00+f01+f10+f11)
return SCM
def fun_J(a, b):
f00, f01, f10, f11 = fun_count(a, b)
J = (f11*1.0)/(f01+f10+f11)
return J
a = np.random.randint(0, 2, 100)
b = np.random.randint(0, 2, 100)
print("Vector a is:")
print(a)
print("Vector b is:")
print(b)
print("SCM coefficient:")
print(fun_SCM(a, b))
print("Jaccard coefficient")
print(fun_J(a, b))
七、余弦相似度是衡量数值向量相似度的一种方法，设x和y是两个n元向量，则余弦相似度定义为
编写函数实现余弦相似度的计算，并调用函数计算下面两个向量的余弦相似度。
参考输出结果：
code
def fun_xy(a, b):
sum = 0.0
for i in range(len(a)):
sum += (a[i]*b[i])
return sum
def fun_len(x):
sum = 0.0
for i in x:
sum += i**2
return sum**(1/2)
def fun_cosxy(x, y):
return (fun_xy(x, y)*1.0)/(fun_len(x)*fun_len(y))
x = (3, 2, 0, 5, 0, 0, 0, 2, 0, 0)
y = (1, 0, 0, 0, 0, 0, 0, 1, 0, 2)
print("Vetor Angle Consine:{}".format(fun_cosxy(x, y)))
八、从统计学的观点看，计算多维向量的距离时，应该使用mahalanobis距离，其定义为：
属性间协方差的定义如下所示：
实现欧几里得距离、使用方差逆矩阵调整的距离、使用协方差逆矩阵调整的距离三类距离的计算函数，并构造两个随机向量验证三种距离的不同。
参考输出结果：
code
from math import *
import re
import numpy as np
# 计算两个矩阵的乘积
# numpy的对应函数 np.linalg.dot(X,Y)
def i_matrix_mul(X, Y):
rrow = np.size(X, axis=0)
rcol = np.size(Y, axis=1)
rmid = np.size(X, axis=1)
rel = np.linspace(0.0, 0.0, rrow*rcol).reshape(rrow, rcol)
for r in range(rrow):
for c in range(rcol):
rel[r, c] = 0.0
for mid in range(rmid):
rel[r, c] += X[r, mid] * Y[mid, c]
return rel
# 计算期望值（均值），对应np.average
def i_mean(x):
return sum(x)/len(x)
# 计算样本方差
# 注意，numpy.var计算的是总体的，使用时需要调整
def i_var(x):
x_mean = i_mean(x)
return sum((x-x_mean)**2)/(len(x)-1)
# 计算样本标准差
# 注意，numpy.std计算的是总体的，使用时需要调整
def i_std(x):
x_mean = i_mean(x)
return sqrt(i_var(x))
# 序列标准化
def i_norm_vec(x):
x_mean = i_mean(x)
x_std = i_std(x)
return (x-x_mean)/x_std
# 计算两个向量的协方差
# numpy的对应方法 np.cov(x,y),但是返回协方差矩阵
def i_cov(x, y):
x_av = np.average(x)
y_av = np.average(y)
return sum((x-x_av)*(y-y_av))/(len(x)-1)
# 计算一个numpy矩阵的协方差（每行是一个变量）
# numpy的对应方法 np.cov(X)
def i_cov_matrix(X):
vc = np.size(X, axis=0)
cm = np.linspace(0.0, 0.0, vc**2).reshape(vc, vc)
for i in range(vc):
for j in range(i, vc):
cm[j, i] = cm[i, j] = i_cov(X[i, ...], X[j, ...])
return cm
# 计算方差矩阵，（一般不用这个函数，这里为实验演示添加）
def i_var_matrix(X):
vc = np.size(X, axis=0)
cm = np.linspace(0.0, 0.0, vc**2).reshape(vc, vc)
for i in range(vc):
cm[i, i] = i_cov(X[i, ...], X[i, ...])
return cm
# 计算两个向量的相关系数（先标准化再计算协方差）
# numpy的对应方法 np.corrcoef(x,y)
def i_cor(x, y):
x_s = i_norm_vec(x)
y_s = i_norm_vec(y)
return i_cov(x_s, y_s)
# 计算一个numpy矩阵的相关系数阵（每行是一个变量）
# numpy的对应方法 np.corrcoef(X)
def i_cor_matrix(X):
vc = np.size(X, axis=0)
cm = np.linspace(0.0, 0.0, vc**2).reshape(vc, vc)
for i in range(vc):
for j in range(i, vc):
cm[i, j] = i_cov(X[i, ...], X[i, ...])
return cm
# 计算矩阵的逆，注意，未检查奇异矩阵
# numpy的对应函数np.lialq.inv(Y)
# 其实当不能从下面找到交换行时，就可以判断奇异矩阵了
def i_matrx_inv(Y):
X = np.copy(Y) # 不破换原来的矩阵
n = np.size(X, axis=0)
rel = np.identity(n)
""" 先变为右上三角举证 """
for i in range(n):
""" 如果（i，i)位置为0则从下面行向上交换 """
if abs(X[i, i]) < 1e-5:
for r in range(i+1, n):
if X[r, i] > 1e-5:
X[i, ...], X[r, ...] = X[r, ...], X[i, ...]
rel[i, ...], rel[r, ...] = rel[r, ...], rel[i, ...] # 伴随矩阵要一块变化
break
""" 如果找不到，可以触发异常，是奇异矩阵 """
""" 当前行（i，j）位置变为1 """
t = float(X[i, i])
for j in range(i, n):
""" 只需处理i后面的列就可以，因为前面的已经都变成0了 """
X[j, j] /= t
for j in range(n):
"""伴随矩阵要一块变化且要变化所有的列 """
rel[i, j] /= t
""" 将所有下面行的相应位置变为0 """
for row in range(i+1, n):
rh = float(X[row, i])
for col in range(i, n):
""" 这里不必处理所有的列，因为前面已经都变为0了 """
X[row, col] -= rh*X[i, col]
for col in range(n):
""" 伴随矩阵要一块变化并且要处理所有的列 """
rel[row, col] -= rh*rel[i, col]
for i in range(n-1, 0, -1):
for row in range(i-1, -1, -1):
rh = float(X[row, i])
for col in range(n):
""" 伴随矩阵变化，要处理所有的列 """
rel[row, col] -= rh*rel[i, col]
return rel
# 计算欧几里得距离
# 也可以用numpy的计算模的函数完成计算np.linalg.norm(x-y)
def dist(x, y):
return sqrt(sum([(x[i]-y[i]) ** 2 for i in range(len(x))]))
# 距离计算的例子,生成两个相关序列，观察其不同的距离
np.random.seed(856479)
a = np.random.randn(100) # 标准正态分布的序列
b = 50*a+10000 # 期望值为10000,方差为2500的正态分布，且与a完全线性相关
m = np.vstack((a, b))
# 首尾两个点的距离
fp = np.array([a[0], b[0]])
lp = np.array([a[99], b[99]])
print("Two data points:")
print(fp)
print(lp)
print()
# 完全未调整的距离
print("Unadjusted diatance :{}".format(dist(fp, lp)))
print()
# 计算相关矩阵
dv = (fp-lp).reshape(1, 2)
dv_t = dv.T
print("Difference vector")
print(dv)
print("Transposed difference vecyor:")
print(dv_t)
print()
# 经过方差调整的距离
print("Adjusted only by variance:")
m_var = i_var_matrix(m)
m_var_inv = i_matrx_inv(m_var)
print("Variance Matrix:")
print(m_var)
print("Inver variance matrix:")
print(m_var_inv)
print("Distance by variance:")
print(np.sqrt(i_matrix_mul(i_matrix_mul(dv, m_var_inv), dv_t)))
print()
# 经协方差矩阵调整的距离
print("Adjusted by covariance：")
m_cov = i_cov_matrix(m)
m_cov_inv = i_matrx_inv(m_cov)
print("Covariance matrix：")
print(m_cov)
print("Inverse covariance matrix:")
print(m_cov_inv)
print("Distance by covariance:")
print(np.sqrt(i_matrix_mul(i_matrix_mul(dv, m_cov_inv), dv_t)))