y=c(160,260,210,265,240,220,275,160,275,250)

x1=c(70,75,65,74,72,68,78,66,70,65)

x2=c(35,40,40,42,38,45,42,36,44,42)

x3=c(1,2.4,2,3,1.2,1.5,4,2,3.2,3)

shuju<-data.frame(y,x1,x2,x3)

shuju



cor(shuju)   #(1)变量间相关系数阵



shuju.reg=lm(y~.,data=shuju)#(2)求三元线性回归方程(3)拟合优度检验

summary(shuju.reg)#(3)拟合优度检验(4)(5)回归方程、回归系数显著性检验



shuju.reg1=lm(y~x1+x2,data=shuju)#(6)剔除不显著的回归系数

summary(shuju.reg1)





stdshuju<-scale(shuju, center=T,scale=T)#对各列数据进行标准化

shuju.stdreg <- lm(y~x1+x2,data=data.frame(stdshuju))

summary(shuju.stdreg)#(8)求标准化回归方程



confint(shuju.reg1)   #(7)每一个回归系数的区间估计



shuju.pred1<-predict(shuju.reg1,newdata=data.frame(x1=75,x2=42))  #(9)预测y0

shuju.pred1   #(9)预测y0

shuju.pred2<-predict(shuju.reg1,newdata=data.frame(x1=75,x2=42),interval='prediction') 

shuju.pred2    #(9)y0的精确置信区间

267.829-2*24.08；267.829+2*24.08   #(9)手工计算y0的近似预测区间







y.res=resid(shuju.reg)    #残差值

y.fit=predict(shuju.reg)  #y拟合值

plot(y.res~y.fit)         #残差序列图



m = c(-15.47480708,12.82498933,5.34434481,-0.09088064,33.22548824,-25.19759251,-17.55449753,-20.00684234,8.23434759,18.69545013 )

acf(m)   #残差的各阶自相关系数



library(lmtest)

dwtest(shuju.reg)   #做DW检验

