1.4 局部加权回归的正规方程推导及实现

import numpy as np
import matplotlib.pyplot as plt

#定义维度
m=40
n=2

#定义一个正态分布，参数分别为均值，方差以及X的行向量
def guassianDistribution(mean,var,x):
    return 1/np.sqrt( 2 * np.pi * var )*np.exp( - (x-mean) ** 2 / (2*var) )

#传入一个m*n维矩阵，一个1*n维的矩阵，默认带宽参数为0.07
def getWeight(X,input_x,bd=0.07):
    W=np.identity(m)
    diff=X-input_x
    for i in range(m):
        W[i,i]=np.exp( diff[i]*diff[i].T/(-2* bd**2) )
    return W

def getTheta(X,Y,W):
    A=X.T * W * X
    B=X.T * W * Y
    # 求秩，需要判断下是否为奇异矩阵
    rank = np.linalg.matrix_rank(A)
    # 不是直接求解就行，是的话求最小二范数二乘解
    if rank == min(m, n):
        theta = np.linalg.solve(A, B)
    else:
        theta = np.linalg.lstsq(A, B)
    return theta

#定义输入矩阵，输出矩阵
X=np.mat([[1,x1] for x1 in np.random.rand(m)])
Y=np.mat(np.array([guassianDistribution(0.5,0.5,x) for x in X[:,1]])).T
input_x=np.mat([1,0.7])

#计算
theta=getTheta(X,Y,getWeight(X,input_x))

#画图
plt.scatter(np.array(X[:,1]),np.array(Y.T))
plt.plot(X[:,1],X*theta)
plt.show()
1.4 局部加权回归的正规方程推导及实现

相关推荐
