python实现简单线性回归

用python实现R的线性模型(lm)中一元线性回归的简单方法，使用R的women示例数据，R的运行结果：

> summary(fit)Call:lm(formula = weight ~ height, data = women)Residuals:    Min      1Q  Median      3Q     Max -1.7333 -1.1333 -0.3833  0.7417  3.1167 Coefficients:             Estimate Std. Error t value Pr(>|t|)    (Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***height        3.45000    0.09114   37.85 1.09e-14 ***---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 1.525 on 13 degrees of freedomMultiple R-squared:  0.991, Adjusted R-squared:  0.9903 F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

python实现的功能包括：
1. 计算pearson相关系数
2. 使用最小二乘法计算回归系数
3. 计算拟合优度判定系数R2
4. 计算估计标准误差Se
5. 计算显著性检验的F和P值

import numpy as npimport scipy.stats as ssclass Lm:    """简单一元线性模型，计算回归系数、拟合优度的判定系数和    估计标准误差，显著性水平"""    def __init__(self, data_source, separator):        self.beta = np.matrix(np.zeros(2))        self.yhat = np.matrix(np.zeros(2))        self.r2 = 0.0        self.se = 0.0        self.f = 0.0        self.msr = 0.0        self.mse = 0.0        self.p = 0.0        data_mat = np.genfromtxt(data_source, delimiter=separator)        self.xarr = data_mat[:, :-1]        self.yarr = data_mat[:, -1]        self.ybar = np.mean(self.yarr)        self.dfd = len(self.yarr) - 2  # 自由度n-2        return    # 计算协方差    @staticmethod    def cov_custom(x, y):        result = sum((x - np.mean(x)) * (y - np.mean(y))) / (len(x) - 1)        return result    # 计算相关系数    @staticmethod    def corr_custom(x, y):        return Lm.cov_custom(x, y) / (np.std(x, ddof=1) * np.std(y, ddof=1))    # 计算回归系数    def simple_regression(self):        xmat = np.mat(self.xarr)        ymat = np.mat(self.yarr).T        xtx = xmat.T * xmat        if np.linalg.det(xtx) == 0.0:            print('Can not resolve the problem')            return        self.beta = np.linalg.solve(xtx, xmat.T * ymat)  # xtx.I * (xmat.T * ymat)        self.yhat = (xmat * self.beta).flatten().A[0]        return    # 计算拟合优度的判定系数R方，即相关系数corr的平方    def r_square(self):        y = np.mat(self.yarr)        ybar = np.mean(y)        self.r2 = np.sum((self.yhat - ybar) ** 2) / np.sum((y.A - ybar) ** 2)        return    # 计算估计标准误差    def estimate_deviation(self):        y = np.array(self.yarr)        self.se = np.sqrt(np.sum((y - self.yhat) ** 2) / self.dfd)        return    # 显著性检验F    def sig_test(self):        ybar = np.mean(self.yarr)        self.msr = np.sum((self.yhat - ybar) ** 2)        self.mse = np.sum((self.yarr - self.yhat) ** 2) / self.dfd        self.f = self.msr / self.mse        self.p = ss.f.sf(self.f, 1, self.dfd)        return    def summary(self):        self.simple_regression()        corr_coe = Lm.corr_custom(self.xarr[:, -1], self.yarr)        self.r_square()        self.estimate_deviation()        self.sig_test()        print('The Pearson\'s correlation coefficient: %.3f' % corr_coe)        print('The Regression Coefficient: %s' % self.beta.flatten().A[0])        print('R square: %.3f' % self.r2)        print('The standard error of estimate: %.3f' % self.se)        print('F-statistic:  %d on %s and %s DF,  p-value: %.3e' % (self.f, 1, self.dfd, self.p))

python执行结果：

The Regression Coefficient: [-87.51666667   3.45      ]R square: 0.991The standard error of estimate: 1.525F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

其中求回归系数时用矩阵转置求逆再用numpy内置的解线性方程组的方法是最快的：

a = np.mat(women.xarr); b = np.mat(women.yarr).Ttimeit (a.I * b)99.9 µs ± 941 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)timeit ata.I * (a.T*b)64.9 µs ± 717 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)timeit np.linalg.solve(ata, a.T*b)15.1 µs ± 126 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)