我正在使用在
scipy Cookbook下载的ols.py代码(下载在第一段用粗体OLS)但我需要了解而不是使用ols函数的随机数据来进行多元线性回归.
我有一个特定的因变量y和三个解释变量.每当我尝试将变量放入随机变量时,它就会给出错误：
TypeError: this constructor takes no arguments.
有人可以帮忙吗？这可能吗？
这是我尝试使用的ols代码的副本以及我尝试输入的变量
from __future__ import division
from scipy import c_, ones, dot, stats, diff
from scipy.linalg import inv, solve, det
from numpy import log, pi, sqrt, square, diagonal
from numpy.random import randn, seed
import time
class ols:
"""
Author: Vincent Nijs (+ ?)
Email: v-nijs at kellogg.northwestern.edu
Last Modified: Mon Jan 15 17:56:17 CST 2007
Dependencies: See import statement at the top of this file
Doc: Class for multi-variate regression using OLS
Input:
dependent variable
y_varnm = string with the variable label for y
x = independent variables, note that a constant is added by default
x_varnm = string or list of variable labels for the independent variables
Output:
There are no values returned by the class. Summary provides printed output.
All other measures can be accessed as follows:
Step 1: Create an OLS instance by passing data to the class
m = ols(y,x,y_varnm = 'y',x_varnm = ['x1','x2','x3','x4'])
Step 2: Get specific metrics
To print the coefficients:
>>> print m.b
To print the coefficients p-values:
>>> print m.p
"""
y = [29.4, 29.9, 31.4, 32.8, 33.6, 34.6, 35.5, 36.3, 37.2, 37.8, 38.5, 38.8,
38.6, 38.8, 39, 39.7, 40.6, 41.3, 42.5, 43.9, 44.9, 45.3, 45.8, 46.5,
77.1, 48.2, 48.8, 50.5, 51, 51.3, 50.7, 50.7, 50.6, 50.7, 50.6, 50.7]
#tuition
x1 = [376, 407, 438, 432, 433, 479, 512, 543, 583, 635, 714, 798, 891,
971, 1045, 1106, 1218, 1285, 1356, 1454, 1624, 1782, 1942, 2057, 2179,
2271, 2360, 2506, 2562, 2700, 2903, 3319, 3629, 3874, 4102, 4291]
#research and development
x2 = [28740.00, 30952.00, 33359.00, 35671.00, 39435.00, 43338.00, 48719.00, 55379.00, 63224.00,
72292.00, 80748.00, 89950.00, 102244.00, 114671.00, 120249.00, 126360.00, 133881.00, 141891.00,
151993.00, 160876.00, 165350.00, 165730.00, 169207.00, 183625.00, 197346.00, 212152.00, 226402.00,
267298.00, 277366.00, 276022.00, 288324.00, 299201.00, 322104.00, 347048.00, 372535.00,
397629.00]
#one/none parents
x3 = [11610, 12143, 12486, 13015, 13028, 13327, 14074, 14094, 14458, 14878, 15610, 15649,
15584, 16326, 16379, 16923, 17237, 17088, 17634, 18435, 19327, 19712, 21424, 21978,
22684, 22597, 22735, 22217, 22214, 22655, 23098, 23602, 24013, 24003, 21593, 22319]
def __init__(self,y,x1,y_varnm = 'y',x_varnm = ''):
"""
Initializing the ols class.
"""
self.y = y
#self.x1 = c_[ones(x1.shape[0]),x1]
self.y_varnm = y_varnm
if not isinstance(x_varnm,list):
self.x_varnm = ['const'] + list(x_varnm)
else:
self.x_varnm = ['const'] + x_varnm
# Estimate model using OLS
self.estimate()
def estimate(self):
# estimating coefficients, and basic stats
self.inv_xx = inv(dot(self.x.T,self.x))
xy = dot(self.x.T,self.y)
self.b = dot(self.inv_xx,xy)                    # estimate coefficients
self.nobs = self.y.shape[0]                     # number of observations
self.ncoef = self.x.shape[1]                    # number of coef.
self.df_e = self.nobs - self.ncoef              # degrees of freedom, error
self.df_r = self.ncoef - 1                      # degrees of freedom, regression
self.e = self.y - dot(self.x,self.b)            # residuals
self.sse = dot(self.e,self.e)/self.df_e         # SSE
self.se = sqrt(diagonal(self.sse*self.inv_xx))  # coef. standard errors
self.t = self.b / self.se                       # coef. t-statistics
self.p = (1-stats.t.cdf(abs(self.t), self.df_e)) * 2    # coef. p-values
self.R2 = 1 - self.e.var()/self.y.var()         # model R-squared
self.R2adj = 1-(1-self.R2)*((self.nobs-1)/(self.nobs-self.ncoef))   # adjusted R-square
self.F = (self.R2/self.df_r) / ((1-self.R2)/self.df_e)  # model F-statistic
self.Fpv = 1-stats.f.cdf(self.F, self.df_r, self.df_e)  # F-statistic p-value
def dw(self):
"""
Calculates the Durbin-Waston statistic
"""
de = diff(self.e,1)
dw = dot(de,de) / dot(self.e,self.e);
return dw
def omni(self):
"""
Omnibus test for normality
"""
return stats.normaltest(self.e)
def JB(self):
"""
Calculate residual skewness, kurtosis, and do the JB test for normality
"""
# Calculate residual skewness and kurtosis
skew = stats.skew(self.e)
kurtosis = 3 + stats.kurtosis(self.e)
# Calculate the Jarque-Bera test for normality
JB = (self.nobs/6) * (square(skew) + (1/4)*square(kurtosis-3))
JBpv = 1-stats.chi2.cdf(JB,2);
return JB, JBpv, skew, kurtosis
def ll(self):
"""
Calculate model log-likelihood and two information criteria
"""
# Model log-likelihood, AIC, and BIC criterion values
ll = -(self.nobs*1/2)*(1+log(2*pi)) - (self.nobs/2)*log(dot(self.e,self.e)/self.nobs)
aic = -2*ll/self.nobs + (2*self.ncoef/self.nobs)
bic = -2*ll/self.nobs + (self.ncoef*log(self.nobs))/self.nobs
return ll, aic, bic
def summary(self):
"""
Printing model output to screen
"""
# local time & date
t = time.localtime()
# extra stats
ll, aic, bic = self.ll()
JB, JBpv, skew, kurtosis = self.JB()
omni, omnipv = self.omni()
# printing output to screen
print                                                                                         '\n=============================================================================='
print "Dependent Variable: " + self.y_varnm
print "Method: Least Squares"
print "Date: ", time.strftime("%a, %d %b %Y",t)
print "Time: ", time.strftime("%H:%M:%S",t)
print '# obs:               %5.0f' % self.nobs
print '# variables:     %5.0f' % self.ncoef
print '=============================================================================='
print 'variable     coefficient     std. Error      t-statistic     prob.'
print '=============================================================================='
for i in range(len(self.x_varnm)):
print '''% -5s          % -5.6f     % -5.6f     % -5.6f     % -5.6f'''  %            tuple([self.x_varnm[i],self.b[i],self.se[i],self.t[i],self.p[i]])
print '=============================================================================='
print 'Models stats                         Residual stats'
print '=============================================================================='
print 'R-squared            % -5.6f         Durbin-Watson stat  % -5.6f' %              tuple([self.R2, self.dw()])
print 'Adjusted R-squared   % -5.6f         Omnibus stat        % -5.6f' % tuple([self.R2adj, omni])
print 'F-statistic          % -5.6f         Prob(Omnibus stat)  % -5.6f' % tuple([self.F, omnipv])
print 'Prob (F-statistic)   % -5.6f         JB stat             % -5.6f' % tuple([self.Fpv, JB])
print 'Log likelihood       % -5.6f         Prob(JB)            % -5.6f' % tuple([ll, JBpv])
print 'AIC criterion        % -5.6f         Skew                % -5.6f' % tuple([aic, skew])
print 'BIC criterion        % -5.6f         Kurtosis            % -5.6f' % tuple([bic, kurtosis])
print '=============================================================================='
if __name__ == '__main__':
##########################
### testing the ols class
##########################
# intercept is added, by default
m = ols(y,x1,y_varnm = 'y',x_varnm = ['x1','x2','x3'])
m.summary()

