多元线性回归最小二乘法及其应用
下面的Java代码由网友设计,采用第一种方法计算参数。
1 /** 2 * 多元线性回归分析 3 * 4 * @param x[m][n] 5 * 每一列存放m个自变量的观察值 6 * @param y[n] 7 * 存放随即变量y的n个观察值 8 * @param m 9 * 自变量的个数 10 * @param n 11 * 观察数据的组数 12 * @param a 13 * 返回回归系数a0,...,am 14 * @param dt[4] 15 * dt[0]偏差平方和q,dt[1] 平均标准偏差s dt[2]返回复相关系数r dt[3]返回回归平方和u 16 * @param v[m] 17 * 返回m个自变量的偏相关系数 18 */ 19 public static void sqt2(double[][] x, double[] y, int m, int n, double[] a, 20 double[] dt, double[] v) { 21 int i, j, k, mm; 22 double q, e, u, p, yy, s, r, pp; 23 double[] b = new double[(m + 1) * (m + 1)]; 24 mm = m + 1; 25 b[mm * mm - 1] = n; 26 for (j = 0; j <= m - 1; j++) { 27 p = 0.0; 28 for (i = 0; i <= n - 1; i++) 29 p = p + x[j][i]; 30 b[m * mm + j] = p; 31 b[j * mm + m] = p; 32 } 33 for (i = 0; i <= m - 1; i++) 34 for (j = i; j <= m - 1; j++) { 35 p = 0.0; 36 for (k = 0; k <= n - 1; k++) 37 p = p + x[i][k] * x[j][k]; 38 b[j * mm + i] = p; 39 b[i * mm + j] = p; 40 } 41 a[m] = 0.0; 42 for (i = 0; i <= n - 1; i++) 43 a[m] = a[m] + y[i]; 44 for (i = 0; i <= m - 1; i++) { 45 a[i] = 0.0; 46 for (j = 0; j <= n - 1; j++) 47 a[i] = a[i] + x[i][j] * y[j]; 48 } 49 chlk(b, mm, 1, a); 50 yy = 0.0; 51 for (i = 0; i <= n - 1; i++) 52 yy = yy + y[i] / n; 53 q = 0.0; 54 e = 0.0; 55 u = 0.0; 56 for (i = 0; i <= n - 1; i++) { 57 p = a[m]; 58 for (j = 0; j <= m - 1; j++) 59 p = p + a[j] * x[j][i]; 60 q = q + (y[i] - p) * (y[i] - p); 61 e = e + (y[i] - yy) * (y[i] - yy); 62 u = u + (yy - p) * (yy - p); 63 } 64 s = Math.sqrt(q / n); 65 r = Math.sqrt(1.0 - q / e); 66 for (j = 0; j <= m - 1; j++) { 67 p = 0.0; 68 for (i = 0; i <= n - 1; i++) { 69 pp = a[m]; 70 for (k = 0; k <= m - 1; k++) 71 if (k != j) 72 pp = pp + a[k] * x[k][i]; 73 p = p + (y[i] - pp) * (y[i] - pp); 74 } 75 v[j] = Math.sqrt(1.0 - q / p); 76 } 77 dt[0] = q; 78 dt[1] = s; 79 dt[2] = r; 80 dt[3] = u; 81 } 82 private static int chlk(double[] a, int n, int m, double[] d) { 83 int i, j, k, u, v; 84 if ((a[0] + 1.0 == 1.0) || (a[0] < 0.0)) { 85 System.out.println("fail\n"); 86 return (-2); 87 } 88 a[0] = Math.sqrt(a[0]); 89 for (j = 1; j <= n - 1; j++) 90 a[j] = a[j] / a[0]; 91 for (i = 1; i <= n - 1; i++) { 92 u = i * n + i; 93 for (j = 1; j <= i; j++) { 94 v = (j - 1) * n + i; 95 a[u] = a[u] - a[v] * a[v]; 96 } 97 if ((a[u] + 1.0 == 1.0) || (a[u] < 0.0)) { 98 System.out.println("fail\n"); 99 return (-2); 100 } 101 a[u] = Math.sqrt(a[u]); 102 if (i != (n - 1)) { 103 for (j = i + 1; j <= n - 1; j++) { 104 v = i * n + j; 105 for (k = 1; k <= i; k++) 106 a[v] = a[v] - a[(k - 1) * n + i] * a[(k - 1) * n + j]; 107 a[v] = a[v] / a[u]; 108 } 109 } 110 } 111 for (j = 0; j <= m - 1; j++) { 112 d[j] = d[j] / a[0]; 113 for (i = 1; i <= n - 1; i++) { 114 u = i * n + i; 115 v = i * m + j; 116 for (k = 1; k <= i; k++) 117 d[v] = d[v] - a[(k - 1) * n + i] * d[(k - 1) * m + j]; 118 d[v] = d[v] / a[u]; 119 } 120 } 121 for (j = 0; j <= m - 1; j++) { 122 u = (n - 1) * m + j; 123 d[u] = d[u] / a[n * n - 1]; 124 for (k = n - 1; k >= 1; k--) { 125 u = (k - 1) * m + j; 126 for (i = k; i <= n - 1; i++) { 127 v = (k - 1) * n + i; 128 d[u] = d[u] - a[v] * d[i * m + j]; 129 } 130 v = (k - 1) * n + k - 1; 131 d[u] = d[u] / a[v]; 132 } 133 } 134 return (2); 135 } 136 /** 137 * @param args 138 */ 139 public static void main(String[] args) { 140 // TODO Auto-generated method stub 141 /** 142 * 一元回归 143 */ 144 // int i; 145 // double[] dt=new double[6]; 146 // double[] a=new double[2]; 147 // double[] x={ 0.0,0.1,0.2,0.3,0.4,0.5, 148 // 0.6,0.7,0.8,0.9,1.0}; 149 // double[] y={ 2.75,2.84,2.965,3.01,3.20, 150 // 3.25,3.38,3.43,3.55,3.66,3.74}; 151 // SPT.SPT1(x,y,11,a,dt); 152 // System.out.println(""); 153 // System.out.println("a="+a[1]+" b="+a[0]); 154 // System.out.println("q="+dt[0]+" s="+dt[1]+" p="+dt[2]); 155 // System.out.println(" umax="+dt[3]+" umin="+dt[4]+" u="+dt[5]); 156 /** 157 * 多元回归 158 */ 159 int i; 160 double[] a = new double[4]; 161 double[] v = new double[3]; 162 double[] dt = new double[4]; 163 double[][] x = { { 1.1, 1.0, 1.2, 1.1, 0.9 }, 164 { 2.0, 2.0, 1.8, 1.9, 2.1 }, { 3.2, 3.2, 3.0, 2.9, 2.9 } }; 165 double[] y = { 10.1, 10.2, 10.0, 10.1, 10.0 }; 166 SPT.sqt2(x, y, 3, 5, a, dt, v); 167 for (i = 0; i <= 3; i++) 168 System.out.println("a(" + i + ")=" + a[i]); 169 System.out.println("q=" + dt[0] + " s=" + dt[1] + " r=" + dt[2]); 170 for (i = 0; i <= 2; i++) 171 System.out.println("v(" + i + ")=" + v[i]); 172 System.out.println("u=" + dt[3]); 173 }
View Code
下面的C++代码由网友提供,采用第二中方法计算系数。
#include<iostream>
#include<fstream>
#include<iomanip>
using namespace std;
void transpose(double **p1,double **p2,int m,int n);
void multipl(double **p1,double **p2,double **p3,int m,int n,int p);
void Inver(double **p1,double **p2,int n);
double SD(double **p1,double **p2,double **p3,double **p4,int m,int n);
double ST(double **p1,int m);
void de_allocate(double **data,int m);
int main() {
int row,col;
char filename[30];
double SDsum,STsum,F,R2;
cout<<"Input original data file: \n";
ifstream infile; //打开文件
cin>>filename;
infile.open(filename);
if(!infile) {
cout<<"Opening the file failed!\n";
exit(1);
}
infile>>row>>col; //读入文件中的行数和列数
double **matrix=new double*[row]; //为动态二维数组分配内存
double **X=new double*[row];
double **Y=new double*[row];
double **XT=new double*[col];
double **XTX=new double*[col];
double **XTXInv=new double*[col];
double **XTXInvXT=new double*[col];
double **B=new double*[col];
double **YE=new double*[row];
for(int i=0;i<row;i++) {matrix[i]=new double[col];X[i]=new double[col];Y[i]=new double[1];Y[i]=new double[1];YE[i]=new double[1];
}
for(int i=0;i<col;i++) {XT[i]=new double[row];XTX[i]=new double[2*col];/为什么必须分配2*col列空间而不是col?在矩阵求逆时,XTX变增广矩阵,列数变为原来2吧倍,跟求逆算法有关。XTXInv[i]=new double[col];XTXInvXT[i]=new double[row];B[i]=new double[1];
}
for(int i=0;i<row;i++)for(int j=0;j<col;j++)infile>>matrix[i][j];
infile.close();
for(int i=0;i<row;i++) { //提取1X和Y数组列X[i][0]=1;Y[i][0]=matrix[i][col-1];for(int j=0;j<col-1;j++)X[i][j+1]=matrix[i][j];
}
transpose(X,XT,row,col);
multipl(XT,X,XTX,col,row,col);
Inver(XTX,XTXInv,col);
multipl(XTXInv,XT,XTXInvXT,col,col,row);
multipl(XTXInvXT,Y,B,col,row,1);
SDsum=SD(Y,X,B,YE,row,col);
STsum=ST(Y,row);
F=((STsum-SDsum)/(col-1))/(SDsum/(row-col));
R2=1/(1+(row-col)/F/(col-1));
cout<<"输出B:\n"; //屏幕输出结果B,SD,ST,F,R2
for(int i=0;i<col;i++)cout<<setiosflags(ios::fixed)<<setprecision(4)<<B[i][0]<<' ';cout<<endl;
cout<<"SD="<<SDsum<<';'<<"ST="<<STsum<<';'<<"F="<<F<<';'<<"R2="<<R2<<endl;
ofstream outfile; // 结果写入文件
cout<<"Output file'name:\n";
cin>>filename;
outfile.open(filename);
if(!outfile) {cout<<"Opening the file failed!\n";exit(1);
}
outfile<<"输出B:\n";
for(int i=0;i<col;i++)outfile<<B[i][0]<<' ';outfile<<endl;
outfile<<setiosflags(ios::fixed)<<setprecision(4)<<"SD="<<SDsum<<';'<<"ST="<<STsum<<';'<<"F="<<F<<';'<<"R2="<<R2<<endl;
outfile<<"Y and YE and Y-YE's value are:\n";
for(int i=0;i<row;i++)outfile<<Y[i][0]<<" "<<YE[i][0]<<" "<<Y[i][0]-YE[i][0]<<endl;
outfile.close();
de_allocate(matrix,row);
de_allocate(X,row);
de_allocate(Y,row);
de_allocate(XT,col);
de_allocate(XTX,col);
de_allocate(XTXInv,col);
de_allocate(XTXInvXT,col);
de_allocate(B,col);
de_allocate(YE,row);
system("pause");
return(0);
}
void de_allocate(double **data,int m) { //释放内存单元for(int i=0;i<m;i++)delete []data[i];delete []data;
}
double ST(double **p1,int m) { //求总离差平方和STdouble sum1=0,sum2=0,Yave=0;for(int i=0;i<m;i++)sum1+=p1[i][0];Yave=sum1/m;for(int i=0;i<m;i++)sum2+=(p1[i][0]-Yave)*(p1[i][0]-Yave);return sum2;
}
double SD(double **p1,double **p2,double **p3,double **p4,int m,int n) { //求偏差平方和SDdouble sum1=0,sum2=0;for(int i=0;i<m;i++) {sum1=0;for(int k=0;k<n;k++)sum1+=p2[i][k]*p3[k][0];p4[i][0]=sum1;}for(int i=0;i<m;i++)sum2+=(p1[i][0]-p4[i][0])*(p1[i][0]-p4[i][0]);return sum2;
}
void transpose(double **p1,double **p2,int m,int n) { //矩阵转置
for(int i=0;i<n;i++)for(int j=0;j<m;j++)p2[i][j]=p1[j][i];
}
void multipl(double **p1,double **p2,double **p3,int m,int n,int p) { //矩阵相乘
double sum;
for(int i=0;i<m;i++) {for(int j=0;j<p;j++) {sum=0;for(int k=0;k<n;k++)sum+=p1[i][k]*p2[k][j];p3[i][j]=sum;}
}
}
void Inver(double **p1,double **p2,int n) { //求逆矩阵
//初始化矩阵在右侧加入单位阵for(int i=0;i<n;i++) {for(int j=0;j<n;j++) {p1[i][j+n]=0;p1[i][i+n]=1;}}
//对于对角元素为0的进行换行操作for(int i=0;i<n;i++){while(p1[i][i]==0){for(int j=i+1;j<n;j++){if (p1[j][i]!=0){ double temp=0;for(int r=i;r<2*n;r++){temp=p1[j][r];p1[j][r]=p1[i][r];p1[i][r]=temp;}}break;}}//if (p1[i][i]==0) return 0;}//行变换为上三角矩阵double k=0;for(int i=0;i<n;i++) {for(int j=i+1;j<n;j++) {k=(-1)*p1[j][i]/p1[i][i];for(int r=i;r<2*n;r++)p1[j][r]+=k*p1[i][r];}}//行变换为下三角矩阵//double k=0;for(int i=n-1;i>=0;i--) {for(int j=i-1;j>=0;j--) {k=(-1)*p1[j][i]/p1[i][i];for(int r=0;r<2*n;r++)p1[j][r]+=k*p1[i][r];}}//化为单位阵for(int i=n-1;i>=0;i--) {k=p1[i][i];for(int j=0;j<2*n;j++)p1[i][j]/=k;}//拆分出逆矩阵for(int i=0;i<n;i++) {for(int j=0;j<n;j++)p2[i][j]=p1[i][n+j];}
}
二元线性回归最小二乘法拟合空间直线。网友提供
#include "stdafx.h"
using namespace std;
using std::cout;
using std::cin;
using std::endl;
#include<math.h>
#include <windows.h>
/*
这是一个控制台程序,任何一个空间直线方程都能以如下的方式表示
x=az+b
y=cz+d
z=z
即
(x-b)/a=(y-d)/c=z/1
*/
#include <vector>
using std::vector;
/* the fitting vaule of this line are: a=2 b=3 c=3 d=0
double z[10]={0.5,0.7,1,1.2,1.5,1.8,2,2.5,2.8,3};
double y[10]={1.5,2.1,3,3.6,4.5,5.4,6,7.5,8.4,9};
double x[10]={4,4.4,5,5.4,6,6.6,7,8,8.6,9};
*/
double x[]={1,1.5,2,2.5,3,3.5,4,4.5,5};
double y[]={-8.1,-7.2,-6.2,-5.5,-4.8,-3.8,-3,-2.2,-1.3};
double z[]={-12,-11.8,-10.7,-9.5,-8.2,-7,-6,-4.5,-3.5};
int _tmain(int argc, _TCHAR* argv[])
{vector<double>position_z;vector<double>position_y;vector<double>position_x;if ((sizeof(z)/sizeof(double))!=(sizeof(x)/sizeof(double))||(sizeof(z)/sizeof(double))!=(sizeof(y)/sizeof(double))||(sizeof(x)/sizeof(double))!=(sizeof(y)/sizeof(double))){::MessageBox(NULL,"请检查输入数组的长度是否相等","方程无解!",MB_OK);}//int ss=sizeof(z)/sizeof(double);for (int i=0;i<(sizeof(z)/sizeof(double));++i){position_z.push_back(z[i]);position_y.push_back(y[i]);position_x.push_back(x[i]);}//double m = position_z.size();//caculate o1,p1,c1,o2,p2,c2double o1 = 0.0;double p1 = 0.0;double x1 = 0.0;double o2 = 0.0;double p2 = 0.0;double x2 = 0.0;double y1 = 0.0;double y2 = 0.0;for (vector<double>::size_type t=0;t< position_z.size();++t){o1 += position_z[t]*position_z[t];p1 += position_z[t];x1 += position_x[t]*position_z[t];o2 += position_z[t];p2 = position_z.size();x2 += position_x[t] ;y2 += position_y[t];y1 += position_y[t]*position_z[t];}//caculate a bdouble a = 0.0 ,b = 0.0, c = 0.0, d = 0.0 ;if ((o1*p2-o2*p1)!= 0){a = (x1*p2 - x2*p1) /(o1*p2-o2*p1);if (a<1e-10){a=0;}b = (x2*o1 - x1*o2) /(o1*p2-o2*p1);if (b<1e-10){b=0;}c = (y1*p2 - y2*p1) /(o1*p2-o2*p1);if (c<1e-10){c=0;}d = (y2*o1 - y1*o2) /(o1*p2-o2*p1);if (d<1e-10){d=0;}}else{::MessageBox(NULL,"OK","方程无解!",MB_OK);}cout<<a<<" " <<b<<" "<<c<<" "<<d<<endl;system ("pause");return 0;
}
转载于:https://www.cnblogs.com/star91/p/4750102.html
多元线性回归最小二乘法及其应用相关推荐
- 波士顿房价的三种预测方式(模型预测,最小二乘法,多元线性回归)
首先导入库并避免显示错误FutureWorning import numpy as np import matplotlib.pyplot as plt import pandas as pd imp ...
- c++多元线性回归_五种优化算法实现多元线性回归
实现多元线性回归的要求及假设条件: '''线性回归的假设条件:1.样本独立,即每个预测样本之间没有依赖关系:2.残差e要服从正态分布,即y_true-y_pred的残差需要服从高斯分布:3.特征之间独 ...
- UA MATH571A 多元线性回归V 自相关与非线性模型简介
UA MATH571A 多元线性回归V 自相关与非线性模型简介 一阶误差自相关模型 Durbin-Watson检验 一阶自相关的消去 Cochrane-Orcutt方法 Hildreth-Lu方法 非 ...
- UA MATH571A 多元线性回归I 模型设定与推断
UA MATH571A 多元线性回归I 模型设定与推断 模型设定 最小二乘法(Method of Least Square) 系数 Mean Response and Residual 多元回归的AN ...
- 机器学习 回归篇(1)——多元线性回归
机器学习 回归篇(1)--多元线性回归 摘要 线性回归简介 python实现 运行结果及可视化 摘要 本文介绍了最基础的回归问题--多元线性回归,并通过python进行实现及可视化展示运行结果. 线性 ...
- 机器学习基础-多元线性回归-02
矩阵运算 多元线性回归 梯度下降法-多元线性回归 import numpy as np from numpy import genfromtxt import matplotlib.pyplot as ...
- 线性模型(1) —— 多元线性回归
提纲: 线性模型的基本形式 多元线性回归的损失函数 最小二乘法求多元线性回归的参数 最小二乘法和随机梯度下降的区别 疑问 学习和参考资料 1.线性模型的基本形式 线性模型是一种形式简单,易于建模,且可 ...
- 《计量经济学》学习笔记之多元线性回归模型
导航 上一章:一元线性回归模型 下一章:放宽基本假定的模型 文章目录 导航 3.1多元线性回归模型 一.多元线性回归模型 二.多元线性回归的基本假设 3.2多元线性回归模型的参数估计 四.参数统计量的 ...
- 第十二章 多元线性回归
1 多元线性回归模型 1 多元回归模型与回归方程 多元回归模型: y=β 0 +β 1 x 1 +β 2 x 2 +...+β k x k +ε y=\beta_0+\beta_1x_1+\beta_ ...
最新文章
- [转载]PhotoShop性能优化
- Kaggle心得(一)
- AE教程:学会这个,你做的Logo就可以单独出道了
- 2020-03-21
- 计算机应用基础简单实操,浅谈《计算机应用基础》实操课的教学管理
- 如何跨服务器复制表中数据
- freeswitch 用户配置详解_FreeSwitch安装配置记录-阿里云开发者社区
- jquery的extend和fn.extend
- 树、森林和二叉树之间的转换
- 微信公众号管理平台使用教程
- vscode运行C程序
- windows7最简单最快速解决“此windows副本不是正版”(“This copy of Windows is not genuine”)方法
- 数组和ArrayList的区别
- 关于vscode中输入的中文变繁体的问题
- 关于Windows命令提示符中的 xxx > nul 2 > nul
- 国际期刊预警名单网址
- 用马悦凌的养生方法--减肥
- 一个简单的登录注册界面流程介绍
- 无线pda是快递员随身携带的设备
- Python安装配置: python install python安装
热门文章
- 阿里大师总结的Web安全超全知识点,看这一篇就够了
- 数据库中表的连接方式详解
- 数字图像处理之图像的梯度倒数加权平滑:用Java+OpenCV实现,附源码
- PCSE_NASA【wofost】 天气数据获取(5.5新版本)
- 【OpenGL】GLFW创建三角形
- 第三方支付宝支付(非真实金额支付)
- 如何写一个python程序浏览淘宝_一篇文章教会你用Python爬取淘宝评论数据(写在记事本)...
- 1096 大美数 (15 分)(测试点有坑)
- 使用KALI破解WIFI(wpa/wpa2)密码
- overwrite java_java中的overwrite怎么用?最好是有代码的