機械学習/線形回帰/ソース - AIとか機械学習とか、それらのコンテストとかのwiki

#include <map>
#include <set>
#include <cmath>
#include <stack>
#include <queue>
#include <string>
#include <vector>
#include <bitset>
#include <fstream>
#include <sstream>
#include <stdio.h>
#include <ctype.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <sys/time.h>
using namespace std;
#define li        long long int
#define rep(i,to) for(li i=0;i<((li)(to));++i)
#define pb        push_back
#define sz(v)     ((li)(v).size())
#define bit(n)    (1ll<<(li)(n))
#define all(vec)  (vec).begin(),(vec).end()
#define each(i,c) for(__typeof((c).begin()) i=(c).begin();i!=(c).end();i++)
#define MP        make_pair
#define F         first
#define S         second
 
 
//	p.180
//
#define MAT_TYPE double
//#define MAT_MOD 
//
 
class MAT{
public:
	long long int n,m;
	vector<vector<MAT_TYPE> > mat;
 
	inline MAT();
	inline MAT(int _n);
	inline MAT(int _n,int _m,int num);
	inline void print();
	inline long long int size();
	inline MAT operator+(MAT mat0);
	inline MAT operator-(MAT mat0);
	inline MAT operator*(MAT mat0);
	inline MAT operator*(MAT_TYPE num);
	inline vector<MAT_TYPE>& operator[](int num){ return mat[num]; }
	friend MAT transpose(MAT mat0);
	friend MAT pow(MAT mat0,long long int a);
};
 
long long int MAT::size(){ return n; }
 
MAT transpose(MAT mat0){
	MAT res(mat0.m,mat0.n,0);
	for(int i=0;i<mat0.n;i++)for(int j=0;j<mat0.m;j++) res[j][i]=mat0[i][j];
	return res;
}
 
MAT pow(MAT mat0,long long int a){
	MAT res(mat0.n);
	for(;a;a>>=1){
		if(a&1) res=res*mat0;
		mat0=mat0*mat0;
	}
	return res;
}
 
MAT::MAT(){
	n=m=0;
}
 
MAT::MAT(int _n){
	n=m=_n;
	mat=vector<vector<MAT_TYPE> >(n,vector<MAT_TYPE>(n,0));
	for(int i=0;i<n;i++) mat[i][i]=1;
}
 
MAT::MAT(int _n,int _m,int num){
	n=_n;
	m=_m;
	mat=vector<vector<MAT_TYPE> >(n,vector<MAT_TYPE>(m,num));
}
 
inline MAT MAT::operator+(MAT mat0){
	MAT res=*this;
	for(int i=0;i<n;i++)for(int j=0;j<m;j++) res[i][j]=res[i][j]+mat0[i][j];
#ifdef MAT_MOD
	for(int i=0;i<n;i++)for(int j=0;j<m;j++)if(MAT_MOD<=res[i][j]) res[i][j]-=MAT_MOD;
#endif
	return res;
}
 
inline MAT MAT::operator-(MAT mat0){
	MAT res=*this;
	for(int i=0;i<n;i++)for(int j=0;j<m;j++) res[i][j]=res[i][j]-mat0[i][j];
#ifdef MAT_MOD
	for(int i=0;i<n;i++)for(int j=0;j<m;j++)if(res[i][j]<0) res[i][j]+=MAT_MOD;
#endif
	return res;
}
 
inline MAT MAT::operator*(MAT mat0){
	MAT res=MAT(n,mat0.m,0);
	for(int i=0;i<n;i++)for(int j=0;j<mat0.m;j++){
#ifdef MAT_MOD
		for(int k=0;k<m;k++) res[i][j]=(res[i][j]+mat[i][k]*mat0[k][j])%MAT_MOD;
#else
		for(int k=0;k<m;k++) res[i][j]=(res[i][j]+mat[i][k]*mat0[k][j]);
#endif
	}
	return res;
}
 