% 多元线性回归
% b 回归系数
% bint 回归系数的置信度为95％的置信区间
% r 残差
% rint 各残差的置信区间
% stats 用于检验回归模型的统计量，有四个数值：相关系数R^2、F值、与F对应的概率p，误差方差。相关系数R^2越接近1，说明回归方程越显著；F > F1-α(k，n-k-1)时拒绝H0，F越大，说明回归方程越显著；与F对应的概率p 时拒绝H0，回归模型成立。p值在0.01-0.05之间，越小越好。
clc
clear
x = xlsread('D:\数学建模\8.28\data1.xlsx','Sheet1','B3:I35');
y = xlsread('D:\数学建模\8.28\data1.xlsx','Sheet2','B3:B35');
X = [ones(length(y),1),x]; % 把行向量转秩为列向量(此处不用转)
Y = y; % 把行向量转秩为列向量(此处不用转)
[b,bint,r,rint,stats] = regress(Y,X); %求回归系数的点估计和区间估计，并检验回归模型
b
bint
stats
% 残差分析
subplot(1,2,1);
rcoplot(r,rint);
%画图
subplot(1,2,2);
b2=[b(2),b(3),b(4),b(5),b(6),b(7),b(8),b(9)];
z = b(1)+b2*x';
plot(X,Y,'k+',X,z,'r');
grid on
% 预测数据
tx = [120731808012512581.190];
ty=b(1)+b2*tx';
ty


前言
本文的实现主要是参照了Andrew NG的机器学习课程所讲的内容。理论知识上一篇博文已经介绍。由于刚接触C++，代码写得比较粗糙，望见谅。
实验环境
Visual Studio 2013
OpenCV 2.4
实验代码
1.定义一个线性回归的类Regression：
头文件 Regression.h
#ifndef  _Regression_H_
#define  _Regression_H_
#include <vector>
#include <string>
#include <iostream>
#include <strstream>
#include <fstream>
#include <opencv2\core\core.hpp>
using namespace cv;
using namespace std;
class Regression
{
public:
	//加载数据
	void loadData(string filePath, string regx);
	//数据返回
	Mat& getX();
	Mat& getEX();
	Mat& getTheta();
	Mat& getY();
	//特征归一化
	void featureNormalize(Mat& x);
	//预测
	float predict(Mat& x);
	//假设函数，此处设成静态函数，是方便在其他类中直接调用
    static void calculateHx(Mat& x, Mat& theta, Mat& Hx);
private:
	Mat originalX;
	Mat extendX;
	Mat y;
	int dim;//数据维数
	Mat theta;//列向量
};
#endif
实现文件 Regression.cpp
#include "Regression.h"

//把字符串转换成数值类型
template<typename T>
T convertStringData(string data)
{
	strstream ss;
	T y;
	ss <<data;
	ss >> y;
	return y;
}

//字符串分割函数 
inline vector<string> split(string str, string pattern)
{
	string::size_type pos;
	vector<string> result;
	str += pattern;//扩展字符串以方便操作
	int size = str.size();

	for (int i = 0; i<size; i++)
	{
		pos = str.find(pattern, i);
		if (pos<size)
		{
			std::string s = str.substr(i, pos - i);
			result.push_back(s);
			i = pos + pattern.size() - 1;
		}
	}
	return result;
}

