利用python手动写最小二乘估计
2015-10-13 21:25
525 查看
今天朋友请我吃完晚饭回来,刚刚一个python群问起最小二乘回归估计的问题,他不知道在python里面怎么实现,我告诉他有很多方法去实现,比如说sklean库里面有关于最小二乘估计现成模块,利用pyper调用R软件里面lm,再或者利用sas里面的proc reg过程,最后还是不知道怎么弄,不得不说从原理到实现跟他讲了一遍,说最小二乘估计最后就是矩阵求逆和求转置相乘估计出来,最后利用numpy库根据原理自己写code重新实现一遍,数据是用的鸢尾花的数据(iris,R里面自带的),直接上程序:
估计出来的参数是不一模一样?
当存在共线性时可以用到岭回归和lasso回归来处理,如果不用到岭回归和lasso回归,当变量共线性时矩阵可能不存在逆,这时候参数根本就无法估计正常情况下,还好学过高等代数的都知道,矩阵里面存在一种求伪逆的方法,在Python中是np.linalg.pinv方法求矩阵的逆,如果矩阵是奇异矩阵,则可以通过np.linalg.pinv(x.T*x)的形式求出。
哈哈哈 希望能帮你们理解最小二乘估计原理, 上面过程是无截距项估计,要是要用到有截距项估计,在上面矩阵加入一列全部是1即可, 在lm方法中加入formula = Sepal.Length ~ 1+ Sepal.Width + Petal.Length, 即可,话说回来,在工业界做模型更多用的无截距项,R里面估计出来的结果信息比较多,t检验、f检验的结果都全部出来,但是自己写的最小二乘估计就差远了,洗澡。。。。睡觉去了,不说了
#I是矩阵求逆 ,T是表示矩阵转置在numpy库中 import numpy as np import pandas as pd test=pd.read_csv("C:\\Users\\Administrator\\Desktop\\iris.csv") ols_matrix_x=np.matrix(test.iloc[:,2:4],dtype=np.float64) ols_matrix_y=np.matrix(test.iloc[:,1],dtype=np.float64).T b=(ols_matrix_x.T*ols_matrix_x).I*ols_matrix_x.T*ols_matrix_y print u"参数项矩阵为{0}".format(b) i=0 parameter=[] while i<2: parameter.append(b[i,0]) i+=1 temp_e=ols_matrix_y-ols_matrix_x*b mye=temp_e.sum()/temp_e.size e=np.matrix([mye,mye,mye]).T print "%f*Sepal.Width+%f*Petal.Length"%(parameter[0],parameter[1]) 结果: 最小二乘法估计的参数为[[ 1.20291972] [ 0.56905789]] 1.202920*Sepal.Width+0.569058*Petal.Length 别急 利用pyper库调用R里面的lm方法进行验证比较,直接上代码: import pyper as pr import pandas as pd mpg=pd.read_csv("C:\\Users\\Administrator\\Desktop\\iris.csv") r=pr.R(RCMD="D:\\Program Files\\R\\R-3.1.3\\bin\\R",use_dict=True,use_pandas=True,use_numpy=True) r.assign("rmpg",mpg) # print r("summary(rmpg)") print r("colnames(rmpg)") r("result_lm<-lm(Sepal.Length~ 0+Sepal.Width + Petal.Length ,data=rmpg)") print r("summary(result_lm)") 结果是: try({colnames(rmpg)}) [1] "Unnamed..0" "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" [6] "Species" try({summary(result_lm)}) Call: lm(formula = Sepal.Length ~ 0 + Sepal.Width + Petal.Length, data = rmpg) Residuals: Min 1Q Median 3Q Max -1.08398 -0.26510 0.05398 0.34020 1.07735 Coefficients: Estimate Std. Error t value Pr(>|t|) Sepal.Width 1.20292 0.02233 53.86 <2e-16 *** Petal.Length 0.56906 0.01662 34.24 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.4148 on 148 degrees of freedom Multiple R-squared: 0.9951, Adjusted R-squared: 0.9951 F-statistic: 1.51e+04 on 2 and 148 DF, p-value: < 2.2e-16
估计出来的参数是不一模一样?
当存在共线性时可以用到岭回归和lasso回归来处理,如果不用到岭回归和lasso回归,当变量共线性时矩阵可能不存在逆,这时候参数根本就无法估计正常情况下,还好学过高等代数的都知道,矩阵里面存在一种求伪逆的方法,在Python中是np.linalg.pinv方法求矩阵的逆,如果矩阵是奇异矩阵,则可以通过np.linalg.pinv(x.T*x)的形式求出。
哈哈哈 希望能帮你们理解最小二乘估计原理, 上面过程是无截距项估计,要是要用到有截距项估计,在上面矩阵加入一列全部是1即可, 在lm方法中加入formula = Sepal.Length ~ 1+ Sepal.Width + Petal.Length, 即可,话说回来,在工业界做模型更多用的无截距项,R里面估计出来的结果信息比较多,t检验、f检验的结果都全部出来,但是自己写的最小二乘估计就差远了,洗澡。。。。睡觉去了,不说了
相关文章推荐
- 利用python中的pyquery库简单的抓取数据
- python中pandas库学习笔记
- python 中 exec、 eval、 execfile 和 compile 用法
- python中的pandas包的数据清洗能力
- python中knn算法实现
- 搜索引擎关键词抓取 以百度为例 python
- python中结巴分词快速入门
- Windows Opencv-3.0 + Python-2.7.10 配置(numpy-1.8.1-64位)
- python 中ggplot画图
- python调用R
- java中调用python的方法
- Python3笔记——IDE的选择
- python if __name__ == '__main__' 详解
- Caffe 多爱Python一丢丢,Cifar-10
- 在Anaconda中安装python包seaborn
- Pythonic到底是什么玩意儿?
- python学习--异常
- python列表和分片
- Ubuntu环境下安装python的flask
- Python3笔记——定义变量