inline MAT MAT::operator*(MAT_TYPE num){
	MAT res=*this;
#ifdef MAT_MOD
	for(int i=0;i<n;i++)for(int j=0;j<m;j++) res[i][j]=(res[i][j]*num)%MAT_MOD;
#else
	for(int i=0;i<n;i++)for(int j=0;j<m;j++) res[i][j]=(res[i][j]*num);
#endif
	return res;
}
 
void MAT::print(){
	for(int i=0;i<n;i++){
		for(int j=0;j<m;j++) printf("%0.6lf ",(double)mat[i][j]);
		cout<<endl;
	}
}
 
 
class LU{
public:
	int n;
	MAT base,L,U;
	vector<int> swaped;
 
	LU(MAT mat);
	MAT solve(MAT mat);
	MAT solve(vector<MAT_TYPE> vec);
 
};
 
LU::LU(MAT mat){
	n=mat.n;
	base=mat;
	L=MAT(n);
	U=MAT(n,n,0);
	for(int i=0;i<n;i++) swaped.pb(i);
	for(int i=0;i<n;i++){
		for(int j=0;j<n;j++){
			double sum=0.0;
			for(int k=0;k<min(i,j);k++) sum+=L[i][k]*U[k][j];
			if(j<i) L[i][j]=(base[i][j]-sum)/U[j][j];
			else U[i][j]=base[i][j]-sum;
		}
		pair<double,int> maxi=MP(abs(U[i][i]),i);
		for(int j=i;j<n;j++) maxi=max(maxi,MP(abs(U[i][j]),j));
		swap(swaped[i],swaped[maxi.S]);
		for(int j=0;j<n;j++) swap(U[j][i],U[j][maxi.S]);
		for(int j=0;j<n;j++) swap(base[j][i],base[j][maxi.S]);
	}
}
 
MAT LU::solve(MAT mat){
	MAT ans(n,1,0);
	MAT res(n,1,0);
	for(int i=0;i<n;i++){
		double sum=mat[i][0];
		for(int j=0;j<i;j++) sum-=L[i][j]*ans[j][0];
		ans[i][0]=sum/L[i][i];
	}
	for(int i=n-1;0<=i;i--){
		double sum=ans[i][0];
		for(int j=n-1;i<j;j--) sum-=U[i][j]*ans[j][0];
		ans[i][0]=sum/U[i][i];
	}
	for(int i=0;i<n;i++) res[swaped[i]]=ans[i];
	return res;
}
 
MAT LU::solve(vector<MAT_TYPE> vec){
	MAT mat(vec.size(),1,0);
	for(int i=0;i<mat.n;i++) mat[i][0]=vec[i];
	return solve(mat);
}
 
 
/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
 
 
 
 
 
 
 
 
 
vector<double> x,y;
 
double func(double var){
	return 5.0*var*var*var*var+4.0*var*var*var-3.0*var*var+2.0*var+1.0;
}
 
int main(){
	rep(i,100){
		x.pb(rand()%1000);
		y.pb(func(x.back()));
	}
	MAT X(sz(x),5,0),Y(sz(x),1,0);
	rep(i,sz(x)){
		X[i][0]=1;
		rep(j,4) X[i][j+1]=X[i][j]*x[i];
	}
	rep(i,sz(y)) Y[i][0]=y[i];
	MAT TX=transpose(X);
	MAT RIGHT=TX*Y;
	LU LEFT(TX*X);
	MAT ans=LEFT.solve(RIGHT);
	rep(i,sz(ans)) printf("%lld : %0.6lf\n",i,ans[i][0]);
	double sum=0,SUM=0;
	rep(i,sz(x)) sum+=abs(X[i][4]*ans[4][0]+X[i][3]*ans[3][0]+X[i][2]*ans[2][0]+X[i][1]*ans[1][0]+X[i][0]*ans[0][0]-y[i]);
	rep(i,sz(y)) SUM+=y[i];
	printf("誤差 %0.20lf%%\n",sum/SUM);
}