//加载数据
void Regression::loadData(string filePath,string regx)
{
	ifstream fin(filePath);
	while (!fin.eof())
	{
		string temp;
		getline(fin, temp);
		if (strcmp(temp.c_str(), "") == 0) continue;
		vector<string> resStr = split(temp, regx);
		//设置x的维数
		dim = resStr.size() - 1;
		float yi=convertStringData<float>(resStr[resStr.size() - 1]);
		y.push_back(yi);
		Mat data = Mat::zeros(1, dim,CV_32FC1);//临时存放每一行的数据
		for (int i = 0; i < dim; i++)
		{
			data.row(0).col(i) = convertStringData<float>(resStr[i]);;
		}
		originalX.push_back(data);
	}
	extendX = Mat::ones(originalX.rows, dim + 1, CV_32FC1);
	originalX.copyTo(extendX.colRange(1, extendX.cols));
	extendX.col(0) = Mat::ones(extendX.rows, 1, extendX.type());
	theta = Mat::zeros(dim+1,1,originalX.type());
	fin.close();
}


//数据返回
Mat& Regression::getX()
{
	return originalX;
}

Mat& Regression::getY()
{
	return y;
}

Mat& Regression::getEX()
{
	return extendX;
}

Mat& Regression::getTheta()
{
	return theta;
}

//特征归一化
void Regression::featureNormalize(Mat& ex)
{
	if (ex.cols <= 2) return;//只有一个特征的时候不用归一化
	for (int col = 1; col < ex.cols; col++)
	{
		Mat mean;
		Mat stddev;
		meanStdDev(ex.col(col), mean, stddev);
		//归一化
		ex.col(col) = ex.col(col) - mean;
		ex.col(col) = ex.col(col) / stddev;
	}
}


//预测
//x:一个样本
//theta 训练得到的参数
float Regression::predict(Mat& x)
{
	Mat hx;
	calculateHx(x, this->getTheta(), hx);
	return hx.at<float>(0,0);
}

//计算假设函数
void Regression::calculateHx(Mat& x, Mat& theta, Mat& Hx)
{
	Hx = x*theta;
}

2.定义一个梯度下降的类GradientDescent：

头文件 GradientDescent.h
#ifndef _GRADIENTDESCENT_H
#define _GRADIENTDESCENT_H
#include <opencv2\core\core.hpp>
using namespace cv;
class GradientDescent
{
public:
	//批量梯度下降
	void gradientDescent(Mat& x, Mat& y, Mat &theta, int num_iters, float alpha);
	//计算代价
	float computeCost(Mat& x, Mat& y, Mat &theta);
	void normalizeEquation(Mat& x, Mat& y, Mat &theta);
};
#endif
实现文件 GradientDescent.cpp
#include "GradientDescent.h"
#include "Regression.h"
#include <iostream>
using namespace std;
void GradientDescent::gradientDescent(Mat& x, Mat& y, Mat &theta, int num_iters, float alpha)
{
	int dataSize = x.rows;//数据总数
	int dim = theta.rows;//theta的维数
	for (int i = 0; i < num_iters; i++)
	{
		//矢量化编程
		Mat hx;
		Regression::calculateHx(x, theta, hx);
		Mat thetaNew = theta-alpha*x.t()*(hx - y)/dataSize;
		thetaNew.copyTo(theta);
	}
}

float GradientDescent::computeCost(Mat& x, Mat& y, Mat &theta)
{
	float error = .0;
	Mat hx;
	Mat j;
	hx = x*theta;
	j = (hx - y).t()*(hx - y);
	error = j.at<float>(0, 0) / (2 * x.rows);
	return error;
}

void GradientDescent::normalizeEquation(Mat& x, Mat& y, Mat &theta)
{
	theta=(x.t()*x).inv()*x.t()*y;
}

3.代码测试：
样本集由Andrew NG机器学习课程提供，可在网上下载
main.cpp内容如下：
#include "Regression.h"
#include "GradientDescent.h"
#include <iostream>
#include <opencv2\opencv.hpp>
using namespace cv;
using namespace std;
int main()
{
	Regression regression;
	regression.loadData("test_data/ex1data1.txt",",");
	regression.featureNormalize(regression.getEX());
	GradientDescent gd;
	float error=gd.computeCost(regression.getEX(), regression.getY(), regression.getTheta());
	cout << "cost:" << error << endl;
	gd.gradientDescent(regression.getEX(), regression.getY(), regression.getTheta(), 1500, 0.01);
	//gd.normalizeEquation(regression.getEX(), regression.getY(), regression.getTheta());
	cout << regression.getTheta() << endl;
	getchar();
	return 0;
}
输出结果：

与matlab运行结果一致