模仿R语言c++ 向量类c 矩阵类matrix等（持续更新 欢迎指点）
2017-05-17 10:49
573 查看
头文件：
[cpp] view plain copy print?
class_of_R.h#ifndef CLASS_OF_R_H_
#define CLASS_OF_R_H_
#include<iostream>
using std::cout;
using std::endl;
using std::cin;
using std::istream;
using std::ios_base;
#include<vector>
using std::vector;
#include<initializer_list>
using std::initializer_list;
#include<algorithm>
using std::for_each;
using std::sort;
using std::copy;
using std::find;
#include<stdexcept>
using std::runtime_error;
using std::out_of_range;
#include<cstdlib>
using std::srand;
using std::rand;
using std::system;
#include<ctime>
using std::time;
#include<cmath>
using std::sqrt;
using std::log;
using std::exp;
#include<utility>
using std::pair;
using std::make_pair;
#include<iomanip>
using std::setw;
#include<set>
using std::multiset;
using std::set;
#include<unordered_map>
using std::unordered_map;
using std::unordered_multimap;
#include<map>
using std::multimap;
using std::map;
#include<string>
using std::string;
#include<fstream>
using std::ifstream;
#include<memory>
using std::shared_ptr;
#include<list>
using std::list;
class matrix;
//类c声明
class c
{
private:
vector<double> c_val;
size_t dim;
public:
c():c_val(vector<double>()),dim(0){}
c(initializer_list<double>lst);
c(const vector<double>&c_v):c_val(c_v),dim(c_v.size()){}
c(const matrix&m);
void out(size_t n=1)const{cout<<“The vector is:\n”;for(size_t i=0;i<dim;i++)cout<<setw(n)<<c_val[i]<<“ ”;cout<<endl<<endl;}
double operator[](size_t i)const{return c_val[i];}
double& operator[](size_t i){return c_val[i];}
size_t re_dim()const{return dim;}
friend c operator+(const c&a,const c&b);
friend c operator-(const c&a);
friend c operator-(const c&a,const c&b);
friend c rep(const double&d,size_t n);
friend c operator+(const c&a,const double&d);
friend c operator-(const c&a,const double&d);
friend c operator+(const double&d,const c&a);
friend c operator-(const double&d,const c&a);
friend c operator*(const double&d,const c&a);
friend c operator*(const c&a,const double&d);
friend c operator/(const double&d,const c&a);
friend c operator/(const c&a,const double&d);
friend c operator*(const c&a,const c&b);
friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}
friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}
vector<double> const_re_c_val()const{return c_val;}
//在使用const_re_c_val()时可进行拷贝后使用
};
//类matrix声明
class matrix
{
private:
vector<c> m_val;
size_t rdim;
size_t cdim;
public:
matrix():m_val(vector<c>()),rdim(0),cdim(0){}
matrix(const vector<c>&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}
matrix(size_t m,size_t n,const double&d):m_val(vector<c>(m,vector<double>(n,d))),rdim(m),cdim(n){}
matrix(const c&v,size_t r,size_t c);
matrix(const c&v);
c operator[](size_t i)const{return m_val[i];}
c& operator[](size_t i){return m_val[i];}
void out(size_t n=3)const;
c col(size_t i)const;
void col_change(size_t i,const c&a);
friend matrix operator+(const matrix&a,const matrix&b);
friend matrix operator-(const matrix&a);
friend matrix operator-(const matrix&a,const matrix&b);
friend matrix operator+(const matrix&a,const double&d);
friend matrix operator+(const double&d,const matrix&a);
friend matrix operator-(const matrix&a,const double&d);
friend matrix operator-(const double&b,const matrix&a);
friend matrix operator*(const matrix&a,const double&d);
friend matrix operator*(const double&d,const matrix&a);
friend matrix operator/(const matrix&a,const double&d);
friend matrix operator/(const double&d,const matrix&a);
void Swap(size_t i,size_t j)
{c temp=this->operator[](i);this->operator[](i)=this->operator[](j);this->operator[](j)=temp;}
void col_Swap(size_t i,size_t j)
{
auto temp=this->col(i);
this->col_change(i,this->col(j));
this->col_change(j,temp);
}
friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}
friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}
pair<size_t,size_t> re_dim()const{return make_pair(rdim,cdim);}
friend c operator%(const matrix&m,const c&v);
friend matrix operator%(const matrix&a,const matrix&b);
friend matrix operator%(const c&v,const matrix&m);
};
//类factor声明
class factor
{
private:
vector<pair<double,size_t>>f_val;
map<double,size_t>f_val_to_levels;
vector<size_t>i_val;
set<double>Levels;
public:
factor(const c&v=c());
void out(size_t n=1);
friend c levels(const factor&f);
friend size_t nlevels(const factor&f);
friend factor as_integer(const factor&f);
friend factor as_factor(const c&v);
};
#endif // CLASS_OF_R_H_
类成员函数定义文件：
class_of_R.cpp
#include”class_of_R.h”
//类c的有关定义
c::c(initializer_list<double>lst)
{
vector<double>out;
for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});
c_val=out;
dim=out.size();
}
c::c(const matrix&m)
{
vector<double>temp;
for(size_t j=0;j<m.re_dim().second;j++)
for(size_t i=0;i<m.re_dim().first;i++)
temp.push_back(m[i][j]);
c_val=temp;
dim=temp.size();
}
//类matrix的有关定义
void matrix::out(size_t n)const
{
cout<<”The matrix is:\n”;
for(size_t i=0;i<rdim;i++)
{
for(size_t j=0;j<cdim;j++)
cout<<setw(n)<<(*this)[i][j]<<“ ”;
cout<<endl;
}
cout<<endl;
}
c matrix::col(size_t i)const
{
vector<double>temp;
for(size_t i=0;i<rdim;i++)
temp.push_back(0);
c out(temp);
for(size_t j=0;j<rdim;j++)
out[j]=(*this)[j][i];
return out;
}
void matrix::col_change(size_t i,const c&a)
{
for(size_t j=0;j<a.re_dim();j++)
(*this)[j][i]=a[j];
}
matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c)
{
matrix out(r,c,0);
for(size_t i=0;i<v.re_dim();i++)
out[i%r][i/r]=v[i];
(this->m_val)=out.m_val;
}
matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1)
{
matrix out(v.re_dim(),1,0);
for(size_t i=0;i<out.rdim;i++)
out[i][0]=v[i];
this->m_val=out.m_val;
}
//类factor有关定义
factor::factor(const c&v)
{
vector<double>temp=v.const_re_c_val();
for_each(temp.begin(),temp.end(),
[&](const double&d){Levels.insert(d);}
);
vector<double>s_to_v(Levels.begin(),Levels.end());
vector<pair<double,size_t>>ttemp;
for(size_t i=0;i<s_to_v.size();i++)
ttemp.push_back(make_pair(s_to_v[i],i));
for(size_t i=0;i<temp.size();i++)
{
for(size_t j=0;j<ttemp.size();j++)
{
if(temp[i]==ttemp[j].first)
{
f_val.push_back(ttemp[j]);
}
}
}
for(size_t i=0;i<ttemp.size();i++)
f_val_to_levels.insert(ttemp[i]);
auto a=f_val.begin();
for(size_t i=0;i<f_val.size();i++)
{
i_val.push_back(a->second);
a++;
}
}
void factor::out(size_t n)
{
vector<double>temp(f_val.size(),0);
c out0(temp);
auto a=f_val.begin();
for(size_t i=0;i<out0.re_dim();i++)
{
out0[i]=a->first;
a++;
}
vector<double>temp0(Levels.size(),0);
c out1(temp0);
set<double>::iterator b=Levels.begin();
for(size_t i=0;i<out1.re_dim();i++)
{
out1[i]=*b;
b++;
}
out0.out(n);
cout<<”Levels: ”;
cout<<endl;
out1.out(n);
}
其他函数即测试程序文件：
use_R.cpp
#include”class_of_R.h”
//若干友元函数
c operator+(const c&a,const c&b)
{
size_t MAX=b.re_dim();
if(a.re_dim()>b.re_dim())
MAX=a.re_dim();
vector<double>temp0(MAX,0);
vector<double>temp1(MAX,0);
vector<double>temp2(MAX,0);
copy(a.c_val.begin(),a.c_val.end(),temp0.begin());
copy(b.c_val.begin(),b.c_val.end(),temp1.begin());
for(size_t i=0;i<MAX;i++)
temp2[i]=temp0[i]+temp1[i];
return c(temp2);
}
c operator-(const c&a)
{
c out(a);
for(size_t i=0;i<out.re_dim();i++)
out[i]=-out[i];
return out;
}
c operator-(const c&a,const c&b)
{
return (a+(-b));
}
c rep(const double&d,size_t n)
{
return c(vector<double>(n,d));
}
c operator+(const c&a,const double&d)
{
return a+rep(d,a.re_dim());
}
c operator-(const c&a,const double&d)
{
return (a+(-d));
}
c operator+(const double&d,const c&a)
{
return a+d;
}
c operator-(const double&d,const c&a)
{
return (-a)+d;
}
c operator*(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]*=d;
return out;
}
c operator*(const c&a,const double&d)
{
return d*a;
}
c operator/(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]=d/a[i];
return out;
}
c operator/(const c&a,const double&d)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]=a[i]/d;
return out;
}
c operator*(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error(“dim error!\n”);
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]*=b[i];
return out;
}
matrix operator+(const matrix&a,const matrix&b)
{
size_t r_MAX=a.rdim,c_MAX=a.cdim;
if(a.rdim<b.rdim)
r_MAX=b.rdim;
if(a.cdim<b.cdim)
c_MAX=b.cdim;
matrix tempa(r_MAX,c_MAX,0);
for(size_t i=0;i<a.rdim;i++)
for(size_t j=0;j<a.cdim;j++)
tempa[i][j]=a[i][j];
matrix tempb(r_MAX,c_MAX,0);
for(size_t i=0;i<b.rdim;i++)
for(size_t j=0;j<b.cdim;j++)
tempb[i][j]=b[i][j];
matrix tempout(r_MAX,c_MAX,0);
for(size_t i=0;i<tempout.rdim;i++)
for(size_t j=0;j<tempout.cdim;j++)
tempout[i][j]=tempa[i][j]+tempb[i][j];
return tempout;
}
matrix operator-(const matrix&a)
{
matrix out(a);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]=-out[i][j];
return out;
}
matrix operator-(const matrix&a,const matrix&b)
{
return ((-b)+a);
}
matrix operator+(const matrix&a,const double&d)
{
return a+matrix(a.rdim,a.cdim,d);
}
matrix operator+(const double&d,const matrix&a)
{
return a+d;
}
matrix operator-(const matrix&a,const double&d)
{
return a+(-d);
}
matrix operator-(const double&b,const matrix&a)
{
return (-a)+b;
}
matrix operator*(const matrix&a,const double&d)
{
matrix out(a);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]*=d;
return out;
}
matrix operator*(const double&d,const matrix&a)
{
return a*d;
}
matrix operator/(const matrix&a,const double&d)
{
return a*(1/d);
}
matrix operator/(const double&d,const matrix&a)
{
matrix out(a.rdim,a.cdim,0);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]=d/a[i][j];
return out;
}
c levels(const factor&f)
{
vector<double> temp(f.Levels.begin(),f.Levels.end());
c done(temp);
done.out();
return done;
}
size_t nlevels(const factor&f)
{
size_t out=f.Levels.size();
cout<<out<<endl;
return out;
}
//返回输入因子相应在存储中的下标构成的因子
factor as_integer(const factor&f)
{
vector<double> temp(f.i_val.begin(),f.i_val.end());
factor temp0=factor(c(temp));
for_each(temp0.f_val.begin(),temp0.f_val.end(),
[](const pair<double,size_t>&p)
{
cout<<p.first<<” ”;
}
);
cout<<endl;
return temp0;
}
factor as_factor(const c&v)
{
factor temp0=factor(v);
temp0.out();
return temp0;
}
//基本运算函数
matrix diag(const c&v)
{
matrix out(v.re_dim(),v.re_dim(),0);
for(size_t i=0;i<v.re_dim();i++)
out[i][i]=v[i];
return out;
}
matrix diag(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
out[i][i]=m[i][i];
return out;
}
double inter_product(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error(“dim error!\n”);
double out=0;
for(size_t i=0;i<a.re_dim();i++)
out+=(a[i]*b[i]);
return out;
}
c operator%(const matrix&m,const c&v)
{
if(m.re_dim().second!=v.re_dim())
throw runtime_error(“dim error!\n”);
c out(rep(0,m.re_dim().first));
for(size_t i=0;i<out.re_dim();i++)
out[i]=inter_product(m[i],v);
return out;
}
//同下
double prod(const c&v)
{
double out=1;
for(size_t i=0;i<v.re_dim();i++)
out*=v[i];
return out;
}
//对于matrix的求和由c复制构造函数的重载函数定义：c(const matrix&m)
double sum(const c&v)
{
double out=v[0];
for(size_t i=1;i<v.re_dim();i++)
out+=v[i];
return out;
}
matrix operator%(const c&v,const matrix&m)
{
if(v.re_dim()!=m.re_dim().first)
throw runtime_error(“dim error!\n”);
matrix out(1,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().second;i++)
out[0][i]=inter_product(v,m.col(i));
return out;
}
matrix outer_product(const c&a,const c&b)
{
matrix out(a.re_dim(),b.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=a[i]*b[j];
return out;
}
size_t length(const c&v)
{
return v.re_dim();
}
c sort(const c&v)
{
vector<double>out;
for(size_t i=0;i<v.re_dim();i++)
out.push_back(v[i]);
sort(out.begin(),out.end());
return c(out);
}
double mean(const c&v)
{
return sum(v)/length(v);
}
//其转换版本返回矩阵转为向量的方差
double var(const c&v)
{
return mean((v-mean(v))*(v-mean(v)));
}
c connect(initializer_list<c> lst)
{
vector<double>temp;
for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});
return c(temp);
}
//标准差
template<typename any>
double STD(const any&a)
{
return sqrt(var(a)*(length(a))/(length(a)-1));
}
matrix rbind(initializer_list<c>lst)
{
size_t r=0;
auto a=lst.begin();
while(a!=lst.end())
{
a++;
r++;
}
size_t co=0;
for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});
matrix out(r,co,0);
size_t i=0;
a=lst.begin();
while(a!=lst.end())
{
vector<double>temp(co,0);
for(size_t i=0;i<a->const_re_c_val().size();i++)
temp[i]=a->operator[](i);
out[i]=c(temp);
i++;
a++;
}
return out;
}
size_t dim(const c&v)
{
return v.re_dim();
}
pair<size_t,size_t> dim(const matrix&m)
{
return m.re_dim();
}
double min(const c&v)
{
double temp=v[0];
for(size_t i=1;i<length(v);i++)
if(v[i]<temp)
temp=v[i];
return temp;
}
double max(const c&v)
{
double temp=v[0];
for(size_t i=1;i<length(v);i++)
if(v[i]>temp)
temp=v[i];
return temp;
}
c pmin(initializer_list<c>lst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i<temp.re_dim();i++)
{
temp[i]=min(*a);
a++;
}
return temp;
}
c pmax(initializer_list<c>lst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i<temp.re_dim();i++)
{
temp[i]=max(*a);
a++;
}
return temp;
}
template<typename any>
pair<double,double> range(const any&v)
{
cout<<”MIN: ”<<min(v)<<“ MAX: ”<<max(v)<<endl;
return make_pair(min(v),max(v));
}
double median(const c&v)
{
return sort(v)[length(v)/2];
}
double var(const c&a,const c&b)
{
if(length(a)!=length(b))
throw runtime_error(“dim error!\n”);
c t1(a-mean(a)),t2(b-mean(b));
return inter_product(t1,t2)/(length(a)-1);
}
double cov(const c&a,const c&b)
{
return var(a,b);
}
double cor(const c&a,const c&b)
{
return var(a,b)/STD(a)/STD(b);
}
matrix var(const matrix&a,const matrix&b)
{
if(a.re_dim().first!=b.re_dim().first)
throw runtime_error(“dim error!\n”);
matrix temp(a.re_dim().second,b.re_dim().second,0);
for(size_t i=0;i<temp.re_dim().first;i++)
for(size_t j=0;j<temp.re_dim().second;j++)
temp[i][j]=var(a.col(i),b.col(j));
return temp;
}
matrix cov(const matrix&a,const matrix&b)
{
return var(a,b);
}
matrix cor(const matrix&a,const matrix&b)
{
matrix temp=var(a,b);
for(size_t i=0;i<temp.re_dim().first;i++)
for(size_t j=0;j<temp.re_dim().second;j++)
temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));
return temp;
}
matrix cov(const matrix&a)
{
return var(a,a);
}
matrix cor(const matrix&a)
{
return cor(a,a);
}
//由于关联容器的删除方法利用误差排序的方法
//利用顺序容器迭代器的删除方法有困难
c rank(const c&v)
{
vector<double>temp00(v.const_re_c_val());
for(size_t i=0;i<temp00.size();i++)
temp00[i]+=(0.00000001*i);
vector<double>temp0(temp00);
multiset<double>temp2(temp0.begin(),temp0.end());
vector<double>temp1(v.re_dim());
size_t i=0,j=0;
double pre_val=min(c(vector<double>(temp2.begin(),temp2.end())));
vector<double>::iterator s0;
while(i<v.re_dim())
{
double val=min(c(vector<double>(temp2.begin(),temp2.end())));
s0=find(temp0.begin(),temp0.end(),val);
if(sqrt((pre_val-val)*(pre_val-val))>0.000001)
j++;
temp1[s0-temp0.begin()]=j;
temp2.erase(val);
i++;
pre_val=val;
}
return c(temp1);
}
c rev(const c&v)
{
vector<double> temp;
auto a=v.const_re_c_val().rbegin();
while(a!=v.const_re_c_val().rend())
{
temp.push_back(*a);
a++;
}
return c(temp);
}
//向量先排序再输出排序前的下标值
c order(const c&v)
{
unordered_multimap<double,double>m0;
for(size_t i=0;i<length(v);i++)
m0.insert(make_pair(v[i],i));
multimap<double,double>m1(m0.begin(),m0.end());
c temp(rep(0,length(v)));
auto a=m1.begin();
size_t i=0;
while(a!=m1.end())
{
temp[i]=a->second;
i++;
a++;
}
return temp;
}
c cumsum(const c&v)
{
c temp(v);
for(size_t i=1;i<dim(v);i++)
temp[i]+=temp[i-1];
return temp;
}
c cumprod(const c&v)
{
c temp(v);
for(size_t i=1;i<length(v);i++)
temp[i]*=temp[i-1];
return temp;
}
matrix cumsum(const matrix&m)
{
c temp(m);
for(size_t i=1;i<length(m);i++)
temp[i]+=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
matrix cumprod(const matrix&m)
{
c temp(m);
for(size_t i=1;i<length(m);i++)
temp[i]*=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
//矩阵运算函数
//初等行变换化上三角求解行列式
double det(const matrix&m)
{
auto f=[](const matrix&mm)->pair<matrix,double>
{
if(mm.re_dim().first!=mm.re_dim().second)
throw runtime_error(“dim error!\n”);
size_t MIN=mm.re_dim().first;
double mult=1;
size_t col=0;
size_t bar=0;
size_t b_c=0;
matrix t(mm);
while(b_c<MIN)
{
for(size_t i=bar+1;i<MIN;)
{
if(t[i][col]==0)
i++;
else
{
t.Swap(i,bar);
mult*=(-1);
break;
}
}
if(bar!=MIN+1)
{
if(t[bar][col]!=0)
{
for(size_t i=bar+1;i<MIN;i++)
{
t[bar][i]=t[bar][i]/t[bar][col];
if(i==bar+1)
mult*=t[bar][col];
}
if(col!=MIN-1)
t[bar][col]=1;
}
for(size_t i=bar+1;i<MIN;i++)
{
t[i]=t[i]-t[i][col]*t[bar];
}
}
col++;
bar++;
b_c=bar;
if(col>bar);
b_c=col;
}
return make_pair(t,mult);
};
pair<matrix,double> tt=f(m);
return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));
}
//Cramer法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0
c solve(const matrix&m,const c&v)
{
c temp(rep(0,v.re_dim()));
double low=det(m);
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
if(low==0)
throw runtime_error(“singular!\n”);
for(size_t i=0;i<v.re_dim();i++)
{
matrix temp0(m);
temp0.col_change(i,v);
temp[i]=det(temp0)/low;
}
return temp;
}
double tr(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));
}
matrix T(const matrix&m)
{
matrix out(m.re_dim().second,m.re_dim().first,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=m[j][i];
return out;
}
matrix cbind(initializer_list<c>lst)
{
return T(rbind(lst));
}
matrix operator%(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first,b.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=inter_product(a[i],b.col(j));
return out;
}
double F_norm(const matrix&m)
{
return sqrt(tr(T(m)%m));
}
//a,b为区间 d为步长 e为精度 默认为求整数精度：以1为精度特征值 如所求特征值太少可通过减少d改善
//只能求实值
//可以用F_norm作为a、b的界（-F_norm(m),F_norm(m)）可确保值在范围中
//可用于求解但效率低
c eigen(const matrix&m,const double&a,const double &b,const double &d=0.0001,const double &e=1)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
vector<double>temp;
double pre_beg=100000000;
for(double beg=a;beg<b;beg+=d)
{
matrix mm0(diag(rep(beg,m.re_dim().first))-m);
double s0=det(mm0);
if(sqrt(s0*s0)<e)
{
if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)
{
temp.push_back(beg);
pre_beg=beg;
cout<<”Done!\n”;
cout<<beg<<endl;//eigenvalue
}
}
}
return c(temp);
}
//重载求eign的函数 d:步长 可保证绝对求出所有特征值（利用循环增大精度）
c eigen(const matrix&m,const double&d=0.1)
{
double dd=d;
double ra=F_norm(m);
double trr=tr(m);
c temp;
double r=0;
do
{
dd*=0.1;//原来设为0.1
temp=eigen(m,-ra,ra,dd);
if(temp.re_dim()==m.re_dim().first)
break;
r=sum(temp)-trr;
}
while(sqrt((r*r)>1));
//与迹的差距（精度）设定为1 绝对值小于此值的特征值被忽略
//此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握
//只适用于所有特征值为实值的情况
//如有复特征值，只输出实数的并循环不能结束
return temp;
}
//矩阵的余子阵
matrix M(const matrix&m,size_t k,size_t l)
{
matrix temp(m);
set<size_t> s1;
for(size_t i=0;i<m.re_dim().first;i++)
s1.insert(i);
set<size_t> s2;
for(size_t i=0;i<m.re_dim().second;i++)
s2.insert(i);
s1.erase(k);
s2.erase(l);
vector<double>temp00;
for(auto j:s2)
for(auto i:s1)
temp00.push_back(temp[i][j]);
return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);
}
//方阵的伴随矩阵
matrix A_accompany(const matrix&m)
{
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);
//代数余子式
return T(out);
}
//求方阵的逆矩阵（利用伴随矩阵方法）
matrix solve(const matrix&m)
{
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error(“singular!\n”);
return A_accompany(m)/de;
}
//利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆
c solve_eq(const matrix&m,const c&v)
{
c out;
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error(“singular!\n”);
return solve(m)%v;
}
//化列向量组构成的矩阵为正交化后列向量组构成的矩阵（Schmidt 正交化）
matrix orth_normal(const matrix &m)
{
matrix out(rep(0,m.re_dim().first*m.re_dim().second),m.re_dim().first,m.re_dim().second);
vector<c>dec;
dec.push_back(m.col(0)/sqrt(inter_product(m.col(0),m.col(0))));
for(size_t i=1;i<m.re_dim().second;i++)
{
c temp=m.col(i);
for(size_t j=1;j<=i;j++)
temp=temp-inter_product(dec[j-1],m.col(i))*dec[j-1];
temp=temp/sqrt(inter_product(temp,temp));
dec.push_back(temp);
}
for(size_t i=0;i<m.re_dim().second;i++)
out.col_change(i,dec[i]);
return out;
}
//生成相应位数的由0/1字符构成的字符串
vector<string> generator0_1(size_t n)
{
list<string>out;
list<string>require;
out.push_back(”0”);
out.push_back(”1”);
size_t i=0;
while(i<n-1)
{
if(i>0)
out=require;
size_t s0=out.size();
list<string>::iterator pre_begin=out.begin();
auto pre_end=out.end();
size_t k=0;
while(pre_begin!=pre_end)
{
pre_begin++;
k++;
}
pre_begin=out.begin();
for(size_t i=0;i<k;i++)
{
string temp0=*pre_begin+”0”;
string temp1=*pre_begin+”1”;
out.push_back(temp0);
out.push_back(temp1);
pre_begin++;
}
require.clear();
list<string>::iterator a=out.begin();
for(size_t i=0;i<s0;i++)
a++;
while(a!=out.end())
{
require.push_back(*a);
a++;
}
i++;
}
return vector<string>(require.begin(),require.end());
}
//利用矩阵向量方法生成上述数组
matrix matrix_generator0_1(size_t n)
{
matrix temp0(2,1,0);
matrix temp1;
temp0[0][0]=0;
temp0[1][0]=1;
size_t i=0;
while(i<n-1)
{
if(i>0)
temp0=temp1;
// size_t s0=temp0.re_dim().second;
size_t pre_begin=0;
auto pre_end=temp0.re_dim().first;
temp1=matrix(temp0.re_dim().first*2,temp0.re_dim().second+1,0);
size_t temp1_index=0;
while(pre_begin!=pre_end)
{
vector<double> tt0=temp0[pre_begin].const_re_c_val();
vector<double> tt1=tt0;
tt0.push_back(0);
tt1.push_back(1);
temp1[temp1_index]=c(tt0);
temp1_index++;
temp1[temp1_index]=c(tt1);
temp1_index++;
pre_begin++;
}
i++;
}
return temp1;
}
//将固定位数的下标值与0/1数值构成的map以vector的形式返回
//每个map都是一种组合
vector<map<size_t,size_t>> index_map_generator(size_t n)
{
vector<string>s0(generator0_1(n));
vector<map<size_t,size_t>>outer;
for(size_t i=0;i<s0.size();i++)
{
map<size_t,size_t>tp;
for(size_t j=0;j<s0[i].size();j++)
{
size_t temp=1;
if(s0[i][j]==‘0’)
temp=0;
tp.insert(make_pair(j,temp));
}
outer.push_back(tp);
}
return outer;
}
//上述产生方式的矩阵版本
vector<map<size_t,size_t>> index_map_matrix_generator(size_t n)
{
matrix temp=matrix_generator0_1(n);
vector<map<size_t,size_t>> out(temp.re_dim().first);
for(size_t i=0;i<out.size();i++)
{
map<size_t,size_t>m;
for(size_t j=0;j<temp.re_dim().second;j++)
m.insert(make_pair(j,temp[i][j]));
out[i]=m;
}
return out;
}
//返回上述index_map_generator(m)中指向0_1和为n的子集
vector<map<size_t,size_t>> require_dim(size_t m,size_t n)
{
vector<map<size_t,size_t>>v0=index_map_generator(m);
vector<map<size_t,size_t>>out;
auto a=v0.begin();
while(a!=v0.end())
{
size_t temp=0;
for_each(a->begin(),a->end(),[&](const pair<size_t,size_t>&p){
if(p.second==1)
temp++;
});
if(temp==n)
out.push_back(*a);
a++;
}
return out;
}
//matrix 所有n阶子式构成的向量
vector<double> seq_of_subdet(const matrix&m,size_t n)
{
vector<double>out;
vector<map<size_t,size_t>>r_map=require_dim(m.re_dim().first,n);&
4000
nbsp;
vector<map<size_t,size_t>>c_map=require_dim(m.re_dim().second,n);
for(size_t i=0;i<r_map.size();i++)
for(size_t j=0;j<c_map.size();j++)
{
set<size_t>r_require;
for_each(r_map[i].begin(),r_map[i].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)r_require.insert(p.first);});
set<size_t>c_require;
for_each(c_map[j].begin(),c_map[j].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)c_require.insert(p.first);});
matrix r_matrix(n,m.re_dim().second,0);
auto a=r_require.begin();
for(size_t i=0;i<r_matrix.re_dim().first;i++)
{
r_matrix[i]=m[*a];
a++;
}
matrix c_matrix(r_matrix.re_dim().first,n,0);
auto b=c_require.begin();
for(size_t i=0;i<n;i++)
{
c_matrix.col_change(i,r_matrix.col(*b));
b++;
}
out.push_back(det(c_matrix));
}
return out;
}
//利用子式的方法求出矩阵的秩
size_t rank(const matrix&m)
{
size_t MAX=m.re_dim().first;
if(m.re_dim().second<m.re_dim().first)
MAX=m.re_dim().second;
size_t out=MAX;
size_t pre_num_of_non0=0;
for(size_t i=MAX;i>0;)
{
vector<double>temp=seq_of_subdet(m,i);
auto a=temp.begin();
while(a!=temp.end())
{
*a=sqrt((*a)*(*a));
a++;
}
size_t num_of_non0=0;
for_each(temp.begin(),temp.end(),[&](const double&d){if(d>0.00000001)num_of_non0++;});
if(pre_num_of_non0==0)
i–;
if(num_of_non0!=0)
return i+1;
pre_num_of_non0=num_of_non0;
}
return out;
}
//回归系数的求解 要求增广后的数据阵（含常系数）可逆
c lm(const c&y,const matrix&m)
{
matrix X=matrix(m.re_dim().first,m.re_dim().second+1,0);
X.col_change(0,rep(1,X.re_dim().first));
for(size_t i=0;i<m.re_dim().second;i++)
X.col_change(i+1,m.col(i));
if(det(T(X)%X)<0.00000001)
throw runtime_error(“singular!\n”);
c out(solve(T(X)%X)%T(X)%y);
cout<<”The intercept is: ”<<out[0]<<endl;
cout<<”Others are:\n”;
for(size_t i=1;i<out.re_dim();i++)
cout<<”x”<<i<<“ is ”<<out[i]<<endl;
cout<<endl;
return out;
}
//抽样函数 参数为向量v 抽取的个数n 是否有放回replace 抽取的概率prob
//使用此方法前必须在函数调用外设定随机数种子
//要求没有0概率点
c sample(const c&v,size_t n,bool replace=1,const c&prob=c())
{
if(!replace&&v.re_dim()<n)
throw out_of_range(“n is too large\n”);
//抽样概率
vector<double>probb;
if(prob==c())
probb=rep(1/(double)v.re_dim(),v.re_dim()).const_re_c_val();
else
{
probb=prob.const_re_c_val();
double total=0;
for(auto a:probb)
total+=a;
for(auto &a:probb)
a/=total;
}
ifstream fin(”h://runif0.txt”);
vector<double>runif;
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
//c(probb).out();
c cprobb=cumsum(c(probb));
//cprobb.out();
c require(rep(0,n));
if(replace)
{
size_t i=0;
while(i<n)
{
double uniform=runif[rand()%30000];
//cout<<uniform<<endl;
double conclusion=v[v.re_dim()-1];
for(size_t i=0;cprobb[i]<1;i++)
{
if(uniform<cprobb[i])
{
&n
1c143
bsp; conclusion=v[i];
break;
}
}
require[i]=conclusion;
i++;
}
}
else
{
size_t i=0;
map<double,double>mdd;
for(size_t i=0;i<v.re_dim();i++)
mdd.insert(make_pair(v[i],c(probb)[i]));
while(i<n)
{
vector<double> temp_v;
vector<double> temp_prob;
for_each(mdd.begin(),mdd.end(),[&](const pair<double,double>&p){
temp_v.push_back(p.first);
temp_prob.push_back(p.second);
});
double d=sample(c(temp_v),1,1,c(temp_prob)).const_re_c_val()[0];
require[i]=d;
mdd.erase(d);
i++;
}
}
return require;
}
//产生几何分布随机数 n 产生个数 p成功概率
c rgeom(const size_t& n,const double& p)
{
vector<double>out;
size_t i=0;
srand(0);
ifstream fin(”h://runif0.txt”);
vector<double>runif;
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
while(i<n)
{
double temp=runif[rand()%30000];
out.push_back(size_t(log(temp)/log(1-p))+1);
i++;
}
return c(out);
}
//possion分布随机数 n 生成个数 d 分布参数
c rpois(const size_t&n,const double&d)
{
size_t i=0;
vector<double>conclusion;
srand(time(0));
ifstream fin(”h://runif0.txt”);
vector<double>runif;//储存30000个均匀随机数的vector
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
while(i<n)
{
double p=exp(-1*d);
size_t out=0;
double F=p;
double temp=runif[rand()%30000];
while(temp>=F)
{
p=d/(double)(out+1)*p;
F+=p;
out++;
}
conclusion.emplace_back(out);
i++;
}
//for_each(conclusion.begin(),conclusion.end(),[](const double&s){cout<<s<<” ”;});
//cout<<endl;
return c(conclusion);
}
//利用间距生成向量
c seq(const double&from,const double&to,const double&by)
{
if(from>to)
throw runtime_error(“from > to\n”);
vector<double>out;
auto begin=from;
while(begin<=to)
{
out.push_back(begin);
begin+=by;
}
return c(out);
}
//生成列表向量元素的所有组合并将结果以矩阵形式输出
//前面的generate_grid_ 函数为铺垫而设 expand_grid函数接受c为参数给出输出
matrix generate_grid_0(pair<double,c>p0)
{
matrix out(p0.second.re_dim(),2,0);
for(size_t i=0;i<out.re_dim().first;i++)
{
out[i][0]=p0.first;
out[i][1]=p0.second[i];
}
return out;
}
matrix generate_grid_1(const c&c0,const c&c1)
{
vector<double> temp0=c0.const_re_c_val();
size_t len=temp0.size();
matrix out(len*c1.re_dim(),2,0);
for(size_t i=0;i<len;i++)
{
matrix temp00=generate_grid_0(make_pair(temp0[i],c1));
for(size_t j=0;j<c1.re_dim();j++)
{
out[i*c1.re_dim()+j]=temp00[j];
}
}
return out;
}
matrix generate_grid_2(const matrix&m,const c&v)
{
matrix temp=generate_grid_1(m.col(m.re_dim().second-1),v);
//temp.out();
//system(“pause”);
matrix out(temp.re_dim().first,m.re_dim().second+1,0);
//out.out();
//system(“pause”);
matrix add(m.re_dim().first,m.re_dim().second-1,0);
for(size_t i=0;i<add.re_dim().second;i++)
add.col_change(i,m.col(i));
//add.out();
//system(“pause”);
for(size_t i=0;i<m.re_dim().first;i++)
for(size_t j=0;j<v.re_dim();j++)
{
for(size_t k=0;k<add.re_dim().second;k++)
out[i*v.re_dim()+j][k]=add[i][k];
}
out.col_change(out.re_dim().second-2,temp.col(0));
out.col_change(out.re_dim().second-1,temp.col(1));
return out;
}
//若干向量其元素的所有组合
matrix expand_grid(initializer_list<c>i0)
{
vector<c>temp(i0);
if(temp.size()==1)
return matrix(temp[0]);
matrix beg0=generate_grid_1(temp[0],temp[1]);
if(temp.size()==2)
return beg0;
for(size_t i=2;i<temp.size();i++)
{
//beg0.out();
beg0=generate_grid_2(beg0,temp[i]);
//cout<<”Done!\n”;
//beg0.out();
//cout<<”Done!\n”;
}
return beg0;
}
//得到向量v的所有n个元素的组合 输出以矩阵排列 每列为一个组合
//调用了求得下标组合的 require_dim()
//当组合项数太多时主张利用转置按行查看
matrix combn(const c&v,size_t n)
{
auto a=require_dim(v.re_dim(),n);
auto b=vector<set<size_t>>(a.size());
auto a0=a.begin();
auto b0=b.begin();
while(a0!=a.end())
{
for_each(a0->begin(),a0->end(),
[&](const pair<size_t,size_t>&p)
{
if(p.second==1)
b0->insert(p.first);
});
a0++;
b0++;
}
auto cc=vector<vector<double>>(b.size());
for(size_t i=0;i<cc.size();i++)
{
cc[i]=vector<double>(b[i].begin(),b[i].end());
}
matrix out(cc[0].size(),cc.size(),0);
for(size_t i=0;i<out.re_dim().second;i++)
{
out.col_change(i,cc[i]);
}
matrix oout(out.re_dim().first,out.re_dim().second,0);
for(size_t i=0;i<oout.re_dim().first;i++)
for(size_t j=0;j<oout.re_dim().second;j++)
oout[i][j]=v[out[i][j]];
return oout;
}
//上述函数的重载版本 返回0:m的相应n组合
matrix combn(size_t m,size_t n)
{
vector<double>temp;
for(size_t i=0;i<m+1;i++)
temp.push_back(i);
return combn(temp,n);
}
//向量相应元素在前面是否出现
//出现显示1
c duplicated(const c&v)
{
//函数a返回向量cc的前n个元素构成的向量
auto a=[](const c&cc,size_t n)->c
{
c out(rep(0,n));
for(size_t i=0;i<n;i++)
out[i]=cc[i];
return out;
};
vector<c>vc0(v.re_dim());
for(size_t i=0;i<vc0.size();i++)
vc0[i]=a(v,i+1);
vector<vector<double>>vs0(vc0.size());
for(size_t i=0;i<vs0.size();i++)
vs0[i]=vector<double>(vc0[i].const_re_c_val());
c out(rep(0,v.re_dim()));
for(size_t i=1;i<v.re_dim();i++)
{
if(find(vs0[i-1].begin(),vs0[i-1].end(),v[i])!=vs0[i-1].end())
out[i]=1;
}
return out;
}
//实现unique的操作看似与factor的levels相同
//实际在序上是不同的
//levels保持double的序（由set有序决定）
//unique应保持原向量的序
//unique由duplicated定义
c unique(const c&v)
{
c d0=duplicated(v);
vector<pair<double,double>>vm0;
for(size_t i=0;i<v.re_dim();i++)
vm0.push_back(make_pair(v[i],d0[i]));
vector<double>out;
for(size_t i=0;i<vm0.size();i++)
{
if(vm0[i].second==0)
out.push_back(vm0[i].first);
}
return c(out);
}
//矩阵叉乘 T(a)%b (向量内积的矩阵版本)
matrix crossprod(const matrix&a,const matrix&b)
{
return T(a)%b;
}
//outer 生成外积 相应的可将外机作用于函数
matrix outer(const c&v1,const c&v2,function<double(const double&,const double&)>p=function<double(const double&,const double&)>())
{
if(!p)
{
matrix out(v1.re_dim(),v2.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=v1[i]*v2[j];
return out;
}
else
{
matrix out(v1.re_dim(),v2.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=p(v1[i],v2[j]);
return out;
}
}
//矩阵Kronecker积
matrix kronecker(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first*b.re_dim().first,a.re_dim().second*b.re_dim().second,0);
vector<vector<matrix>>vvm0(a.re_dim().first,vector<matrix>(a.re_dim().second));
for(size_t i=0;i<a.re_dim().first;i++)
for(size_t j=0;j<a.re_dim().second;j++)
{
vvm0[i][j]=a[i][j]*b;
}
for(size_t i=0;i<vvm0.size();i++)
for(size_t j=0;j<vvm0[0].size();j++)
{
matrix temp=vvm0[i][j];
for(size_t k=0;k<temp.re_dim().first;k++)
for(size_t l=0;l<temp.re_dim().second;l++)
out[i*b.re_dim().first+k][j*b.re_dim().second+l]=temp[k][l];
}
return out;
}
//由方阵上下三角部分组成的0-1阵 不包含对角线元素
matrix lower_tri(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=1;i<out.re_dim().first;i++)
for(size_t j=0;j<i;j++)
out[i][j]=1;
return out;
}
matrix upper_tri(const matrix&m)
{
return T(lower_tri(m));
}
//矩阵拉直(按列)
c vec(const matrix&m)
{
c out=m.col(0);
for(size_t i=1;i<m.re_dim().second;i++)
out=connect({out,m.col(i)});
return out;
}
//矩阵半拉直 拉直运算忽略主对角线以上部分
c vech(const matrix&m)
{
c out=m.col(0);
for(size_t i=1;i<m.re_dim().second;i++)
{
vector<double> temp=m.col(i).const_re_c_val();
auto a=temp.begin();
for(size_t j=0;j<i;j++)
a++;
vector<double>temp0(a,temp.end());
out=connect({out,c(temp0)});
}
return out;
}
//在方阵满秩情况下 给出方阵A的QR分解
pair<matrix,matrix> qr_of_full_rank(const matrix&A)
{
if(A.re_dim().first!=A.re_dim().second)
throw runtime_error(“dim error!\n”);
// /*
size_t rk=rank(A);
cout<<”The rank: ”<<rk<<endl;
if(rk!=A.re_dim().first)
throw runtime_error(“singular!\n”);
//*/
cout<<”A=QR”<<endl;
cout<<endl;
cout<<”Q of A:\n”;
matrix Q=orth_normal(A);
Q.out();
cout<<”R of m0:\n”;
matrix R=T(Q)%A;
R.val_to_double(lower_tri(R),0);
R.out();
// /*
cout<<”The F_norm gap between A and QR:\n”;
cout<<F_norm(A-Q%R)<<endl;
//*/
return make_pair(Q,R);
}
//要求自变量矩阵为可逆方阵 利用qr分解 求解线性方程组
//这只是一个形式上的实现 因为求解R逆的方法是用伴随矩阵法
//如能改为初等行变换法应会提高效率
c qr_solve(const matrix&m,const c&v)
{
auto qr=qr_of_full_rank(m);
return solve(qr.second)%T(qr.first)%v;
}
int main()
{
//测试程序：
matrix m2(c({5,0,10,11,0,4,4,3,10,4,3,3,11,3,3,2,10,10,10,10,11,12,13,14,15}),5,5);
m2.out();
c v1({1,1,1,1,2});
c v2=solve(m2,v1);
c v3=qr_solve(m2,v1);
v2.out();
v3.out();
cout<<”The gap of F_norm:\n”;
cout<<F_norm(v2-v3)<<endl;
return 0;
}
using std::endl;
using std::cin;
using std::istream;
using std::ios_base;
using std::sort;
using std::copy;
using std::find;
using std::out_of_range;
using std::rand;
using std::system;
using std::log;
using std::exp;
using std::make_pair;
using std::set;
using std::unordered_multimap;
using std::map;
class matrix;
//类c声明
class c
{
private:
vector<double> c_val;
size_t dim;
public:
c():c_val(vector<double>()),dim(0){}
c(initializer_list<double>lst);
c(const vector<double>&c_v):c_val(c_v),dim(c_v.size()){}
c(const matrix&m);
void out(size_t n=1)const{cout<<"The vector is:\n";for(size_t i=0;i<dim;i++)cout<<setw(n)<<c_val[i]<<" ";cout<<endl<<endl;}
double operator[](size_t i)const{return c_val[i];}
double& operator[](size_t i){return c_val[i];}
size_t re_dim()const{return dim;}
friend c operator+(const c&a,const c&b);
friend c operator-(const c&a);
friend c operator-(const c&a,const c&b);
friend c rep(const double&d,size_t n);
friend c operator+(const c&a,const double&d);
friend c operator-(const c&a,const double&d);
friend c operator+(const double&d,const c&a);
friend c operator-(const double&d,const c&a);
friend c operator*(const double&d,const c&a);
friend c operator*(const c&a,const double&d);
friend c operator/(const double&d,const c&a);
friend c operator/(const c&a,const double&d);
friend c operator*(const c&a,const c&b);
friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}
friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}
vector<double> const_re_c_val()const{return c_val;}
//在使用const_re_c_val()时可进行拷贝后使用
};
//类matrix声明
class matrix
{
private:
vector<c> m_val;
size_t rdim;
size_t cdim;
public:
matrix():m_val(vector<c>()),rdim(0),cdim(0){}
matrix(const vector<c>&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}
matrix(size_t m,size_t n,const double&d):m_val(vector<c>(m,vector<double>(n,d))),rdim(m),cdim(n){}
matrix(const c&v,size_t r,size_t c);
matrix(const c&v);
c operator[](size_t i)const{return m_val[i];}
c& operator[](size_t i){return m_val[i];}
void out(size_t n=3)const;
c col(size_t i)const;
void col_change(size_t i,const c&a);
friend matrix operator+(const matrix&a,const matrix&b);
friend matrix operator-(const matrix&a);
friend matrix operator-(const matrix&a,const matrix&b);
friend matrix operator+(const matrix&a,const double&d);
friend matrix operator+(const double&d,const matrix&a);
friend matrix operator-(const matrix&a,const double&d);
friend matrix operator-(const double&b,const matrix&a);
friend matrix operator*(const matrix&a,const double&d);
friend matrix operator*(const double&d,const matrix&a);
friend matrix operator/(const matrix&a,const double&d);
friend matrix operator/(const double&d,const matrix&a);
void Swap(size_t i,size_t j)
{c temp=this->operator;this->operator=this->operator;this->operator=temp;}
void col_Swap(size_t i,size_t j)
{
auto temp=this->col(i);
this->col_change(i,this->col(j));
this->col_change(j,temp);
}
friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}
friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}
pair<size_t,size_t> re_dim()const{return make_pair(rdim,cdim);}
friend c operator%(const matrix&m,const c&v);
friend matrix operator%(const matrix&a,const matrix&b);
friend matrix operator%(const c&v,const matrix&m);
};
//类factor声明
class factor
{
private:
vector<pair<double,size_t>>f_val;
map<double,size_t>f_val_to_levels;
vector<size_t>i_val;
set<double>Levels;
public:
factor(const c&v=c());
void out(size_t n=1);
friend c levels(const factor&f);
friend size_t nlevels(const factor&f);
friend factor as_integer(const factor&f);
friend factor as_factor(const c&v);
};
class_of_R.cpp
c::c(initializer_list<double>lst)
{
vector<double>out;
for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});
c_val=out;
dim=out.size();
}
c::c(const matrix&m)
{
vector<double>temp;
for(size_t j=0;j<m.re_dim().second;j++)
for(size_t i=0;i<m.re_dim().first;i++)
temp.push_back(m[i][j]);
c_val=temp;
dim=temp.size();
}
//类matrix的有关定义
void matrix::out(size_t n)const
{
cout<<"The matrix is:\n";
for(size_t i=0;i<rdim;i++)
{
for(size_t j=0;j<cdim;j++)
cout<<setw(n)<<(*this)[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
c matrix::col(size_t i)const
{
vector<double>temp;
for(size_t i=0;i<rdim;i++)
temp.push_back(0);
c out(temp);
for(size_t j=0;j<rdim;j++)
out[j]=(*this)[j][i];
return out;
}
void matrix::col_change(size_t i,const c&a)
{
for(size_t j=0;j<a.re_dim();j++)
(*this)[j][i]=a[j];
}
matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c)
{
matrix out(r,c,0);
for(size_t i=0;i<v.re_dim();i++)
out[i%r][i/r]=v[i];
(this->m_val)=out.m_val;
}
matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1)
{
matrix out(v.re_dim(),1,0);
for(size_t i=0;i<out.rdim;i++)
out[i][0]=v[i];
this->m_val=out.m_val;
}
//类factor有关定义
factor::factor(const c&v)
{
vector<double>temp=v.const_re_c_val();
for_each(temp.begin(),temp.end(),
[&](const double&d){Levels.insert(d);}
);
vector<double>s_to_v(Levels.begin(),Levels.end());
vector<pair<double,size_t>>ttemp;
for(size_t i=0;i<s_to_v.size();i++)
ttemp.push_back(make_pair(s_to_v[i],i));
for(size_t i=0;i<temp.size();i++)
{
for(size_t j=0;j<ttemp.size();j++)
{
if(temp[i]==ttemp[j].first)
{
f_val.push_back(ttemp[j]);
}
}
}
for(size_t i=0;i<ttemp.size();i++)
f_val_to_levels.insert(ttemp[i]);
auto a=f_val.begin();
for(size_t i=0;i<f_val.size();i++)
{
i_val.push_back(a->second);
a++;
}
}
void factor::out(size_t n)
{
vector<double>temp(f_val.size(),0);
c out0(temp);
auto a=f_val.begin();
for(size_t i=0;i<out0.re_dim();i++)
{
out0[i]=a->first;
a++;
}
vector<double>temp0(Levels.size(),0);
c out1(temp0);
set<double>::iterator b=Levels.begin();
for(size_t i=0;i<out1.re_dim();i++)
{
out1[i]=*b;
b++;
}
out0.out(n);
cout<<"Levels: ";
cout<<endl;
out1.out(n);
}
其他函数即测试程序文件：
use_R.cpp
c operator+(const c&a,const c&b)
{
size_t MAX=b.re_dim();
if(a.re_dim()>b.re_dim())
MAX=a.re_dim();
vector<double>temp0(MAX,0);
vector<double>temp1(MAX,0);
vector<double>temp2(MAX,0);
copy(a.c_val.begin(),a.c_val.end(),temp0.begin());
copy(b.c_val.begin(),b.c_val.end(),temp1.begin());
for(size_t i=0;i<MAX;i++)
temp2[i]=temp0[i]+temp1[i];
return c(temp2);
}
c operator-(const c&a)
{
c out(a);
for(size_t i=0;i<out.re_dim();i++)
out[i]=-out[i];
return out;
}
c operator-(const c&a,const c&b)
{
return (a+(-b));
}
c rep(const double&d,size_t n)
{
return c(vector<double>(n,d));
}
c operator+(const c&a,const double&d)
{
return a+rep(d,a.re_dim());
}
c operator-(const c&a,const double&d)
{
return (a+(-d));
}
c operator+(const double&d,const c&a)
{
return a+d;
}
c operator-(const double&d,const c&a)
{
return (-a)+d;
}
c operator*(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]*=d;
return out;
}
c operator*(const c&a,const double&d)
{
return d*a;
}
c operator/(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]=d/a[i];
return out;
}
c operator/(const c&a,const double&d)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]=a[i]/d;
return out;
}
c operator*(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error("dim error!\n");
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]*=b[i];
return out;
}
matrix operator+(const matrix&a,const matrix&b)
{
size_t r_MAX=a.rdim,c_MAX=a.cdim;
if(a.rdim<b.rdim)
r_MAX=b.rdim;
if(a.cdim<b.cdim)
c_MAX=b.cdim;
matrix tempa(r_MAX,c_MAX,0);
for(size_t i=0;i<a.rdim;i++)
for(size_t j=0;j<a.cdim;j++)
tempa[i][j]=a[i][j];
matrix tempb(r_MAX,c_MAX,0);
for(size_t i=0;i<b.rdim;i++)
for(size_t j=0;j<b.cdim;j++)
tempb[i][j]=b[i][j];
matrix tempout(r_MAX,c_MAX,0);
for(size_t i=0;i<tempout.rdim;i++)
for(size_t j=0;j<tempout.cdim;j++)
tempout[i][j]=tempa[i][j]+tempb[i][j];
return tempout;
}
matrix operator-(const matrix&a)
{
matrix out(a);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]=-out[i][j];
return out;
}
matrix operator-(const matrix&a,const matrix&b)
{
return ((-b)+a);
}
matrix operator+(const matrix&a,const double&d)
{
return a+matrix(a.rdim,a.cdim,d);
}
matrix operator+(const double&d,const matrix&a)
{
return a+d;
}
matrix operator-(const matrix&a,const double&d)
{
return a+(-d);
}
matrix operator-(const double&b,const matrix&a)
{
return (-a)+b;
}
matrix operator*(const matrix&a,const double&d)
{
matrix out(a);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]*=d;
return out;
}
matrix operator*(const double&d,const matrix&a)
{
return a*d;
}
matrix operator/(const matrix&a,const double&d)
{
return a*(1/d);
}
matrix operator/(const double&d,const matrix&a)
{
matrix out(a.rdim,a.cdim,0);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]=d/a[i][j];
return out;
}
c levels(const factor&f)
{
vector<double> temp(f.Levels.begin(),f.Levels.end());
c done(temp);
done.out();
return done;
}
size_t nlevels(const factor&f)
{
size_t out=f.Levels.size();
cout<<out<<endl;
return out;
}
//返回输入因子相应在存储中的下标构成的因子
factor as_integer(const factor&f)
{
vector<double> temp(f.i_val.begin(),f.i_val.end());
factor temp0=factor(c(temp));
for_each(temp0.f_val.begin(),temp0.f_val.end(),
[](const pair<double,size_t>&p)
{
cout<<p.first<<" ";
}
);
cout<<endl;
return temp0;
}
factor as_factor(const c&v)
{
factor temp0=factor(v);
temp0.out();
return temp0;
}
//基本运算函数
matrix diag(const c&v)
{
matrix out(v.re_dim(),v.re_dim(),0);
for(size_t i=0;i<v.re_dim();i++)
out[i][i]=v[i];
return out;
}
matrix diag(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
out[i][i]=m[i][i];
return out;
}
double inter_product(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error("dim error!\n");
double out=0;
for(size_t i=0;i<a.re_dim();i++)
out+=(a[i]*b[i]);
return out;
}
c operator%(const matrix&m,const c&v)
{
if(m.re_dim().second!=v.re_dim())
throw runtime_error("dim error!\n");
c out(rep(0,m.re_dim().first));
for(size_t i=0;i<out.re_dim();i++)
out[i]=inter_product(m[i],v);
return out;
}
//同下
double prod(const c&v)
{
double out=1;
for(size_t i=0;i<v.re_dim();i++)
out*=v[i];
return out;
}
//对于matrix的求和由c复制构造函数的重载函数定义：c(const matrix&m)
double sum(const c&v)
{
double out=v[0];
for(size_t i=1;i<v.re_dim();i++)
out+=v[i];
return out;
}
matrix operator%(const c&v,const matrix&m)
{
if(v.re_dim()!=m.re_dim().first)
throw runtime_error("dim error!\n");
matrix out(1,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().second;i++)
out[0][i]=inter_product(v,m.col(i));
return out;
}
matrix outer_product(const c&a,const c&b)
{
matrix out(a.re_dim(),b.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=a[i]*b[j];
return out;
}
size_t length(const c&v)
{
return v.re_dim();
}
c sort(const c&v)
{
vector<double>out;
for(size_t i=0;i<v.re_dim();i++)
out.push_back(v[i]);
sort(out.begin(),out.end());
return c(out);
}
double mean(const c&v)
{
return sum(v)/length(v);
}
//其转换版本返回矩阵转为向量的方差
double var(const c&v)
{
return mean((v-mean(v))*(v-mean(v)));
}
c connect(initializer_list<c> lst)
{
vector<double>temp;
for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});
return c(temp);
}
//标准差
template<typename any>
double STD(const any&a)
{
return sqrt(var(a)*(length(a))/(length(a)-1));
}
matrix rbind(initializer_list<c>lst)
{
size_t r=0;
auto a=lst.begin();
while(a!=lst.end())
{
a++;
r++;
}
size_t co=0;
for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});
matrix out(r,co,0);
size_t i=0;
a=lst.begin();
while(a!=lst.end())
{
vector<double>temp(co,0);
for(size_t i=0;i<a->const_re_c_val().size();i++)
temp[i]=a->operator;
out[i]=c(temp);
i++;
a++;
}
return out;
}
size_t dim(const c&v)
{
return v.re_dim();
}
pair<size_t,size_t> dim(const matrix&m)
{
return m.re_dim();
}
double min(const c&v)
{
double temp=v[0];
for(size_t i=1;i<length(v);i++)
if(v[i]<temp)
temp=v[i];
return temp;
}
double max(const c&v)
{
double temp=v[0];
for(size_t i=1;i<length(v);i++)
if(v[i]>temp)
temp=v[i];
return temp;
}
c pmin(initializer_list<c>lst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i<temp.re_dim();i++)
{
temp[i]=min(*a);
a++;
}
return temp;
}
c pmax(initializer_list<c>lst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i<temp.re_dim();i++)
{
temp[i]=max(*a);
a++;
}
return temp;
}
template<typename any>
pair<double,double> range(const any&v)
{
cout<<"MIN: "<<min(v)<<" MAX: "<<max(v)<<endl;
return make_pair(min(v),max(v));
}
double median(const c&v)
{
return sort(v)[length(v)/2];
}
double var(const c&a,const c&b)
{
if(length(a)!=length(b))
throw runtime_error("dim error!\n");
c t1(a-mean(a)),t2(b-mean(b));
return inter_product(t1,t2)/(length(a)-1);
}
double cov(const c&a,const c&b)
{
return var(a,b);
}
double cor(const c&a,const c&b)
{
return var(a,b)/STD(a)/STD(b);
}
matrix var(const matrix&a,const matrix&b)
{
if(a.re_dim().first!=b.re_dim().first)
throw runtime_error("dim error!\n");
matrix temp(a.re_dim().second,b.re_dim().second,0);
for(size_t i=0;i<temp.re_dim().first;i++)
for(size_t j=0;j<temp.re_dim().second;j++)
temp[i][j]=var(a.col(i),b.col(j));
return temp;
}
matrix cov(const matrix&a,const matrix&b)
{
return var(a,b);
}
matrix cor(const matrix&a,const matrix&b)
{
matrix temp=var(a,b);
for(size_t i=0;i<temp.re_dim().first;i++)
for(size_t j=0;j<temp.re_dim().second;j++)
temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));
return temp;
}
matrix cov(const matrix&a)
{
return var(a,a);
}
matrix cor(const matrix&a)
{
return cor(a,a);
}
//由于关联容器的删除方法利用误差排序的方法
//利用顺序容器迭代器的删除方法有困难
c rank(const c&v)
{
vector<double>temp00(v.const_re_c_val());
for(size_t i=0;i<temp00.size();i++)
temp00[i]+=(0.00000001*i);
vector<double>temp0(temp00);
multiset<double>temp2(temp0.begin(),temp0.end());
vector<double>temp1(v.re_dim());
size_t i=0,j=0;
double pre_val=min(c(vector<double>(temp2.begin(),temp2.end())));
vector<double>::iterator s0;
while(i<v.re_dim())
{
double val=min(c(vector<double>(temp2.begin(),temp2.end())));
s0=find(temp0.begin(),temp0.end(),val);
if(sqrt((pre_val-val)*(pre_val-val))>0.000001)
j++;
temp1[s0-temp0.begin()]=j;
temp2.erase(val);
i++;
pre_val=val;
}
return c(temp1);
}
c rev(const c&v)
{
vector<double> temp;
auto a=v.const_re_c_val().rbegin();
while(a!=v.const_re_c_val().rend())
{
temp.push_back(*a);
a++;
}
return c(temp);
}
//向量先排序再输出排序前的下标值
c order(const c&v)
{
unordered_multimap<double,double>m0;
for(size_t i=0;i<length(v);i++)
m0.insert(make_pair(v[i],i));
multimap<double,double>m1(m0.begin(),m0.end());
c temp(rep(0,length(v)));
auto a=m1.begin();
size_t i=0;
while(a!=m1.end())
{
temp[i]=a->second;
i++;
a++;
}
return temp;
}
c cumsum(const c&v)
{
c temp(v);
for(size_t i=1;i<dim(v);i++)
temp[i]+=temp[i-1];
return temp;
}
c cumprod(const c&v)
{
c temp(v);
for(size_t i=1;i<length(v);i++)
temp[i]*=temp[i-1];
return temp;
}
matrix cumsum(const matrix&m)
{
c temp(m);
for(size_t i=1;i<length(m);i++)
temp[i]+=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
matrix cumprod(const matrix&m)
{
c temp(m);
for(size_t i=1;i<length(m);i++)
temp[i]*=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
//矩阵运算函数
//初等行变换化上三角求解行列式
double det(const matrix&m)
{
auto f=[](const matrix&mm)->pair<matrix,double>
{
if(mm.re_dim().first!=mm.re_dim().second)
throw runtime_error("dim error!\n");
size_t MIN=mm.re_dim().first;
double mult=1;
size_t col=0;
size_t bar=0;
size_t b_c=0;
matrix t(mm);
while(b_c<MIN)
{
for(size_t i=bar+1;i<MIN;)
{
if(t[i][col]==0)
i++;
else
{
t.Swap(i,bar);
mult*=(-1);
break;
}
}
if(bar!=MIN+1)
{
if(t[bar][col]!=0)
{
for(size_t i=bar+1;i<MIN;i++)
{
t[bar][i]=t[bar][i]/t[bar][col];
if(i==bar+1)
mult*=t[bar][col];
}
if(col!=MIN-1)
t[bar][col]=1;
}
for(size_t i=bar+1;i<MIN;i++)
{
t[i]=t[i]-t[i][col]*t[bar];
}
}
col++;
bar++;
b_c=bar;
if(col>bar);
b_c=col;
}
return make_pair(t,mult);
};
pair<matrix,double> tt=f(m);
return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));
}
//Cramer法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0
c solve(const matrix&m,const c&v)
{
c temp(rep(0,v.re_dim()));
double low=det(m);
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
if(low==0)
throw runtime_error("singular!\n");
for(size_t i=0;i<v.re_dim();i++)
{
matrix temp0(m);
temp0.col_change(i,v);
temp[i]=det(temp0)/low;
}
return temp;
}
double tr(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));
}
matrix T(const matrix&m)
{
matrix out(m.re_dim().second,m.re_dim().first,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=m[j][i];
return out;
}
matrix cbind(initializer_list<c>lst)
{
return T(rbind(lst));
}
matrix operator%(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first,b.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=inter_product(a[i],b.col(j));
return out;
}
double F_norm(const matrix&m)
{
return sqrt(tr(T(m)%m));
}
//a,b为区间 d为步长 e为精度 默认为求整数精度：以1为精度特征值 如所求特征值太少可通过减少d改善
//只能求实值
//可以用F_norm作为a、b的界（-F_norm(m),F_norm(m)）可确保值在范围中
//可用于求解但效率低
c eigen(const matrix&m,const double&a,const double &b,const double &d=0.0001,const double &e=1)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
vector<double>temp;
double pre_beg=100000000;
for(double beg=a;beg<b;beg+=d)
{
matrix mm0(diag(rep(beg,m.re_dim().first))-m);
double s0=det(mm0);
if(sqrt(s0*s0)<e)
{
if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)
{
temp.push_back(beg);
pre_beg=beg;
cout<<"Done!\n";
cout<<beg<<endl;//eigenvalue
}
}
}
return c(temp);
}
//重载求eign的函数 d:步长 可保证绝对求出所有特征值（利用循环增大精度）
c eigen(const matrix&m,const double&d=0.1)
{
double dd=d;
double ra=F_norm(m);
double trr=tr(m);
c temp;
double r=0;
do
{
dd*=0.1;//原来设为0.1
temp=eigen(m,-ra,ra,dd);
if(temp.re_dim()==m.re_dim().first)
break;
r=sum(temp)-trr;
}
while(sqrt((r*r)>1));
//与迹的差距（精度）设定为1 绝对值小于此值的特征值被忽略
//此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握
//只适用于所有特征值为实值的情况
//如有复特征值，只输出实数的并循环不能结束
return temp;
}
//矩阵的余子阵
matrix M(const matrix&m,size_t k,size_t l)
{
matrix temp(m);
set<size_t> s1;
for(size_t i=0;i<m.re_dim().first;i++)
s1.insert(i);
set<size_t> s2;
for(size_t i=0;i<m.re_dim().second;i++)
s2.insert(i);
s1.erase(k);
s2.erase(l);
vector<double>temp00;
for(auto j:s2)
for(auto i:s1)
temp00.push_back(temp[i][j]);
return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);
}
//方阵的伴随矩阵
matrix A_accompany(const matrix&m)
{
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);
//代数余子式
return T(out);
}
//求方阵的逆矩阵（利用伴随矩阵方法）
matrix solve(const matrix&m)
{
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error("singular!\n");
return A_accompany(m)/de;
}
//利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆
c solve_eq(const matrix&m,const c&v)
{
c out;
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error("singular!\n");
return solve(m)%v;
}
//化列向量组构成的矩阵为正交化后列向量组构成的矩阵（Schmidt 正交化）
matrix orth_normal(const matrix &m)
{
matrix out(rep(0,m.re_dim().first*m.re_dim().second),m.re_dim().first,m.re_dim().second);
vector<c>dec;
dec.push_back(m.col(0)/sqrt(inter_product(m.col(0),m.col(0))));
for(size_t i=1;i<m.re_dim().second;i++)
{
c temp=m.col(i);
for(size_t j=1;j<=i;j++)
temp=temp-inter_product(dec[j-1],m.col(i))*dec[j-1];
temp=temp/sqrt(inter_product(temp,temp));
dec.push_back(temp);
}
for(size_t i=0;i<m.re_dim().second;i++)
out.col_change(i,dec[i]);
return out;
}
//生成相应位数的由0/1字符构成的字符串
vector<string> generator0_1(size_t n)
{
list<string>out;
list<string>require;
out.push_back("0");
out.push_back("1");
size_t i=0;
while(i<n-1)
{
if(i>0)
out=require;
size_t s0=out.size();
list<string>::iterator pre_begin=out.begin();
auto pre_end=out.end();
size_t k=0;
while(pre_begin!=pre_end)
{
pre_begin++;
k++;
}
pre_begin=out.begin();
for(size_t i=0;i<k;i++)
{
string temp0=*pre_begin+"0";
string temp1=*pre_begin+"1";
out.push_back(temp0);
out.push_back(temp1);
pre_begin++;
}
require.clear();
list<string>::iterator a=out.begin();
for(size_t i=0;i<s0;i++)
a++;
while(a!=out.end())
{
require.push_back(*a);
a++;
}
i++;
}
return vector<string>(require.begin(),require.end());
}
//利用矩阵向量方法生成上述数组
matrix matrix_generator0_1(size_t n)
{
matrix temp0(2,1,0);
matrix temp1;
temp0[0][0]=0;
temp0[1][0]=1;
size_t i=0;
while(i<n-1)
{
if(i>0)
temp0=temp1;
// size_t s0=temp0.re_dim().second;
size_t pre_begin=0;
auto pre_end=temp0.re_dim().first;
temp1=matrix(temp0.re_dim().first*2,temp0.re_dim().second+1,0);
size_t temp1_index=0;
while(pre_begin!=pre_end)
{
vector<double> tt0=temp0[pre_begin].const_re_c_val();
vector<double> tt1=tt0;
tt0.push_back(0);
tt1.push_back(1);
temp1[temp1_index]=c(tt0);
temp1_index++;
temp1[temp1_index]=c(tt1);
temp1_index++;
pre_begin++;
}
i++;
}
return temp1;
}
//将固定位数的下标值与0/1数值构成的map以vector的形式返回
//每个map都是一种组合
vector<map<size_t,size_t>> index_map_generator(size_t n)
{
vector<string>s0(generator0_1(n));
vector<map<size_t,size_t>>outer;
for(size_t i=0;i<s0.size();i++)
{
map<size_t,size_t>tp;
for(size_t j=0;j<s0[i].size();j++)
{
size_t temp=1;
if(s0[i][j]=='0')
temp=0;
tp.insert(make_pair(j,temp));
}
outer.push_back(tp);
}
return outer;
}
//上述产生方式的矩阵版本
vector<map<size_t,size_t>> index_map_matrix_generator(size_t n)
{
matrix temp=matrix_generator0_1(n);
vector<map<size_t,size_t>> out(temp.re_dim().first);
for(size_t i=0;i<out.size();i++)
{
map<size_t,size_t>m;
for(size_t j=0;j<temp.re_dim().second;j++)
m.insert(make_pair(j,temp[i][j]));
out[i]=m;
}
return out;
}
//返回上述index_map_generator(m)中指向0_1和为n的子集
vector<map<size_t,size_t>> require_dim(size_t m,size_t n)
{
vector<map<size_t,size_t>>v0=index_map_generator(m);
vector<map<size_t,size_t>>out;
auto a=v0.begin();
while(a!=v0.end())
{
size_t temp=0;
for_each(a->begin(),a->end(),[&](const pair<size_t,size_t>&p){
if(p.second==1)
temp++;
});
if(temp==n)
out.push_back(*a);
a++;
}
return out;
}
//matrix 所有n阶子式构成的向量
vector<double> seq_of_subdet(const matrix&m,size_t n)
{
vector<double>out;
vector<map<size_t,size_t>>r_map=require_dim(m.re_dim().first,n);
vector<map<size_t,size_t>>c_map=require_dim(m.re_dim().second,n);
for(size_t i=0;i<r_map.size();i++)
for(size_t j=0;j<c_map.size();j++)
{
set<size_t>r_require;
for_each(r_map[i].begin(),r_map[i].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)r_require.insert(p.first);});
set<size_t>c_require;
for_each(c_map[j].begin(),c_map[j].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)c_require.insert(p.first);});
matrix r_matrix(n,m.re_dim().second,0);
auto a=r_require.begin();
for(size_t i=0;i<r_matrix.re_dim().first;i++)
{
r_matrix[i]=m[*a];
a++;
}
matrix c_matrix(r_matrix.re_dim().first,n,0);
auto b=c_require.begin();
for(size_t i=0;i<n;i++)
{
c_matrix.col_change(i,r_matrix.col(*b));
b++;
}
out.push_back(det(c_matrix));
}
return out;
}
//利用子式的方法求出矩阵的秩
size_t rank(const matrix&m)
{
size_t MAX=m.re_dim().first;
if(m.re_dim().second<m.re_dim().first)
MAX=m.re_dim().second;
size_t out=MAX;
size_t pre_num_of_non0=0;
for(size_t i=MAX;i>0;)
{
vector<double>temp=seq_of_subdet(m,i);
auto a=temp.begin();
while(a!=temp.end())
{
a=sqrt((*a)(*a));
a++;
}
size_t num_of_non0=0;
for_each(temp.begin(),temp.end(),[&](const double&d){if(d>0.00000001)num_of_non0++;});
if(pre_num_of_non0==0)
i--;
if(num_of_non0!=0)
return i+1;
pre_num_of_non0=num_of_non0;
}
return out;
}
//回归系数的求解 要求增广后的数据阵（含常系数）可逆
c lm(const c&y,const matrix&m)
{
matrix X=matrix(m.re_dim().first,m.re_dim().second+1,0);
X.col_change(0,rep(1,X.re_dim().first));
for(size_t i=0;i<m.re_dim().second;i++)
X.col_change(i+1,m.col(i));
if(det(T(X)%X)<0.00000001)
throw runtime_error("singular!\n");
c out(solve(T(X)%X)%T(X)%y);
cout<<"The intercept is: "<<out[0]<<endl;
cout<<"Others are:\n";
for(size_t i=1;i<out.re_dim();i++)
cout<<"x"<<i<<" is "<<out[i]<<endl;
cout<<endl;
return out;
}
//抽样函数 参数为向量v 抽取的个数n 是否有放回replace 抽取的概率prob
//使用此方法前必须在函数调用外设定随机数种子
//要求没有0概率点
c sample(const c&v,size_t n,bool replace=1,const c&prob=c())
{
if(!replace&&v.re_dim()<n)
throw out_of_range("n is too large\n");
//抽样概率
vector<double>probb;
if(prob==c())
probb=rep(1/(double)v.re_dim(),v.re_dim()).const_re_c_val();
else
{
probb=prob.const_re_c_val();
double total=0;
for(auto a:probb)
total+=a;
for(auto &a:probb)
a/=total;
}
ifstream fin("h://runif0.txt");
vector<double>runif;
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
//c(probb).out();
c cprobb=cumsum(c(probb));
//cprobb.out();
c require(rep(0,n));
if(replace)
{
size_t i=0;
while(i<n)
{
double uniform=runif[rand()%30000];
//cout<<uniform<<endl;
double conclusion=v[v.re_dim()-1];
for(size_t i=0;cprobb[i]<1;i++)
{
if(uniform<cprobb[i])
{
conclusion=v[i];
break;
}
}
require[i]=conclusion;
i++;
}
}
else
{
size_t i=0;
map<double,double>mdd;
for(size_t i=0;i<v.re_dim();i++)
mdd.insert(make_pair(v[i],c(probb)[i]));
while(i<n)
{
vector<double> temp_v;
vector<double> temp_prob;
for_each(mdd.begin(),mdd.end(),[&](const pair<double,double>&p){
temp_v.push_back(p.first);
temp_prob.push_back(p.second);
});
double d=sample(c(temp_v),1,1,c(temp_prob)).const_re_c_val()[0];
require[i]=d;
mdd.erase(d);
i++;
}
}
return require;
}
//产生几何分布随机数 n 产生个数 p成功概率
c rgeom(const size_t& n,const double& p)
{
vector<double>out;
size_t i=0;
srand(0);
ifstream fin("h://runif0.txt");
vector<double>runif;
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
while(i<n)
{
double temp=runif[rand()%30000];
out.push_back(size_t(log(temp)/log(1-p))+1);
i++;
}
return c(out);
}
//possion分布随机数 n 生成个数 d 分布参数
c rpois(const size_t&n,const double&d)
{
size_t i=0;
vector<double>conclusion;
srand(time(0));
ifstream fin("h://runif0.txt");
vector<double>runif;//储存30000个均匀随机数的vector
double val;
&nbs
de9e
p; while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
while(i<n)
{
double p=exp(-1*d);
size_t out=0;
double F=p;
double temp=runif[rand()%30000];
while(temp>=F)
{
p=d/(double)(out+1)*p;
F+=p;
out++;
}
conclusion.emplace_back(out);
i++;
}
//for_each(conclusion.begin(),conclusion.end(),[](const double&s){cout<<s<<" ";});
//cout<<endl;
return c(conclusion);
}
//利用间距生成向量
c seq(const double&from,const double&to,const double&by)
{
if(from>to)
throw runtime_error("from > to\n");
vector<double>out;
auto begin=from;
while(begin<=to)
{
out.push_back(begin);
begin+=by;
}
return c(out);
}
//生成列表向量元素的所有组合并将结果以矩阵形式输出
//前面的generate_grid_ 函数为铺垫而设 expand_grid函数接受c为参数给出输出
matrix generate_grid_0(pair<double,c>p0)
{
matrix out(p0.second.re_dim(),2,0);
for(size_t i=0;i<out.re_dim().first;i++)
{
out[i][0]=p0.first;
out[i][1]=p0.second[i];
}
return out;
}
matrix generate_grid_1(const c&c0,const c&c1)
{
vector<double> temp0=c0.const_re_c_val();
size_t len=temp0.size();
matrix out(len*c1.re_dim(),2,0);
for(size_t i=0;i<len;i++)
{
matrix temp00=generate_grid_0(make_pair(temp0[i],c1));
for(size_t j=0;j<c1.re_dim();j++)
{
out[i*c1.re_dim()+j]=temp00[j];
}
}
return out;
}
matrix generate_grid_2(const matrix&m,const c&v)
{
matrix temp=generate_grid_1(m.col(m.re_dim().second-1),v);
//temp.out();
//system("pause");
matrix out(temp.re_dim().first,m.re_dim().second+1,0);
//out.out();
//system("pause");
matrix add(m.re_dim().first,m.re_dim().second-1,0);
for(size_t i=0;i<add.re_dim().second;i++)
add.col_change(i,m.col(i));
//add.out();
//system("pause");
for(size_t i=0;i<m.re_dim().first;i++)
for(size_t j=0;j<v.re_dim();j++)
{
for(size_t k=0;k<add.re_dim().second;k++)
out[i*v.re_dim()+j][k]=add[i][k];
}
out.col_change(out.re_dim().second-2,temp.col(0));
out.col_change(out.re_dim().second-1,temp.col(1));
return out;
}
//若干向量其元素的所有组合
matrix expand_grid(initializer_list<c>i0)
{
vector<c>temp(i0);
if(temp.size()==1)
return matrix(temp[0]);
matrix beg0=generate_grid_1(temp[0],temp[1]);
if(temp.size()==2)
return beg0;
for(size_t i=2;i<temp.size();i++)
{
//beg0.out();
beg0=generate_grid_2(beg0,temp[i]);
//cout<<"Done!\n";
//beg0.out();
//cout<<"Done!\n";
}
return beg0;
}
//得到向量v的所有n个元素的组合 输出以矩阵排列 每列为一个组合
//调用了求得下标组合的 require_dim()
//当组合项数太多时主张利用转置按行查看
matrix combn(const c&v,size_t n)
{
auto a=require_dim(v.re_dim(),n);
auto b=vector<set<size_t>>(a.size());
auto a0=a.begin();
auto b0=b.begin();
while(a0!=a.end())
{
for_each(a0->begin(),a0->end(),
[&](const pair<size_t,size_t>&p)
{
if(p.second==1)
b0->insert(p.first);
});
a0++;
b0++;
}
auto cc=vector<vector<double>>(b.size());
for(size_t i=0;i<cc.size();i++)
{
cc[i]=vector<double>(b[i].begin(),b[i].end());
}
matrix out(cc[0].size(),cc.size(),0);
for(size_t i=0;i<out.re_dim().second;i++)
{
out.col_change(i,cc[i]);
}
matrix oout(out.re_dim().first,out.re_dim().second,0);
for(size_t i=0;i<oout.re_dim().first;i++)
for(size_t j=0;j<oout.re_dim().second;j++)
oout[i][j]=v[out[i][j]];
return oout;
}
//上述函数的重载版本 返回0:m的相应n组合
matrix combn(size_t m,size_t n)
{
vector<double>temp;
for(size_t i=0;i<m+1;i++)
temp.push_back(i);
return combn(temp,n);
}
//向量相应元素在前面是否出现
//出现显示1
c duplicated(const c&v)
{
//函数a返回向量cc的前n个元素构成的向量
auto a=[](const c&cc,size_t n)->c
{
c out(rep(0,n));
for(size_t i=0;i<n;i++)
out[i]=cc[i];
return out;
};
vector<c>vc0(v.re_dim());
for(size_t i=0;i<vc0.size();i++)
vc0[i]=a(v,i+1);
vector<vector<double>>vs0(vc0.size());
for(size_t i=0;i<vs0.size();i++)
vs0[i]=vector<double>(vc0[i].const_re_c_val());
c out(rep(0,v.re_dim()));
for(size_t i=1;i<v.re_dim();i++)
{
if(find(vs0[i-1].begin(),vs0[i-1].end(),v[i])!=vs0[i-1].end())
out[i]=1;
}
return out;
}
//实现unique的操作看似与factor的levels相同
//实际在序上是不同的
//levels保持double的序（由set有序决定）
//unique应保持原向量的序
//unique由duplicated定义
c unique(const c&v)
{
c d0=duplicated(v);
vector<pair<double,double>>vm0;
for(size_t i=0;i<v.re_dim();i++)
vm0.push_back(make_pair(v[i],d0[i]));
vector<double>out;
for(size_t i=0;i<vm0.size();i++)
{
if(vm0[i].second==0)
out.push_back(vm0[i].first);
}
return c(out);
}
//矩阵叉乘 T(a)%b (向量内积的矩阵版本)
matrix crossprod(const matrix&a,const matrix&b)
{
return T(a)%b;
}
//outer 生成外积 相应的可将外机作用于函数
matrix outer(const c&v1,const c&v2,function<double(const double&,const double&)>p=function<double(const double&,const double&)>())
{
if(!p)
{
matrix out(v1.re_dim(),v2.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=v1[i]*v2[j];
return out;
}
else
{
matrix out(v1.re_dim(),v2.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=p(v1[i],v2[j]);
return out;
}
}
//矩阵Kronecker积
matrix kronecker(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first*b.re_dim().first,a.re_dim().second*b.re_dim().second,0);
vector<vector<matrix>>vvm0(a.re_dim().first,vector<matrix>(a.re_dim().second));
for(size_t i=0;i<a.re_dim().first;i++)
for(size_t j=0;j<a.re_dim().second;j++)
{
vvm0[i][j]=a[i][j]*b;
}
for(size_t i=0;i<vvm0.size();i++)
for(size_t j=0;j<vvm0[0].size();j++)
{
matrix temp=vvm0[i][j];
for(size_t k=0;k<temp.re_dim().first;k++)
for(size_t l=0;l<temp.re_dim().second;l++)
out[i*b.re_dim().first+k][j*b.re_dim().second+l]=temp[k][l];
}
return out;
}
//由方阵上下三角部分组成的0-1阵 不包含对角线元素
matrix lower_tri(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=1;i<out.re_dim().first;i++)
for(size_t j=0;j<i;j++)
out[i][j]=1;
return out;
}
matrix upper_tri(const matrix&m)
{
return T(lower_tri(m));
}
//矩阵拉直(按列)
c vec(const matrix&m)
{
c out=m.col(0);
for(size_t i=1;i<m.re_dim().second;i++)
out=connect({out,m.col(i)});
return out;
}
//矩阵半拉直 拉直运算忽略主对角线以上部分
c vech(const matrix&m)
{
c out=m.col(0);
for(size_t i=1;i<m.re_dim().second;i++)
{
vector<double> temp=m.col(i).const_re_c_val();
auto a=temp.begin();
for(size_t j=0;j<i;j++)
a++;
vector<double>temp0(a,temp.end());
out=connect({out,c(temp0)});
}
return out;
}
//在方阵满秩情况下 给出方阵A的QR分解
pair<matrix,matrix> qr_of_full_rank(const matrix&A)
{
if(A.re_dim().first!=A.re_dim().second)
throw runtime_error("dim error!\n");
// /*
size_t rk=rank(A);
cout<<"The rank: "<<rk<<endl;
if(rk!=A.re_dim().first)
throw runtime_error("singular!\n");
//*/
cout<<"A=QR"<<endl;
cout<<endl;
cout<<"Q of A:\n";
matrix Q=orth_normal(A);
Q.out();
cout<<"R of m0:\n";
matrix R=T(Q)%A;
R.val_to_double(lower_tri(R),0);
R.out();
// /*
cout<<"The F_norm gap between A and QR:\n";
cout<<F_norm(A-Q%R)<<endl;
//*/
return make_pair(Q,R);
}
//要求自变量矩阵为可逆方阵 利用qr分解 求解线性方程组
//这只是一个形式上的实现 因为求解R逆的方法是用伴随矩阵法
//如能改为初等行变换法应会提高效率
c qr_solve(const matrix&m,const c&v)
{
auto qr=qr_of_full_rank(m);
return solve(qr.second)%T(qr.first)%v;
}
int main()
{
//测试程序：
matrix m2(c({5,0,10,11,0,4,4,3,10,4,3,3,11,3,3,2,10,10,10,10,11,12,13,14,15}),5,5);
m2.out();
c v1({1,1,1,1,2});
c v2=solve(m2,v1);
c v3=qr_solve(m2,v1);
v2.out();
v3.out();
cout<<"The gap of F_norm:\n";
cout<<F_norm(v2-v3)<<endl;
return 0;
}
转载自斯温的博客 http://blog.csdn.net/sinat_30665603
[cpp] view plain copy print?
class_of_R.h#ifndef CLASS_OF_R_H_
#define CLASS_OF_R_H_
#include<iostream>
using std::cout;
using std::endl;
using std::cin;
using std::istream;
using std::ios_base;
#include<vector>
using std::vector;
#include<initializer_list>
using std::initializer_list;
#include<algorithm>
using std::for_each;
using std::sort;
using std::copy;
using std::find;
#include<stdexcept>
using std::runtime_error;
using std::out_of_range;
#include<cstdlib>
using std::srand;
using std::rand;
using std::system;
#include<ctime>
using std::time;
#include<cmath>
using std::sqrt;
using std::log;
using std::exp;
#include<utility>
using std::pair;
using std::make_pair;
#include<iomanip>
using std::setw;
#include<set>
using std::multiset;
using std::set;
#include<unordered_map>
using std::unordered_map;
using std::unordered_multimap;
#include<map>
using std::multimap;
using std::map;
#include<string>
using std::string;
#include<fstream>
using std::ifstream;
#include<memory>
using std::shared_ptr;
#include<list>
using std::list;
class matrix;
//类c声明
class c
{
private:
vector<double> c_val;
size_t dim;
public:
c():c_val(vector<double>()),dim(0){}
c(initializer_list<double>lst);
c(const vector<double>&c_v):c_val(c_v),dim(c_v.size()){}
c(const matrix&m);
void out(size_t n=1)const{cout<<“The vector is:\n”;for(size_t i=0;i<dim;i++)cout<<setw(n)<<c_val[i]<<“ ”;cout<<endl<<endl;}
double operator[](size_t i)const{return c_val[i];}
double& operator[](size_t i){return c_val[i];}
size_t re_dim()const{return dim;}
friend c operator+(const c&a,const c&b);
friend c operator-(const c&a);
friend c operator-(const c&a,const c&b);
friend c rep(const double&d,size_t n);
friend c operator+(const c&a,const double&d);
friend c operator-(const c&a,const double&d);
friend c operator+(const double&d,const c&a);
friend c operator-(const double&d,const c&a);
friend c operator*(const double&d,const c&a);
friend c operator*(const c&a,const double&d);
friend c operator/(const double&d,const c&a);
friend c operator/(const c&a,const double&d);
friend c operator*(const c&a,const c&b);
friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}
friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}
vector<double> const_re_c_val()const{return c_val;}
//在使用const_re_c_val()时可进行拷贝后使用
};
//类matrix声明
class matrix
{
private:
vector<c> m_val;
size_t rdim;
size_t cdim;
public:
matrix():m_val(vector<c>()),rdim(0),cdim(0){}
matrix(const vector<c>&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}
matrix(size_t m,size_t n,const double&d):m_val(vector<c>(m,vector<double>(n,d))),rdim(m),cdim(n){}
matrix(const c&v,size_t r,size_t c);
matrix(const c&v);
c operator[](size_t i)const{return m_val[i];}
c& operator[](size_t i){return m_val[i];}
void out(size_t n=3)const;
c col(size_t i)const;
void col_change(size_t i,const c&a);
friend matrix operator+(const matrix&a,const matrix&b);
friend matrix operator-(const matrix&a);
friend matrix operator-(const matrix&a,const matrix&b);
friend matrix operator+(const matrix&a,const double&d);
friend matrix operator+(const double&d,const matrix&a);
friend matrix operator-(const matrix&a,const double&d);
friend matrix operator-(const double&b,const matrix&a);
friend matrix operator*(const matrix&a,const double&d);
friend matrix operator*(const double&d,const matrix&a);
friend matrix operator/(const matrix&a,const double&d);
friend matrix operator/(const double&d,const matrix&a);
void Swap(size_t i,size_t j)
{c temp=this->operator[](i);this->operator[](i)=this->operator[](j);this->operator[](j)=temp;}
void col_Swap(size_t i,size_t j)
{
auto temp=this->col(i);
this->col_change(i,this->col(j));
this->col_change(j,temp);
}
friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}
friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}
pair<size_t,size_t> re_dim()const{return make_pair(rdim,cdim);}
friend c operator%(const matrix&m,const c&v);
friend matrix operator%(const matrix&a,const matrix&b);
friend matrix operator%(const c&v,const matrix&m);
};
//类factor声明
class factor
{
private:
vector<pair<double,size_t>>f_val;
map<double,size_t>f_val_to_levels;
vector<size_t>i_val;
set<double>Levels;
public:
factor(const c&v=c());
void out(size_t n=1);
friend c levels(const factor&f);
friend size_t nlevels(const factor&f);
friend factor as_integer(const factor&f);
friend factor as_factor(const c&v);
};
#endif // CLASS_OF_R_H_
类成员函数定义文件：
class_of_R.cpp
#include”class_of_R.h”
//类c的有关定义
c::c(initializer_list<double>lst)
{
vector<double>out;
for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});
c_val=out;
dim=out.size();
}
c::c(const matrix&m)
{
vector<double>temp;
for(size_t j=0;j<m.re_dim().second;j++)
for(size_t i=0;i<m.re_dim().first;i++)
temp.push_back(m[i][j]);
c_val=temp;
dim=temp.size();
}
//类matrix的有关定义
void matrix::out(size_t n)const
{
cout<<”The matrix is:\n”;
for(size_t i=0;i<rdim;i++)
{
for(size_t j=0;j<cdim;j++)
cout<<setw(n)<<(*this)[i][j]<<“ ”;
cout<<endl;
}
cout<<endl;
}
c matrix::col(size_t i)const
{
vector<double>temp;
for(size_t i=0;i<rdim;i++)
temp.push_back(0);
c out(temp);
for(size_t j=0;j<rdim;j++)
out[j]=(*this)[j][i];
return out;
}
void matrix::col_change(size_t i,const c&a)
{
for(size_t j=0;j<a.re_dim();j++)
(*this)[j][i]=a[j];
}
matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c)
{
matrix out(r,c,0);
for(size_t i=0;i<v.re_dim();i++)
out[i%r][i/r]=v[i];
(this->m_val)=out.m_val;
}
matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1)
{
matrix out(v.re_dim(),1,0);
for(size_t i=0;i<out.rdim;i++)
out[i][0]=v[i];
this->m_val=out.m_val;
}
//类factor有关定义
factor::factor(const c&v)
{
vector<double>temp=v.const_re_c_val();
for_each(temp.begin(),temp.end(),
[&](const double&d){Levels.insert(d);}
);
vector<double>s_to_v(Levels.begin(),Levels.end());
vector<pair<double,size_t>>ttemp;
for(size_t i=0;i<s_to_v.size();i++)
ttemp.push_back(make_pair(s_to_v[i],i));
for(size_t i=0;i<temp.size();i++)
{
for(size_t j=0;j<ttemp.size();j++)
{
if(temp[i]==ttemp[j].first)
{
f_val.push_back(ttemp[j]);
}
}
}
for(size_t i=0;i<ttemp.size();i++)
f_val_to_levels.insert(ttemp[i]);
auto a=f_val.begin();
for(size_t i=0;i<f_val.size();i++)
{
i_val.push_back(a->second);
a++;
}
}
void factor::out(size_t n)
{
vector<double>temp(f_val.size(),0);
c out0(temp);
auto a=f_val.begin();
for(size_t i=0;i<out0.re_dim();i++)
{
out0[i]=a->first;
a++;
}
vector<double>temp0(Levels.size(),0);
c out1(temp0);
set<double>::iterator b=Levels.begin();
for(size_t i=0;i<out1.re_dim();i++)
{
out1[i]=*b;
b++;
}
out0.out(n);
cout<<”Levels: ”;
cout<<endl;
out1.out(n);
}
其他函数即测试程序文件：
use_R.cpp
#include”class_of_R.h”
//若干友元函数
c operator+(const c&a,const c&b)
{
size_t MAX=b.re_dim();
if(a.re_dim()>b.re_dim())
MAX=a.re_dim();
vector<double>temp0(MAX,0);
vector<double>temp1(MAX,0);
vector<double>temp2(MAX,0);
copy(a.c_val.begin(),a.c_val.end(),temp0.begin());
copy(b.c_val.begin(),b.c_val.end(),temp1.begin());
for(size_t i=0;i<MAX;i++)
temp2[i]=temp0[i]+temp1[i];
return c(temp2);
}
c operator-(const c&a)
{
c out(a);
for(size_t i=0;i<out.re_dim();i++)
out[i]=-out[i];
return out;
}
c operator-(const c&a,const c&b)
{
return (a+(-b));
}
c rep(const double&d,size_t n)
{
return c(vector<double>(n,d));
}
c operator+(const c&a,const double&d)
{
return a+rep(d,a.re_dim());
}
c operator-(const c&a,const double&d)
{
return (a+(-d));
}
c operator+(const double&d,const c&a)
{
return a+d;
}
c operator-(const double&d,const c&a)
{
return (-a)+d;
}
c operator*(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]*=d;
return out;
}
c operator*(const c&a,const double&d)
{
return d*a;
}
c operator/(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]=d/a[i];
return out;
}
c operator/(const c&a,const double&d)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]=a[i]/d;
return out;
}
c operator*(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error(“dim error!\n”);
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]*=b[i];
return out;
}
matrix operator+(const matrix&a,const matrix&b)
{
size_t r_MAX=a.rdim,c_MAX=a.cdim;
if(a.rdim<b.rdim)
r_MAX=b.rdim;
if(a.cdim<b.cdim)
c_MAX=b.cdim;
matrix tempa(r_MAX,c_MAX,0);
for(size_t i=0;i<a.rdim;i++)
for(size_t j=0;j<a.cdim;j++)
tempa[i][j]=a[i][j];
matrix tempb(r_MAX,c_MAX,0);
for(size_t i=0;i<b.rdim;i++)
for(size_t j=0;j<b.cdim;j++)
tempb[i][j]=b[i][j];
matrix tempout(r_MAX,c_MAX,0);
for(size_t i=0;i<tempout.rdim;i++)
for(size_t j=0;j<tempout.cdim;j++)
tempout[i][j]=tempa[i][j]+tempb[i][j];
return tempout;
}
matrix operator-(const matrix&a)
{
matrix out(a);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]=-out[i][j];
return out;
}
matrix operator-(const matrix&a,const matrix&b)
{
return ((-b)+a);
}
matrix operator+(const matrix&a,const double&d)
{
return a+matrix(a.rdim,a.cdim,d);
}
matrix operator+(const double&d,const matrix&a)
{
return a+d;
}
matrix operator-(const matrix&a,const double&d)
{
return a+(-d);
}
matrix operator-(const double&b,const matrix&a)
{
return (-a)+b;
}
matrix operator*(const matrix&a,const double&d)
{
matrix out(a);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]*=d;
return out;
}
matrix operator*(const double&d,const matrix&a)
{
return a*d;
}
matrix operator/(const matrix&a,const double&d)
{
return a*(1/d);
}
matrix operator/(const double&d,const matrix&a)
{
matrix out(a.rdim,a.cdim,0);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]=d/a[i][j];
return out;
}
c levels(const factor&f)
{
vector<double> temp(f.Levels.begin(),f.Levels.end());
c done(temp);
done.out();
return done;
}
size_t nlevels(const factor&f)
{
size_t out=f.Levels.size();
cout<<out<<endl;
return out;
}
//返回输入因子相应在存储中的下标构成的因子
factor as_integer(const factor&f)
{
vector<double> temp(f.i_val.begin(),f.i_val.end());
factor temp0=factor(c(temp));
for_each(temp0.f_val.begin(),temp0.f_val.end(),
[](const pair<double,size_t>&p)
{
cout<<p.first<<” ”;
}
);
cout<<endl;
return temp0;
}
factor as_factor(const c&v)
{
factor temp0=factor(v);
temp0.out();
return temp0;
}
//基本运算函数
matrix diag(const c&v)
{
matrix out(v.re_dim(),v.re_dim(),0);
for(size_t i=0;i<v.re_dim();i++)
out[i][i]=v[i];
return out;
}
matrix diag(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
out[i][i]=m[i][i];
return out;
}
double inter_product(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error(“dim error!\n”);
double out=0;
for(size_t i=0;i<a.re_dim();i++)
out+=(a[i]*b[i]);
return out;
}
c operator%(const matrix&m,const c&v)
{
if(m.re_dim().second!=v.re_dim())
throw runtime_error(“dim error!\n”);
c out(rep(0,m.re_dim().first));
for(size_t i=0;i<out.re_dim();i++)
out[i]=inter_product(m[i],v);
return out;
}
//同下
double prod(const c&v)
{
double out=1;
for(size_t i=0;i<v.re_dim();i++)
out*=v[i];
return out;
}
//对于matrix的求和由c复制构造函数的重载函数定义：c(const matrix&m)
double sum(const c&v)
{
double out=v[0];
for(size_t i=1;i<v.re_dim();i++)
out+=v[i];
return out;
}
matrix operator%(const c&v,const matrix&m)
{
if(v.re_dim()!=m.re_dim().first)
throw runtime_error(“dim error!\n”);
matrix out(1,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().second;i++)
out[0][i]=inter_product(v,m.col(i));
return out;
}
matrix outer_product(const c&a,const c&b)
{
matrix out(a.re_dim(),b.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=a[i]*b[j];
return out;
}
size_t length(const c&v)
{
return v.re_dim();
}
c sort(const c&v)
{
vector<double>out;
for(size_t i=0;i<v.re_dim();i++)
out.push_back(v[i]);
sort(out.begin(),out.end());
return c(out);
}
double mean(const c&v)
{
return sum(v)/length(v);
}
//其转换版本返回矩阵转为向量的方差
double var(const c&v)
{
return mean((v-mean(v))*(v-mean(v)));
}
c connect(initializer_list<c> lst)
{
vector<double>temp;
for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});
return c(temp);
}
//标准差
template<typename any>
double STD(const any&a)
{
return sqrt(var(a)*(length(a))/(length(a)-1));
}
matrix rbind(initializer_list<c>lst)
{
size_t r=0;
auto a=lst.begin();
while(a!=lst.end())
{
a++;
r++;
}
size_t co=0;
for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});
matrix out(r,co,0);
size_t i=0;
a=lst.begin();
while(a!=lst.end())
{
vector<double>temp(co,0);
for(size_t i=0;i<a->const_re_c_val().size();i++)
temp[i]=a->operator[](i);
out[i]=c(temp);
i++;
a++;
}
return out;
}
size_t dim(const c&v)
{
return v.re_dim();
}
pair<size_t,size_t> dim(const matrix&m)
{
return m.re_dim();
}
double min(const c&v)
{
double temp=v[0];
for(size_t i=1;i<length(v);i++)
if(v[i]<temp)
temp=v[i];
return temp;
}
double max(const c&v)
{
double temp=v[0];
for(size_t i=1;i<length(v);i++)
if(v[i]>temp)
temp=v[i];
return temp;
}
c pmin(initializer_list<c>lst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i<temp.re_dim();i++)
{
temp[i]=min(*a);
a++;
}
return temp;
}
c pmax(initializer_list<c>lst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i<temp.re_dim();i++)
{
temp[i]=max(*a);
a++;
}
return temp;
}
template<typename any>
pair<double,double> range(const any&v)
{
cout<<”MIN: ”<<min(v)<<“ MAX: ”<<max(v)<<endl;
return make_pair(min(v),max(v));
}
double median(const c&v)
{
return sort(v)[length(v)/2];
}
double var(const c&a,const c&b)
{
if(length(a)!=length(b))
throw runtime_error(“dim error!\n”);
c t1(a-mean(a)),t2(b-mean(b));
return inter_product(t1,t2)/(length(a)-1);
}
double cov(const c&a,const c&b)
{
return var(a,b);
}
double cor(const c&a,const c&b)
{
return var(a,b)/STD(a)/STD(b);
}
matrix var(const matrix&a,const matrix&b)
{
if(a.re_dim().first!=b.re_dim().first)
throw runtime_error(“dim error!\n”);
matrix temp(a.re_dim().second,b.re_dim().second,0);
for(size_t i=0;i<temp.re_dim().first;i++)
for(size_t j=0;j<temp.re_dim().second;j++)
temp[i][j]=var(a.col(i),b.col(j));
return temp;
}
matrix cov(const matrix&a,const matrix&b)
{
return var(a,b);
}
matrix cor(const matrix&a,const matrix&b)
{
matrix temp=var(a,b);
for(size_t i=0;i<temp.re_dim().first;i++)
for(size_t j=0;j<temp.re_dim().second;j++)
temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));
return temp;
}
matrix cov(const matrix&a)
{
return var(a,a);
}
matrix cor(const matrix&a)
{
return cor(a,a);
}
//由于关联容器的删除方法利用误差排序的方法
//利用顺序容器迭代器的删除方法有困难
c rank(const c&v)
{
vector<double>temp00(v.const_re_c_val());
for(size_t i=0;i<temp00.size();i++)
temp00[i]+=(0.00000001*i);
vector<double>temp0(temp00);
multiset<double>temp2(temp0.begin(),temp0.end());
vector<double>temp1(v.re_dim());
size_t i=0,j=0;
double pre_val=min(c(vector<double>(temp2.begin(),temp2.end())));
vector<double>::iterator s0;
while(i<v.re_dim())
{
double val=min(c(vector<double>(temp2.begin(),temp2.end())));
s0=find(temp0.begin(),temp0.end(),val);
if(sqrt((pre_val-val)*(pre_val-val))>0.000001)
j++;
temp1[s0-temp0.begin()]=j;
temp2.erase(val);
i++;
pre_val=val;
}
return c(temp1);
}
c rev(const c&v)
{
vector<double> temp;
auto a=v.const_re_c_val().rbegin();
while(a!=v.const_re_c_val().rend())
{
temp.push_back(*a);
a++;
}
return c(temp);
}
//向量先排序再输出排序前的下标值
c order(const c&v)
{
unordered_multimap<double,double>m0;
for(size_t i=0;i<length(v);i++)
m0.insert(make_pair(v[i],i));
multimap<double,double>m1(m0.begin(),m0.end());
c temp(rep(0,length(v)));
auto a=m1.begin();
size_t i=0;
while(a!=m1.end())
{
temp[i]=a->second;
i++;
a++;
}
return temp;
}
c cumsum(const c&v)
{
c temp(v);
for(size_t i=1;i<dim(v);i++)
temp[i]+=temp[i-1];
return temp;
}
c cumprod(const c&v)
{
c temp(v);
for(size_t i=1;i<length(v);i++)
temp[i]*=temp[i-1];
return temp;
}
matrix cumsum(const matrix&m)
{
c temp(m);
for(size_t i=1;i<length(m);i++)
temp[i]+=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
matrix cumprod(const matrix&m)
{
c temp(m);
for(size_t i=1;i<length(m);i++)
temp[i]*=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
//矩阵运算函数
//初等行变换化上三角求解行列式
double det(const matrix&m)
{
auto f=[](const matrix&mm)->pair<matrix,double>
{
if(mm.re_dim().first!=mm.re_dim().second)
throw runtime_error(“dim error!\n”);
size_t MIN=mm.re_dim().first;
double mult=1;
size_t col=0;
size_t bar=0;
size_t b_c=0;
matrix t(mm);
while(b_c<MIN)
{
for(size_t i=bar+1;i<MIN;)
{
if(t[i][col]==0)
i++;
else
{
t.Swap(i,bar);
mult*=(-1);
break;
}
}
if(bar!=MIN+1)
{
if(t[bar][col]!=0)
{
for(size_t i=bar+1;i<MIN;i++)
{
t[bar][i]=t[bar][i]/t[bar][col];
if(i==bar+1)
mult*=t[bar][col];
}
if(col!=MIN-1)
t[bar][col]=1;
}
for(size_t i=bar+1;i<MIN;i++)
{
t[i]=t[i]-t[i][col]*t[bar];
}
}
col++;
bar++;
b_c=bar;
if(col>bar);
b_c=col;
}
return make_pair(t,mult);
};
pair<matrix,double> tt=f(m);
return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));
}
//Cramer法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0
c solve(const matrix&m,const c&v)
{
c temp(rep(0,v.re_dim()));
double low=det(m);
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
if(low==0)
throw runtime_error(“singular!\n”);
for(size_t i=0;i<v.re_dim();i++)
{
matrix temp0(m);
temp0.col_change(i,v);
temp[i]=det(temp0)/low;
}
return temp;
}
double tr(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));
}
matrix T(const matrix&m)
{
matrix out(m.re_dim().second,m.re_dim().first,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=m[j][i];
return out;
}
matrix cbind(initializer_list<c>lst)
{
return T(rbind(lst));
}
matrix operator%(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first,b.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=inter_product(a[i],b.col(j));
return out;
}
double F_norm(const matrix&m)
{
return sqrt(tr(T(m)%m));
}
//a,b为区间 d为步长 e为精度 默认为求整数精度：以1为精度特征值 如所求特征值太少可通过减少d改善
//只能求实值
//可以用F_norm作为a、b的界（-F_norm(m),F_norm(m)）可确保值在范围中
//可用于求解但效率低
c eigen(const matrix&m,const double&a,const double &b,const double &d=0.0001,const double &e=1)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
vector<double>temp;
double pre_beg=100000000;
for(double beg=a;beg<b;beg+=d)
{
matrix mm0(diag(rep(beg,m.re_dim().first))-m);
double s0=det(mm0);
if(sqrt(s0*s0)<e)
{
if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)
{
temp.push_back(beg);
pre_beg=beg;
cout<<”Done!\n”;
cout<<beg<<endl;//eigenvalue
}
}
}
return c(temp);
}
//重载求eign的函数 d:步长 可保证绝对求出所有特征值（利用循环增大精度）
c eigen(const matrix&m,const double&d=0.1)
{
double dd=d;
double ra=F_norm(m);
double trr=tr(m);
c temp;
double r=0;
do
{
dd*=0.1;//原来设为0.1
temp=eigen(m,-ra,ra,dd);
if(temp.re_dim()==m.re_dim().first)
break;
r=sum(temp)-trr;
}
while(sqrt((r*r)>1));
//与迹的差距（精度）设定为1 绝对值小于此值的特征值被忽略
//此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握
//只适用于所有特征值为实值的情况
//如有复特征值，只输出实数的并循环不能结束
return temp;
}
//矩阵的余子阵
matrix M(const matrix&m,size_t k,size_t l)
{
matrix temp(m);
set<size_t> s1;
for(size_t i=0;i<m.re_dim().first;i++)
s1.insert(i);
set<size_t> s2;
for(size_t i=0;i<m.re_dim().second;i++)
s2.insert(i);
s1.erase(k);
s2.erase(l);
vector<double>temp00;
for(auto j:s2)
for(auto i:s1)
temp00.push_back(temp[i][j]);
return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);
}
//方阵的伴随矩阵
matrix A_accompany(const matrix&m)
{
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);
//代数余子式
return T(out);
}
//求方阵的逆矩阵（利用伴随矩阵方法）
matrix solve(const matrix&m)
{
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error(“singular!\n”);
return A_accompany(m)/de;
}
//利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆
c solve_eq(const matrix&m,const c&v)
{
c out;
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error(“singular!\n”);
return solve(m)%v;
}
//化列向量组构成的矩阵为正交化后列向量组构成的矩阵（Schmidt 正交化）
matrix orth_normal(const matrix &m)
{
matrix out(rep(0,m.re_dim().first*m.re_dim().second),m.re_dim().first,m.re_dim().second);
vector<c>dec;
dec.push_back(m.col(0)/sqrt(inter_product(m.col(0),m.col(0))));
for(size_t i=1;i<m.re_dim().second;i++)
{
c temp=m.col(i);
for(size_t j=1;j<=i;j++)
temp=temp-inter_product(dec[j-1],m.col(i))*dec[j-1];
temp=temp/sqrt(inter_product(temp,temp));
dec.push_back(temp);
}
for(size_t i=0;i<m.re_dim().second;i++)
out.col_change(i,dec[i]);
return out;
}
//生成相应位数的由0/1字符构成的字符串
vector<string> generator0_1(size_t n)
{
list<string>out;
list<string>require;
out.push_back(”0”);
out.push_back(”1”);
size_t i=0;
while(i<n-1)
{
if(i>0)
out=require;
size_t s0=out.size();
list<string>::iterator pre_begin=out.begin();
auto pre_end=out.end();
size_t k=0;
while(pre_begin!=pre_end)
{
pre_begin++;
k++;
}
pre_begin=out.begin();
for(size_t i=0;i<k;i++)
{
string temp0=*pre_begin+”0”;
string temp1=*pre_begin+”1”;
out.push_back(temp0);
out.push_back(temp1);
pre_begin++;
}
require.clear();
list<string>::iterator a=out.begin();
for(size_t i=0;i<s0;i++)
a++;
while(a!=out.end())
{
require.push_back(*a);
a++;
}
i++;
}
return vector<string>(require.begin(),require.end());
}
//利用矩阵向量方法生成上述数组
matrix matrix_generator0_1(size_t n)
{
matrix temp0(2,1,0);
matrix temp1;
temp0[0][0]=0;
temp0[1][0]=1;
size_t i=0;
while(i<n-1)
{
if(i>0)
temp0=temp1;
// size_t s0=temp0.re_dim().second;
size_t pre_begin=0;
auto pre_end=temp0.re_dim().first;
temp1=matrix(temp0.re_dim().first*2,temp0.re_dim().second+1,0);
size_t temp1_index=0;
while(pre_begin!=pre_end)
{
vector<double> tt0=temp0[pre_begin].const_re_c_val();
vector<double> tt1=tt0;
tt0.push_back(0);
tt1.push_back(1);
temp1[temp1_index]=c(tt0);
temp1_index++;
temp1[temp1_index]=c(tt1);
temp1_index++;
pre_begin++;
}
i++;
}
return temp1;
}
//将固定位数的下标值与0/1数值构成的map以vector的形式返回
//每个map都是一种组合
vector<map<size_t,size_t>> index_map_generator(size_t n)
{
vector<string>s0(generator0_1(n));
vector<map<size_t,size_t>>outer;
for(size_t i=0;i<s0.size();i++)
{
map<size_t,size_t>tp;
for(size_t j=0;j<s0[i].size();j++)
{
size_t temp=1;
if(s0[i][j]==‘0’)
temp=0;
tp.insert(make_pair(j,temp));
}
outer.push_back(tp);
}
return outer;
}
//上述产生方式的矩阵版本
vector<map<size_t,size_t>> index_map_matrix_generator(size_t n)
{
matrix temp=matrix_generator0_1(n);
vector<map<size_t,size_t>> out(temp.re_dim().first);
for(size_t i=0;i<out.size();i++)
{
map<size_t,size_t>m;
for(size_t j=0;j<temp.re_dim().second;j++)
m.insert(make_pair(j,temp[i][j]));
out[i]=m;
}
return out;
}
//返回上述index_map_generator(m)中指向0_1和为n的子集
vector<map<size_t,size_t>> require_dim(size_t m,size_t n)
{
vector<map<size_t,size_t>>v0=index_map_generator(m);
vector<map<size_t,size_t>>out;
auto a=v0.begin();
while(a!=v0.end())
{
size_t temp=0;
for_each(a->begin(),a->end(),[&](const pair<size_t,size_t>&p){
if(p.second==1)
temp++;
});
if(temp==n)
out.push_back(*a);
a++;
}
return out;
}
//matrix 所有n阶子式构成的向量
vector<double> seq_of_subdet(const matrix&m,size_t n)
{
vector<double>out;
vector<map<size_t,size_t>>r_map=require_dim(m.re_dim().first,n);&
4000
nbsp;
vector<map<size_t,size_t>>c_map=require_dim(m.re_dim().second,n);
for(size_t i=0;i<r_map.size();i++)
for(size_t j=0;j<c_map.size();j++)
{
set<size_t>r_require;
for_each(r_map[i].begin(),r_map[i].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)r_require.insert(p.first);});
set<size_t>c_require;
for_each(c_map[j].begin(),c_map[j].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)c_require.insert(p.first);});
matrix r_matrix(n,m.re_dim().second,0);
auto a=r_require.begin();
for(size_t i=0;i<r_matrix.re_dim().first;i++)
{
r_matrix[i]=m[*a];
a++;
}
matrix c_matrix(r_matrix.re_dim().first,n,0);
auto b=c_require.begin();
for(size_t i=0;i<n;i++)
{
c_matrix.col_change(i,r_matrix.col(*b));
b++;
}
out.push_back(det(c_matrix));
}
return out;
}
//利用子式的方法求出矩阵的秩
size_t rank(const matrix&m)
{
size_t MAX=m.re_dim().first;
if(m.re_dim().second<m.re_dim().first)
MAX=m.re_dim().second;
size_t out=MAX;
size_t pre_num_of_non0=0;
for(size_t i=MAX;i>0;)
{
vector<double>temp=seq_of_subdet(m,i);
auto a=temp.begin();
while(a!=temp.end())
{
*a=sqrt((*a)*(*a));
a++;
}
size_t num_of_non0=0;
for_each(temp.begin(),temp.end(),[&](const double&d){if(d>0.00000001)num_of_non0++;});
if(pre_num_of_non0==0)
i–;
if(num_of_non0!=0)
return i+1;
pre_num_of_non0=num_of_non0;
}
return out;
}
//回归系数的求解 要求增广后的数据阵（含常系数）可逆
c lm(const c&y,const matrix&m)
{
matrix X=matrix(m.re_dim().first,m.re_dim().second+1,0);
X.col_change(0,rep(1,X.re_dim().first));
for(size_t i=0;i<m.re_dim().second;i++)
X.col_change(i+1,m.col(i));
if(det(T(X)%X)<0.00000001)
throw runtime_error(“singular!\n”);
c out(solve(T(X)%X)%T(X)%y);
cout<<”The intercept is: ”<<out[0]<<endl;
cout<<”Others are:\n”;
for(size_t i=1;i<out.re_dim();i++)
cout<<”x”<<i<<“ is ”<<out[i]<<endl;
cout<<endl;
return out;
}
//抽样函数 参数为向量v 抽取的个数n 是否有放回replace 抽取的概率prob
//使用此方法前必须在函数调用外设定随机数种子
//要求没有0概率点
c sample(const c&v,size_t n,bool replace=1,const c&prob=c())
{
if(!replace&&v.re_dim()<n)
throw out_of_range(“n is too large\n”);
//抽样概率
vector<double>probb;
if(prob==c())
probb=rep(1/(double)v.re_dim(),v.re_dim()).const_re_c_val();
else
{
probb=prob.const_re_c_val();
double total=0;
for(auto a:probb)
total+=a;
for(auto &a:probb)
a/=total;
}
ifstream fin(”h://runif0.txt”);
vector<double>runif;
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
//c(probb).out();
c cprobb=cumsum(c(probb));
//cprobb.out();
c require(rep(0,n));
if(replace)
{
size_t i=0;
while(i<n)
{
double uniform=runif[rand()%30000];
//cout<<uniform<<endl;
double conclusion=v[v.re_dim()-1];
for(size_t i=0;cprobb[i]<1;i++)
{
if(uniform<cprobb[i])
{
&n
1c143
bsp; conclusion=v[i];
break;
}
}
require[i]=conclusion;
i++;
}
}
else
{
size_t i=0;
map<double,double>mdd;
for(size_t i=0;i<v.re_dim();i++)
mdd.insert(make_pair(v[i],c(probb)[i]));
while(i<n)
{
vector<double> temp_v;
vector<double> temp_prob;
for_each(mdd.begin(),mdd.end(),[&](const pair<double,double>&p){
temp_v.push_back(p.first);
temp_prob.push_back(p.second);
});
double d=sample(c(temp_v),1,1,c(temp_prob)).const_re_c_val()[0];
require[i]=d;
mdd.erase(d);
i++;
}
}
return require;
}
//产生几何分布随机数 n 产生个数 p成功概率
c rgeom(const size_t& n,const double& p)
{
vector<double>out;
size_t i=0;
srand(0);
ifstream fin(”h://runif0.txt”);
vector<double>runif;
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
while(i<n)
{
double temp=runif[rand()%30000];
out.push_back(size_t(log(temp)/log(1-p))+1);
i++;
}
return c(out);
}
//possion分布随机数 n 生成个数 d 分布参数
c rpois(const size_t&n,const double&d)
{
size_t i=0;
vector<double>conclusion;
srand(time(0));
ifstream fin(”h://runif0.txt”);
vector<double>runif;//储存30000个均匀随机数的vector
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
while(i<n)
{
double p=exp(-1*d);
size_t out=0;
double F=p;
double temp=runif[rand()%30000];
while(temp>=F)
{
p=d/(double)(out+1)*p;
F+=p;
out++;
}
conclusion.emplace_back(out);
i++;
}
//for_each(conclusion.begin(),conclusion.end(),[](const double&s){cout<<s<<” ”;});
//cout<<endl;
return c(conclusion);
}
//利用间距生成向量
c seq(const double&from,const double&to,const double&by)
{
if(from>to)
throw runtime_error(“from > to\n”);
vector<double>out;
auto begin=from;
while(begin<=to)
{
out.push_back(begin);
begin+=by;
}
return c(out);
}
//生成列表向量元素的所有组合并将结果以矩阵形式输出
//前面的generate_grid_ 函数为铺垫而设 expand_grid函数接受c为参数给出输出
matrix generate_grid_0(pair<double,c>p0)
{
matrix out(p0.second.re_dim(),2,0);
for(size_t i=0;i<out.re_dim().first;i++)
{
out[i][0]=p0.first;
out[i][1]=p0.second[i];
}
return out;
}
matrix generate_grid_1(const c&c0,const c&c1)
{
vector<double> temp0=c0.const_re_c_val();
size_t len=temp0.size();
matrix out(len*c1.re_dim(),2,0);
for(size_t i=0;i<len;i++)
{
matrix temp00=generate_grid_0(make_pair(temp0[i],c1));
for(size_t j=0;j<c1.re_dim();j++)
{
out[i*c1.re_dim()+j]=temp00[j];
}
}
return out;
}
matrix generate_grid_2(const matrix&m,const c&v)
{
matrix temp=generate_grid_1(m.col(m.re_dim().second-1),v);
//temp.out();
//system(“pause”);
matrix out(temp.re_dim().first,m.re_dim().second+1,0);
//out.out();
//system(“pause”);
matrix add(m.re_dim().first,m.re_dim().second-1,0);
for(size_t i=0;i<add.re_dim().second;i++)
add.col_change(i,m.col(i));
//add.out();
//system(“pause”);
for(size_t i=0;i<m.re_dim().first;i++)
for(size_t j=0;j<v.re_dim();j++)
{
for(size_t k=0;k<add.re_dim().second;k++)
out[i*v.re_dim()+j][k]=add[i][k];
}
out.col_change(out.re_dim().second-2,temp.col(0));
out.col_change(out.re_dim().second-1,temp.col(1));
return out;
}
//若干向量其元素的所有组合
matrix expand_grid(initializer_list<c>i0)
{
vector<c>temp(i0);
if(temp.size()==1)
return matrix(temp[0]);
matrix beg0=generate_grid_1(temp[0],temp[1]);
if(temp.size()==2)
return beg0;
for(size_t i=2;i<temp.size();i++)
{
//beg0.out();
beg0=generate_grid_2(beg0,temp[i]);
//cout<<”Done!\n”;
//beg0.out();
//cout<<”Done!\n”;
}
return beg0;
}
//得到向量v的所有n个元素的组合 输出以矩阵排列 每列为一个组合
//调用了求得下标组合的 require_dim()
//当组合项数太多时主张利用转置按行查看
matrix combn(const c&v,size_t n)
{
auto a=require_dim(v.re_dim(),n);
auto b=vector<set<size_t>>(a.size());
auto a0=a.begin();
auto b0=b.begin();
while(a0!=a.end())
{
for_each(a0->begin(),a0->end(),
[&](const pair<size_t,size_t>&p)
{
if(p.second==1)
b0->insert(p.first);
});
a0++;
b0++;
}
auto cc=vector<vector<double>>(b.size());
for(size_t i=0;i<cc.size();i++)
{
cc[i]=vector<double>(b[i].begin(),b[i].end());
}
matrix out(cc[0].size(),cc.size(),0);
for(size_t i=0;i<out.re_dim().second;i++)
{
out.col_change(i,cc[i]);
}
matrix oout(out.re_dim().first,out.re_dim().second,0);
for(size_t i=0;i<oout.re_dim().first;i++)
for(size_t j=0;j<oout.re_dim().second;j++)
oout[i][j]=v[out[i][j]];
return oout;
}
//上述函数的重载版本 返回0:m的相应n组合
matrix combn(size_t m,size_t n)
{
vector<double>temp;
for(size_t i=0;i<m+1;i++)
temp.push_back(i);
return combn(temp,n);
}
//向量相应元素在前面是否出现
//出现显示1
c duplicated(const c&v)
{
//函数a返回向量cc的前n个元素构成的向量
auto a=[](const c&cc,size_t n)->c
{
c out(rep(0,n));
for(size_t i=0;i<n;i++)
out[i]=cc[i];
return out;
};
vector<c>vc0(v.re_dim());
for(size_t i=0;i<vc0.size();i++)
vc0[i]=a(v,i+1);
vector<vector<double>>vs0(vc0.size());
for(size_t i=0;i<vs0.size();i++)
vs0[i]=vector<double>(vc0[i].const_re_c_val());
c out(rep(0,v.re_dim()));
for(size_t i=1;i<v.re_dim();i++)
{
if(find(vs0[i-1].begin(),vs0[i-1].end(),v[i])!=vs0[i-1].end())
out[i]=1;
}
return out;
}
//实现unique的操作看似与factor的levels相同
//实际在序上是不同的
//levels保持double的序（由set有序决定）
//unique应保持原向量的序
//unique由duplicated定义
c unique(const c&v)
{
c d0=duplicated(v);
vector<pair<double,double>>vm0;
for(size_t i=0;i<v.re_dim();i++)
vm0.push_back(make_pair(v[i],d0[i]));
vector<double>out;
for(size_t i=0;i<vm0.size();i++)
{
if(vm0[i].second==0)
out.push_back(vm0[i].first);
}
return c(out);
}
//矩阵叉乘 T(a)%b (向量内积的矩阵版本)
matrix crossprod(const matrix&a,const matrix&b)
{
return T(a)%b;
}
//outer 生成外积 相应的可将外机作用于函数
matrix outer(const c&v1,const c&v2,function<double(const double&,const double&)>p=function<double(const double&,const double&)>())
{
if(!p)
{
matrix out(v1.re_dim(),v2.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=v1[i]*v2[j];
return out;
}
else
{
matrix out(v1.re_dim(),v2.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=p(v1[i],v2[j]);
return out;
}
}
//矩阵Kronecker积
matrix kronecker(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first*b.re_dim().first,a.re_dim().second*b.re_dim().second,0);
vector<vector<matrix>>vvm0(a.re_dim().first,vector<matrix>(a.re_dim().second));
for(size_t i=0;i<a.re_dim().first;i++)
for(size_t j=0;j<a.re_dim().second;j++)
{
vvm0[i][j]=a[i][j]*b;
}
for(size_t i=0;i<vvm0.size();i++)
for(size_t j=0;j<vvm0[0].size();j++)
{
matrix temp=vvm0[i][j];
for(size_t k=0;k<temp.re_dim().first;k++)
for(size_t l=0;l<temp.re_dim().second;l++)
out[i*b.re_dim().first+k][j*b.re_dim().second+l]=temp[k][l];
}
return out;
}
//由方阵上下三角部分组成的0-1阵 不包含对角线元素
matrix lower_tri(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error(“dim error!\n”);
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=1;i<out.re_dim().first;i++)
for(size_t j=0;j<i;j++)
out[i][j]=1;
return out;
}
matrix upper_tri(const matrix&m)
{
return T(lower_tri(m));
}
//矩阵拉直(按列)
c vec(const matrix&m)
{
c out=m.col(0);
for(size_t i=1;i<m.re_dim().second;i++)
out=connect({out,m.col(i)});
return out;
}
//矩阵半拉直 拉直运算忽略主对角线以上部分
c vech(const matrix&m)
{
c out=m.col(0);
for(size_t i=1;i<m.re_dim().second;i++)
{
vector<double> temp=m.col(i).const_re_c_val();
auto a=temp.begin();
for(size_t j=0;j<i;j++)
a++;
vector<double>temp0(a,temp.end());
out=connect({out,c(temp0)});
}
return out;
}
//在方阵满秩情况下 给出方阵A的QR分解
pair<matrix,matrix> qr_of_full_rank(const matrix&A)
{
if(A.re_dim().first!=A.re_dim().second)
throw runtime_error(“dim error!\n”);
// /*
size_t rk=rank(A);
cout<<”The rank: ”<<rk<<endl;
if(rk!=A.re_dim().first)
throw runtime_error(“singular!\n”);
//*/
cout<<”A=QR”<<endl;
cout<<endl;
cout<<”Q of A:\n”;
matrix Q=orth_normal(A);
Q.out();
cout<<”R of m0:\n”;
matrix R=T(Q)%A;
R.val_to_double(lower_tri(R),0);
R.out();
// /*
cout<<”The F_norm gap between A and QR:\n”;
cout<<F_norm(A-Q%R)<<endl;
//*/
return make_pair(Q,R);
}
//要求自变量矩阵为可逆方阵 利用qr分解 求解线性方程组
//这只是一个形式上的实现 因为求解R逆的方法是用伴随矩阵法
//如能改为初等行变换法应会提高效率
c qr_solve(const matrix&m,const c&v)
{
auto qr=qr_of_full_rank(m);
return solve(qr.second)%T(qr.first)%v;
}
int main()
{
//测试程序：
matrix m2(c({5,0,10,11,0,4,4,3,10,4,3,3,11,3,3,2,10,10,10,10,11,12,13,14,15}),5,5);
m2.out();
c v1({1,1,1,1,2});
c v2=solve(m2,v1);
c v3=qr_solve(m2,v1);
v2.out();
v3.out();
cout<<”The gap of F_norm:\n”;
cout<<F_norm(v2-v3)<<endl;
return 0;
}
class_of_R.h
ifndef CLASS_OF_R_H_
define CLASS_OF_R_H_
include<iostream>
using std::cout;using std::endl;
using std::cin;
using std::istream;
using std::ios_base;
include<vector>
using std::vector;include<initializer_list>
using std::initializer_list;include<algorithm>
using std::for_each;using std::sort;
using std::copy;
using std::find;
include<stdexcept>
using std::runtime_error;using std::out_of_range;
include<cstdlib>
using std::srand;using std::rand;
using std::system;
include<ctime>
using std::time;include<cmath>
using std::sqrt;using std::log;
using std::exp;
include<utility>
using std::pair;using std::make_pair;
include<iomanip>
using std::setw;include<set>
using std::multiset;using std::set;
include<unordered_map>
using std::unordered_map;using std::unordered_multimap;
include<map>
using std::multimap;using std::map;
include<string>
using std::string;include<fstream>
using std::ifstream;include<memory>
using std::shared_ptr;include<list>
using std::list;class matrix;
//类c声明
class c
{
private:
vector<double> c_val;
size_t dim;
public:
c():c_val(vector<double>()),dim(0){}
c(initializer_list<double>lst);
c(const vector<double>&c_v):c_val(c_v),dim(c_v.size()){}
c(const matrix&m);
void out(size_t n=1)const{cout<<"The vector is:\n";for(size_t i=0;i<dim;i++)cout<<setw(n)<<c_val[i]<<" ";cout<<endl<<endl;}
double operator[](size_t i)const{return c_val[i];}
double& operator[](size_t i){return c_val[i];}
size_t re_dim()const{return dim;}
friend c operator+(const c&a,const c&b);
friend c operator-(const c&a);
friend c operator-(const c&a,const c&b);
friend c rep(const double&d,size_t n);
friend c operator+(const c&a,const double&d);
friend c operator-(const c&a,const double&d);
friend c operator+(const double&d,const c&a);
friend c operator-(const double&d,const c&a);
friend c operator*(const double&d,const c&a);
friend c operator*(const c&a,const double&d);
friend c operator/(const double&d,const c&a);
friend c operator/(const c&a,const double&d);
friend c operator*(const c&a,const c&b);
friend bool operator==(const c&a,const c&b){return a.c_val==b.c_val;}
friend bool operator!=(const c&a,const c&b){return a.c_val!=b.c_val;}
vector<double> const_re_c_val()const{return c_val;}
//在使用const_re_c_val()时可进行拷贝后使用
};
//类matrix声明
class matrix
{
private:
vector<c> m_val;
size_t rdim;
size_t cdim;
public:
matrix():m_val(vector<c>()),rdim(0),cdim(0){}
matrix(const vector<c>&m):m_val(m),rdim(m_val.size()),cdim(m_val[0].re_dim()){}
matrix(size_t m,size_t n,const double&d):m_val(vector<c>(m,vector<double>(n,d))),rdim(m),cdim(n){}
matrix(const c&v,size_t r,size_t c);
matrix(const c&v);
c operator[](size_t i)const{return m_val[i];}
c& operator[](size_t i){return m_val[i];}
void out(size_t n=3)const;
c col(size_t i)const;
void col_change(size_t i,const c&a);
friend matrix operator+(const matrix&a,const matrix&b);
friend matrix operator-(const matrix&a);
friend matrix operator-(const matrix&a,const matrix&b);
friend matrix operator+(const matrix&a,const double&d);
friend matrix operator+(const double&d,const matrix&a);
friend matrix operator-(const matrix&a,const double&d);
friend matrix operator-(const double&b,const matrix&a);
friend matrix operator*(const matrix&a,const double&d);
friend matrix operator*(const double&d,const matrix&a);
friend matrix operator/(const matrix&a,const double&d);
friend matrix operator/(const double&d,const matrix&a);
void Swap(size_t i,size_t j)
{c temp=this->operator;this->operator=this->operator;this->operator=temp;}
void col_Swap(size_t i,size_t j)
{
auto temp=this->col(i);
this->col_change(i,this->col(j));
this->col_change(j,temp);
}
friend bool operator==(const matrix&a,const matrix&b){return a.m_val==b.m_val;}
friend bool operator!=(const matrix&a,const matrix&b){return a.m_val!=b.m_val;}
pair<size_t,size_t> re_dim()const{return make_pair(rdim,cdim);}
friend c operator%(const matrix&m,const c&v);
friend matrix operator%(const matrix&a,const matrix&b);
friend matrix operator%(const c&v,const matrix&m);
};
//类factor声明
class factor
{
private:
vector<pair<double,size_t>>f_val;
map<double,size_t>f_val_to_levels;
vector<size_t>i_val;
set<double>Levels;
public:
factor(const c&v=c());
void out(size_t n=1);
friend c levels(const factor&f);
friend size_t nlevels(const factor&f);
friend factor as_integer(const factor&f);
friend factor as_factor(const c&v);
};
endif // CLASS_OF_R_H_
类成员函数定义文件：class_of_R.cpp
include"class_of_R.h"
//类c的有关定义c::c(initializer_list<double>lst)
{
vector<double>out;
for_each(lst.begin(),lst.end(),[&](const double&d){out.push_back(d);});
c_val=out;
dim=out.size();
}
c::c(const matrix&m)
{
vector<double>temp;
for(size_t j=0;j<m.re_dim().second;j++)
for(size_t i=0;i<m.re_dim().first;i++)
temp.push_back(m[i][j]);
c_val=temp;
dim=temp.size();
}
//类matrix的有关定义
void matrix::out(size_t n)const
{
cout<<"The matrix is:\n";
for(size_t i=0;i<rdim;i++)
{
for(size_t j=0;j<cdim;j++)
cout<<setw(n)<<(*this)[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
c matrix::col(size_t i)const
{
vector<double>temp;
for(size_t i=0;i<rdim;i++)
temp.push_back(0);
c out(temp);
for(size_t j=0;j<rdim;j++)
out[j]=(*this)[j][i];
return out;
}
void matrix::col_change(size_t i,const c&a)
{
for(size_t j=0;j<a.re_dim();j++)
(*this)[j][i]=a[j];
}
matrix::matrix(const c&v,size_t r,size_t c):rdim(r),cdim(c)
{
matrix out(r,c,0);
for(size_t i=0;i<v.re_dim();i++)
out[i%r][i/r]=v[i];
(this->m_val)=out.m_val;
}
matrix::matrix(const c&v):rdim(v.re_dim()),cdim(1)
{
matrix out(v.re_dim(),1,0);
for(size_t i=0;i<out.rdim;i++)
out[i][0]=v[i];
this->m_val=out.m_val;
}
//类factor有关定义
factor::factor(const c&v)
{
vector<double>temp=v.const_re_c_val();
for_each(temp.begin(),temp.end(),
[&](const double&d){Levels.insert(d);}
);
vector<double>s_to_v(Levels.begin(),Levels.end());
vector<pair<double,size_t>>ttemp;
for(size_t i=0;i<s_to_v.size();i++)
ttemp.push_back(make_pair(s_to_v[i],i));
for(size_t i=0;i<temp.size();i++)
{
for(size_t j=0;j<ttemp.size();j++)
{
if(temp[i]==ttemp[j].first)
{
f_val.push_back(ttemp[j]);
}
}
}
for(size_t i=0;i<ttemp.size();i++)
f_val_to_levels.insert(ttemp[i]);
auto a=f_val.begin();
for(size_t i=0;i<f_val.size();i++)
{
i_val.push_back(a->second);
a++;
}
}
void factor::out(size_t n)
{
vector<double>temp(f_val.size(),0);
c out0(temp);
auto a=f_val.begin();
for(size_t i=0;i<out0.re_dim();i++)
{
out0[i]=a->first;
a++;
}
vector<double>temp0(Levels.size(),0);
c out1(temp0);
set<double>::iterator b=Levels.begin();
for(size_t i=0;i<out1.re_dim();i++)
{
out1[i]=*b;
b++;
}
out0.out(n);
cout<<"Levels: ";
cout<<endl;
out1.out(n);
}
其他函数即测试程序文件：
use_R.cpp
include"class_of_R.h"
//若干友元函数c operator+(const c&a,const c&b)
{
size_t MAX=b.re_dim();
if(a.re_dim()>b.re_dim())
MAX=a.re_dim();
vector<double>temp0(MAX,0);
vector<double>temp1(MAX,0);
vector<double>temp2(MAX,0);
copy(a.c_val.begin(),a.c_val.end(),temp0.begin());
copy(b.c_val.begin(),b.c_val.end(),temp1.begin());
for(size_t i=0;i<MAX;i++)
temp2[i]=temp0[i]+temp1[i];
return c(temp2);
}
c operator-(const c&a)
{
c out(a);
for(size_t i=0;i<out.re_dim();i++)
out[i]=-out[i];
return out;
}
c operator-(const c&a,const c&b)
{
return (a+(-b));
}
c rep(const double&d,size_t n)
{
return c(vector<double>(n,d));
}
c operator+(const c&a,const double&d)
{
return a+rep(d,a.re_dim());
}
c operator-(const c&a,const double&d)
{
return (a+(-d));
}
c operator+(const double&d,const c&a)
{
return a+d;
}
c operator-(const double&d,const c&a)
{
return (-a)+d;
}
c operator*(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]*=d;
return out;
}
c operator*(const c&a,const double&d)
{
return d*a;
}
c operator/(const double&d,const c&a)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]=d/a[i];
return out;
}
c operator/(const c&a,const double&d)
{
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]=a[i]/d;
return out;
}
c operator*(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error("dim error!\n");
c out(a);
for(size_t i=0;i<a.re_dim();i++)
out[i]*=b[i];
return out;
}
matrix operator+(const matrix&a,const matrix&b)
{
size_t r_MAX=a.rdim,c_MAX=a.cdim;
if(a.rdim<b.rdim)
r_MAX=b.rdim;
if(a.cdim<b.cdim)
c_MAX=b.cdim;
matrix tempa(r_MAX,c_MAX,0);
for(size_t i=0;i<a.rdim;i++)
for(size_t j=0;j<a.cdim;j++)
tempa[i][j]=a[i][j];
matrix tempb(r_MAX,c_MAX,0);
for(size_t i=0;i<b.rdim;i++)
for(size_t j=0;j<b.cdim;j++)
tempb[i][j]=b[i][j];
matrix tempout(r_MAX,c_MAX,0);
for(size_t i=0;i<tempout.rdim;i++)
for(size_t j=0;j<tempout.cdim;j++)
tempout[i][j]=tempa[i][j]+tempb[i][j];
return tempout;
}
matrix operator-(const matrix&a)
{
matrix out(a);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]=-out[i][j];
return out;
}
matrix operator-(const matrix&a,const matrix&b)
{
return ((-b)+a);
}
matrix operator+(const matrix&a,const double&d)
{
return a+matrix(a.rdim,a.cdim,d);
}
matrix operator+(const double&d,const matrix&a)
{
return a+d;
}
matrix operator-(const matrix&a,const double&d)
{
return a+(-d);
}
matrix operator-(const double&b,const matrix&a)
{
return (-a)+b;
}
matrix operator*(const matrix&a,const double&d)
{
matrix out(a);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]*=d;
return out;
}
matrix operator*(const double&d,const matrix&a)
{
return a*d;
}
matrix operator/(const matrix&a,const double&d)
{
return a*(1/d);
}
matrix operator/(const double&d,const matrix&a)
{
matrix out(a.rdim,a.cdim,0);
for(size_t i=0;i<out.rdim;i++)
for(size_t j=0;j<out.cdim;j++)
out[i][j]=d/a[i][j];
return out;
}
c levels(const factor&f)
{
vector<double> temp(f.Levels.begin(),f.Levels.end());
c done(temp);
done.out();
return done;
}
size_t nlevels(const factor&f)
{
size_t out=f.Levels.size();
cout<<out<<endl;
return out;
}
//返回输入因子相应在存储中的下标构成的因子
factor as_integer(const factor&f)
{
vector<double> temp(f.i_val.begin(),f.i_val.end());
factor temp0=factor(c(temp));
for_each(temp0.f_val.begin(),temp0.f_val.end(),
[](const pair<double,size_t>&p)
{
cout<<p.first<<" ";
}
);
cout<<endl;
return temp0;
}
factor as_factor(const c&v)
{
factor temp0=factor(v);
temp0.out();
return temp0;
}
//基本运算函数
matrix diag(const c&v)
{
matrix out(v.re_dim(),v.re_dim(),0);
for(size_t i=0;i<v.re_dim();i++)
out[i][i]=v[i];
return out;
}
matrix diag(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
out[i][i]=m[i][i];
return out;
}
double inter_product(const c&a,const c&b)
{
if(a.re_dim()!=b.re_dim())
throw runtime_error("dim error!\n");
double out=0;
for(size_t i=0;i<a.re_dim();i++)
out+=(a[i]*b[i]);
return out;
}
c operator%(const matrix&m,const c&v)
{
if(m.re_dim().second!=v.re_dim())
throw runtime_error("dim error!\n");
c out(rep(0,m.re_dim().first));
for(size_t i=0;i<out.re_dim();i++)
out[i]=inter_product(m[i],v);
return out;
}
//同下
double prod(const c&v)
{
double out=1;
for(size_t i=0;i<v.re_dim();i++)
out*=v[i];
return out;
}
//对于matrix的求和由c复制构造函数的重载函数定义：c(const matrix&m)
double sum(const c&v)
{
double out=v[0];
for(size_t i=1;i<v.re_dim();i++)
out+=v[i];
return out;
}
matrix operator%(const c&v,const matrix&m)
{
if(v.re_dim()!=m.re_dim().first)
throw runtime_error("dim error!\n");
matrix out(1,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().second;i++)
out[0][i]=inter_product(v,m.col(i));
return out;
}
matrix outer_product(const c&a,const c&b)
{
matrix out(a.re_dim(),b.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=a[i]*b[j];
return out;
}
size_t length(const c&v)
{
return v.re_dim();
}
c sort(const c&v)
{
vector<double>out;
for(size_t i=0;i<v.re_dim();i++)
out.push_back(v[i]);
sort(out.begin(),out.end());
return c(out);
}
double mean(const c&v)
{
return sum(v)/length(v);
}
//其转换版本返回矩阵转为向量的方差
double var(const c&v)
{
return mean((v-mean(v))*(v-mean(v)));
}
c connect(initializer_list<c> lst)
{
vector<double>temp;
for_each(lst.begin(),lst.end(),[&](const c&v){for(auto a:v.const_re_c_val())temp.push_back(a);});
return c(temp);
}
//标准差
template<typename any>
double STD(const any&a)
{
return sqrt(var(a)*(length(a))/(length(a)-1));
}
matrix rbind(initializer_list<c>lst)
{
size_t r=0;
auto a=lst.begin();
while(a!=lst.end())
{
a++;
r++;
}
size_t co=0;
for_each(lst.begin(),lst.end(),[&](const c&v){if(v.re_dim()>co)co=v.re_dim();});
matrix out(r,co,0);
size_t i=0;
a=lst.begin();
while(a!=lst.end())
{
vector<double>temp(co,0);
for(size_t i=0;i<a->const_re_c_val().size();i++)
temp[i]=a->operator;
out[i]=c(temp);
i++;
a++;
}
return out;
}
size_t dim(const c&v)
{
return v.re_dim();
}
pair<size_t,size_t> dim(const matrix&m)
{
return m.re_dim();
}
double min(const c&v)
{
double temp=v[0];
for(size_t i=1;i<length(v);i++)
if(v[i]<temp)
temp=v[i];
return temp;
}
double max(const c&v)
{
double temp=v[0];
for(size_t i=1;i<length(v);i++)
if(v[i]>temp)
temp=v[i];
return temp;
}
c pmin(initializer_list<c>lst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i<temp.re_dim();i++)
{
temp[i]=min(*a);
a++;
}
return temp;
}
c pmax(initializer_list<c>lst)
{
c temp=rep(0,lst.size());
auto a=lst.begin();
for(size_t i=0;i<temp.re_dim();i++)
{
temp[i]=max(*a);
a++;
}
return temp;
}
template<typename any>
pair<double,double> range(const any&v)
{
cout<<"MIN: "<<min(v)<<" MAX: "<<max(v)<<endl;
return make_pair(min(v),max(v));
}
double median(const c&v)
{
return sort(v)[length(v)/2];
}
double var(const c&a,const c&b)
{
if(length(a)!=length(b))
throw runtime_error("dim error!\n");
c t1(a-mean(a)),t2(b-mean(b));
return inter_product(t1,t2)/(length(a)-1);
}
double cov(const c&a,const c&b)
{
return var(a,b);
}
double cor(const c&a,const c&b)
{
return var(a,b)/STD(a)/STD(b);
}
matrix var(const matrix&a,const matrix&b)
{
if(a.re_dim().first!=b.re_dim().first)
throw runtime_error("dim error!\n");
matrix temp(a.re_dim().second,b.re_dim().second,0);
for(size_t i=0;i<temp.re_dim().first;i++)
for(size_t j=0;j<temp.re_dim().second;j++)
temp[i][j]=var(a.col(i),b.col(j));
return temp;
}
matrix cov(const matrix&a,const matrix&b)
{
return var(a,b);
}
matrix cor(const matrix&a,const matrix&b)
{
matrix temp=var(a,b);
for(size_t i=0;i<temp.re_dim().first;i++)
for(size_t j=0;j<temp.re_dim().second;j++)
temp[i][j]/=((STD(a.col(i))*STD(b.col(j))));
return temp;
}
matrix cov(const matrix&a)
{
return var(a,a);
}
matrix cor(const matrix&a)
{
return cor(a,a);
}
//由于关联容器的删除方法利用误差排序的方法
//利用顺序容器迭代器的删除方法有困难
c rank(const c&v)
{
vector<double>temp00(v.const_re_c_val());
for(size_t i=0;i<temp00.size();i++)
temp00[i]+=(0.00000001*i);
vector<double>temp0(temp00);
multiset<double>temp2(temp0.begin(),temp0.end());
vector<double>temp1(v.re_dim());
size_t i=0,j=0;
double pre_val=min(c(vector<double>(temp2.begin(),temp2.end())));
vector<double>::iterator s0;
while(i<v.re_dim())
{
double val=min(c(vector<double>(temp2.begin(),temp2.end())));
s0=find(temp0.begin(),temp0.end(),val);
if(sqrt((pre_val-val)*(pre_val-val))>0.000001)
j++;
temp1[s0-temp0.begin()]=j;
temp2.erase(val);
i++;
pre_val=val;
}
return c(temp1);
}
c rev(const c&v)
{
vector<double> temp;
auto a=v.const_re_c_val().rbegin();
while(a!=v.const_re_c_val().rend())
{
temp.push_back(*a);
a++;
}
return c(temp);
}
//向量先排序再输出排序前的下标值
c order(const c&v)
{
unordered_multimap<double,double>m0;
for(size_t i=0;i<length(v);i++)
m0.insert(make_pair(v[i],i));
multimap<double,double>m1(m0.begin(),m0.end());
c temp(rep(0,length(v)));
auto a=m1.begin();
size_t i=0;
while(a!=m1.end())
{
temp[i]=a->second;
i++;
a++;
}
return temp;
}
c cumsum(const c&v)
{
c temp(v);
for(size_t i=1;i<dim(v);i++)
temp[i]+=temp[i-1];
return temp;
}
c cumprod(const c&v)
{
c temp(v);
for(size_t i=1;i<length(v);i++)
temp[i]*=temp[i-1];
return temp;
}
matrix cumsum(const matrix&m)
{
c temp(m);
for(size_t i=1;i<length(m);i++)
temp[i]+=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
matrix cumprod(const matrix&m)
{
c temp(m);
for(size_t i=1;i<length(m);i++)
temp[i]*=temp[i-1];
return matrix(temp,m.re_dim().first,m.re_dim().second);
}
//矩阵运算函数
//初等行变换化上三角求解行列式
double det(const matrix&m)
{
auto f=[](const matrix&mm)->pair<matrix,double>
{
if(mm.re_dim().first!=mm.re_dim().second)
throw runtime_error("dim error!\n");
size_t MIN=mm.re_dim().first;
double mult=1;
size_t col=0;
size_t bar=0;
size_t b_c=0;
matrix t(mm);
while(b_c<MIN)
{
for(size_t i=bar+1;i<MIN;)
{
if(t[i][col]==0)
i++;
else
{
t.Swap(i,bar);
mult*=(-1);
break;
}
}
if(bar!=MIN+1)
{
if(t[bar][col]!=0)
{
for(size_t i=bar+1;i<MIN;i++)
{
t[bar][i]=t[bar][i]/t[bar][col];
if(i==bar+1)
mult*=t[bar][col];
}
if(col!=MIN-1)
t[bar][col]=1;
}
for(size_t i=bar+1;i<MIN;i++)
{
t[i]=t[i]-t[i][col]*t[bar];
}
}
col++;
bar++;
b_c=bar;
if(col>bar);
b_c=col;
}
return make_pair(t,mult);
};
pair<matrix,double> tt=f(m);
return tt.second*prod(diag(tt.first)%rep(1,m.re_dim().first));
}
//Cramer法则解线性方程组：要求自变量矩阵为方阵且相应行列式不为0
c solve(const matrix&m,const c&v)
{
c temp(rep(0,v.re_dim()));
double low=det(m);
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
if(low==0)
throw runtime_error("singular!\n");
for(size_t i=0;i<v.re_dim();i++)
{
matrix temp0(m);
temp0.col_change(i,v);
temp[i]=det(temp0)/low;
}
return temp;
}
double tr(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
return inter_product(diag(m)%rep(1,m.re_dim().first),rep(1,m.re_dim().first));
}
matrix T(const matrix&m)
{
matrix out(m.re_dim().second,m.re_dim().first,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=m[j][i];
return out;
}
matrix cbind(initializer_list<c>lst)
{
return T(rbind(lst));
}
matrix operator%(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first,b.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=inter_product(a[i],b.col(j));
return out;
}
double F_norm(const matrix&m)
{
return sqrt(tr(T(m)%m));
}
//a,b为区间 d为步长 e为精度 默认为求整数精度：以1为精度特征值 如所求特征值太少可通过减少d改善
//只能求实值
//可以用F_norm作为a、b的界（-F_norm(m),F_norm(m)）可确保值在范围中
//可用于求解但效率低
c eigen(const matrix&m,const double&a,const double &b,const double &d=0.0001,const double &e=1)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
vector<double>temp;
double pre_beg=100000000;
for(double beg=a;beg<b;beg+=d)
{
matrix mm0(diag(rep(beg,m.re_dim().first))-m);
double s0=det(mm0);
if(sqrt(s0*s0)<e)
{
if(sqrt((pre_beg-beg)*(pre_beg-beg))>1)
{
temp.push_back(beg);
pre_beg=beg;
cout<<"Done!\n";
cout<<beg<<endl;//eigenvalue
}
}
}
return c(temp);
}
//重载求eign的函数 d:步长 可保证绝对求出所有特征值（利用循环增大精度）
c eigen(const matrix&m,const double&d=0.1)
{
double dd=d;
double ra=F_norm(m);
double trr=tr(m);
c temp;
double r=0;
do
{
dd*=0.1;//原来设为0.1
temp=eigen(m,-ra,ra,dd);
if(temp.re_dim()==m.re_dim().first)
break;
r=sum(temp)-trr;
}
while(sqrt((r*r)>1));
//与迹的差距（精度）设定为1 绝对值小于此值的特征值被忽略
//此种方法的一个缺陷为异号特征值的抵消 对于同号特征值有把握
//只适用于所有特征值为实值的情况
//如有复特征值，只输出实数的并循环不能结束
return temp;
}
//矩阵的余子阵
matrix M(const matrix&m,size_t k,size_t l)
{
matrix temp(m);
set<size_t> s1;
for(size_t i=0;i<m.re_dim().first;i++)
s1.insert(i);
set<size_t> s2;
for(size_t i=0;i<m.re_dim().second;i++)
s2.insert(i);
s1.erase(k);
s2.erase(l);
vector<double>temp00;
for(auto j:s2)
for(auto i:s1)
temp00.push_back(temp[i][j]);
return matrix (c(temp00),m.re_dim().first-1,m.re_dim().second-1);
}
//方阵的伴随矩阵
matrix A_accompany(const matrix&m)
{
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=det(M(m,i,j))*(((i+j)%2)?-1:1);
//代数余子式
return T(out);
}
//求方阵的逆矩阵（利用伴随矩阵方法）
matrix solve(const matrix&m)
{
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error("singular!\n");
return A_accompany(m)/de;
}
//利用逆矩阵法求解方程组 要求自变量阵为方阵 且可逆
c solve_eq(const matrix&m,const c&v)
{
c out;
double de=det(m);
if(sqrt(de*de)<0.000001)
throw runtime_error("singular!\n");
return solve(m)%v;
}
//化列向量组构成的矩阵为正交化后列向量组构成的矩阵（Schmidt 正交化）
matrix orth_normal(const matrix &m)
{
matrix out(rep(0,m.re_dim().first*m.re_dim().second),m.re_dim().first,m.re_dim().second);
vector<c>dec;
dec.push_back(m.col(0)/sqrt(inter_product(m.col(0),m.col(0))));
for(size_t i=1;i<m.re_dim().second;i++)
{
c temp=m.col(i);
for(size_t j=1;j<=i;j++)
temp=temp-inter_product(dec[j-1],m.col(i))*dec[j-1];
temp=temp/sqrt(inter_product(temp,temp));
dec.push_back(temp);
}
for(size_t i=0;i<m.re_dim().second;i++)
out.col_change(i,dec[i]);
return out;
}
//生成相应位数的由0/1字符构成的字符串
vector<string> generator0_1(size_t n)
{
list<string>out;
list<string>require;
out.push_back("0");
out.push_back("1");
size_t i=0;
while(i<n-1)
{
if(i>0)
out=require;
size_t s0=out.size();
list<string>::iterator pre_begin=out.begin();
auto pre_end=out.end();
size_t k=0;
while(pre_begin!=pre_end)
{
pre_begin++;
k++;
}
pre_begin=out.begin();
for(size_t i=0;i<k;i++)
{
string temp0=*pre_begin+"0";
string temp1=*pre_begin+"1";
out.push_back(temp0);
out.push_back(temp1);
pre_begin++;
}
require.clear();
list<string>::iterator a=out.begin();
for(size_t i=0;i<s0;i++)
a++;
while(a!=out.end())
{
require.push_back(*a);
a++;
}
i++;
}
return vector<string>(require.begin(),require.end());
}
//利用矩阵向量方法生成上述数组
matrix matrix_generator0_1(size_t n)
{
matrix temp0(2,1,0);
matrix temp1;
temp0[0][0]=0;
temp0[1][0]=1;
size_t i=0;
while(i<n-1)
{
if(i>0)
temp0=temp1;
// size_t s0=temp0.re_dim().second;
size_t pre_begin=0;
auto pre_end=temp0.re_dim().first;
temp1=matrix(temp0.re_dim().first*2,temp0.re_dim().second+1,0);
size_t temp1_index=0;
while(pre_begin!=pre_end)
{
vector<double> tt0=temp0[pre_begin].const_re_c_val();
vector<double> tt1=tt0;
tt0.push_back(0);
tt1.push_back(1);
temp1[temp1_index]=c(tt0);
temp1_index++;
temp1[temp1_index]=c(tt1);
temp1_index++;
pre_begin++;
}
i++;
}
return temp1;
}
//将固定位数的下标值与0/1数值构成的map以vector的形式返回
//每个map都是一种组合
vector<map<size_t,size_t>> index_map_generator(size_t n)
{
vector<string>s0(generator0_1(n));
vector<map<size_t,size_t>>outer;
for(size_t i=0;i<s0.size();i++)
{
map<size_t,size_t>tp;
for(size_t j=0;j<s0[i].size();j++)
{
size_t temp=1;
if(s0[i][j]=='0')
temp=0;
tp.insert(make_pair(j,temp));
}
outer.push_back(tp);
}
return outer;
}
//上述产生方式的矩阵版本
vector<map<size_t,size_t>> index_map_matrix_generator(size_t n)
{
matrix temp=matrix_generator0_1(n);
vector<map<size_t,size_t>> out(temp.re_dim().first);
for(size_t i=0;i<out.size();i++)
{
map<size_t,size_t>m;
for(size_t j=0;j<temp.re_dim().second;j++)
m.insert(make_pair(j,temp[i][j]));
out[i]=m;
}
return out;
}
//返回上述index_map_generator(m)中指向0_1和为n的子集
vector<map<size_t,size_t>> require_dim(size_t m,size_t n)
{
vector<map<size_t,size_t>>v0=index_map_generator(m);
vector<map<size_t,size_t>>out;
auto a=v0.begin();
while(a!=v0.end())
{
size_t temp=0;
for_each(a->begin(),a->end(),[&](const pair<size_t,size_t>&p){
if(p.second==1)
temp++;
});
if(temp==n)
out.push_back(*a);
a++;
}
return out;
}
//matrix 所有n阶子式构成的向量
vector<double> seq_of_subdet(const matrix&m,size_t n)
{
vector<double>out;
vector<map<size_t,size_t>>r_map=require_dim(m.re_dim().first,n);
vector<map<size_t,size_t>>c_map=require_dim(m.re_dim().second,n);
for(size_t i=0;i<r_map.size();i++)
for(size_t j=0;j<c_map.size();j++)
{
set<size_t>r_require;
for_each(r_map[i].begin(),r_map[i].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)r_require.insert(p.first);});
set<size_t>c_require;
for_each(c_map[j].begin(),c_map[j].end(),[&](const pair<size_t,size_t>&p){if(p.second==1)c_require.insert(p.first);});
matrix r_matrix(n,m.re_dim().second,0);
auto a=r_require.begin();
for(size_t i=0;i<r_matrix.re_dim().first;i++)
{
r_matrix[i]=m[*a];
a++;
}
matrix c_matrix(r_matrix.re_dim().first,n,0);
auto b=c_require.begin();
for(size_t i=0;i<n;i++)
{
c_matrix.col_change(i,r_matrix.col(*b));
b++;
}
out.push_back(det(c_matrix));
}
return out;
}
//利用子式的方法求出矩阵的秩
size_t rank(const matrix&m)
{
size_t MAX=m.re_dim().first;
if(m.re_dim().second<m.re_dim().first)
MAX=m.re_dim().second;
size_t out=MAX;
size_t pre_num_of_non0=0;
for(size_t i=MAX;i>0;)
{
vector<double>temp=seq_of_subdet(m,i);
auto a=temp.begin();
while(a!=temp.end())
{
a=sqrt((*a)(*a));
a++;
}
size_t num_of_non0=0;
for_each(temp.begin(),temp.end(),[&](const double&d){if(d>0.00000001)num_of_non0++;});
if(pre_num_of_non0==0)
i--;
if(num_of_non0!=0)
return i+1;
pre_num_of_non0=num_of_non0;
}
return out;
}
//回归系数的求解 要求增广后的数据阵（含常系数）可逆
c lm(const c&y,const matrix&m)
{
matrix X=matrix(m.re_dim().first,m.re_dim().second+1,0);
X.col_change(0,rep(1,X.re_dim().first));
for(size_t i=0;i<m.re_dim().second;i++)
X.col_change(i+1,m.col(i));
if(det(T(X)%X)<0.00000001)
throw runtime_error("singular!\n");
c out(solve(T(X)%X)%T(X)%y);
cout<<"The intercept is: "<<out[0]<<endl;
cout<<"Others are:\n";
for(size_t i=1;i<out.re_dim();i++)
cout<<"x"<<i<<" is "<<out[i]<<endl;
cout<<endl;
return out;
}
//抽样函数 参数为向量v 抽取的个数n 是否有放回replace 抽取的概率prob
//使用此方法前必须在函数调用外设定随机数种子
//要求没有0概率点
c sample(const c&v,size_t n,bool replace=1,const c&prob=c())
{
if(!replace&&v.re_dim()<n)
throw out_of_range("n is too large\n");
//抽样概率
vector<double>probb;
if(prob==c())
probb=rep(1/(double)v.re_dim(),v.re_dim()).const_re_c_val();
else
{
probb=prob.const_re_c_val();
double total=0;
for(auto a:probb)
total+=a;
for(auto &a:probb)
a/=total;
}
ifstream fin("h://runif0.txt");
vector<double>runif;
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
//c(probb).out();
c cprobb=cumsum(c(probb));
//cprobb.out();
c require(rep(0,n));
if(replace)
{
size_t i=0;
while(i<n)
{
double uniform=runif[rand()%30000];
//cout<<uniform<<endl;
double conclusion=v[v.re_dim()-1];
for(size_t i=0;cprobb[i]<1;i++)
{
if(uniform<cprobb[i])
{
conclusion=v[i];
break;
}
}
require[i]=conclusion;
i++;
}
}
else
{
size_t i=0;
map<double,double>mdd;
for(size_t i=0;i<v.re_dim();i++)
mdd.insert(make_pair(v[i],c(probb)[i]));
while(i<n)
{
vector<double> temp_v;
vector<double> temp_prob;
for_each(mdd.begin(),mdd.end(),[&](const pair<double,double>&p){
temp_v.push_back(p.first);
temp_prob.push_back(p.second);
});
double d=sample(c(temp_v),1,1,c(temp_prob)).const_re_c_val()[0];
require[i]=d;
mdd.erase(d);
i++;
}
}
return require;
}
//产生几何分布随机数 n 产生个数 p成功概率
c rgeom(const size_t& n,const double& p)
{
vector<double>out;
size_t i=0;
srand(0);
ifstream fin("h://runif0.txt");
vector<double>runif;
double val;
while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
while(i<n)
{
double temp=runif[rand()%30000];
out.push_back(size_t(log(temp)/log(1-p))+1);
i++;
}
return c(out);
}
//possion分布随机数 n 生成个数 d 分布参数
c rpois(const size_t&n,const double&d)
{
size_t i=0;
vector<double>conclusion;
srand(time(0));
ifstream fin("h://runif0.txt");
vector<double>runif;//储存30000个均匀随机数的vector
double val;
&nbs
de9e
p; while(fin>>val)
{
runif.push_back(val);
}
fin.clear();
while(i<n)
{
double p=exp(-1*d);
size_t out=0;
double F=p;
double temp=runif[rand()%30000];
while(temp>=F)
{
p=d/(double)(out+1)*p;
F+=p;
out++;
}
conclusion.emplace_back(out);
i++;
}
//for_each(conclusion.begin(),conclusion.end(),[](const double&s){cout<<s<<" ";});
//cout<<endl;
return c(conclusion);
}
//利用间距生成向量
c seq(const double&from,const double&to,const double&by)
{
if(from>to)
throw runtime_error("from > to\n");
vector<double>out;
auto begin=from;
while(begin<=to)
{
out.push_back(begin);
begin+=by;
}
return c(out);
}
//生成列表向量元素的所有组合并将结果以矩阵形式输出
//前面的generate_grid_ 函数为铺垫而设 expand_grid函数接受c为参数给出输出
matrix generate_grid_0(pair<double,c>p0)
{
matrix out(p0.second.re_dim(),2,0);
for(size_t i=0;i<out.re_dim().first;i++)
{
out[i][0]=p0.first;
out[i][1]=p0.second[i];
}
return out;
}
matrix generate_grid_1(const c&c0,const c&c1)
{
vector<double> temp0=c0.const_re_c_val();
size_t len=temp0.size();
matrix out(len*c1.re_dim(),2,0);
for(size_t i=0;i<len;i++)
{
matrix temp00=generate_grid_0(make_pair(temp0[i],c1));
for(size_t j=0;j<c1.re_dim();j++)
{
out[i*c1.re_dim()+j]=temp00[j];
}
}
return out;
}
matrix generate_grid_2(const matrix&m,const c&v)
{
matrix temp=generate_grid_1(m.col(m.re_dim().second-1),v);
//temp.out();
//system("pause");
matrix out(temp.re_dim().first,m.re_dim().second+1,0);
//out.out();
//system("pause");
matrix add(m.re_dim().first,m.re_dim().second-1,0);
for(size_t i=0;i<add.re_dim().second;i++)
add.col_change(i,m.col(i));
//add.out();
//system("pause");
for(size_t i=0;i<m.re_dim().first;i++)
for(size_t j=0;j<v.re_dim();j++)
{
for(size_t k=0;k<add.re_dim().second;k++)
out[i*v.re_dim()+j][k]=add[i][k];
}
out.col_change(out.re_dim().second-2,temp.col(0));
out.col_change(out.re_dim().second-1,temp.col(1));
return out;
}
//若干向量其元素的所有组合
matrix expand_grid(initializer_list<c>i0)
{
vector<c>temp(i0);
if(temp.size()==1)
return matrix(temp[0]);
matrix beg0=generate_grid_1(temp[0],temp[1]);
if(temp.size()==2)
return beg0;
for(size_t i=2;i<temp.size();i++)
{
//beg0.out();
beg0=generate_grid_2(beg0,temp[i]);
//cout<<"Done!\n";
//beg0.out();
//cout<<"Done!\n";
}
return beg0;
}
//得到向量v的所有n个元素的组合 输出以矩阵排列 每列为一个组合
//调用了求得下标组合的 require_dim()
//当组合项数太多时主张利用转置按行查看
matrix combn(const c&v,size_t n)
{
auto a=require_dim(v.re_dim(),n);
auto b=vector<set<size_t>>(a.size());
auto a0=a.begin();
auto b0=b.begin();
while(a0!=a.end())
{
for_each(a0->begin(),a0->end(),
[&](const pair<size_t,size_t>&p)
{
if(p.second==1)
b0->insert(p.first);
});
a0++;
b0++;
}
auto cc=vector<vector<double>>(b.size());
for(size_t i=0;i<cc.size();i++)
{
cc[i]=vector<double>(b[i].begin(),b[i].end());
}
matrix out(cc[0].size(),cc.size(),0);
for(size_t i=0;i<out.re_dim().second;i++)
{
out.col_change(i,cc[i]);
}
matrix oout(out.re_dim().first,out.re_dim().second,0);
for(size_t i=0;i<oout.re_dim().first;i++)
for(size_t j=0;j<oout.re_dim().second;j++)
oout[i][j]=v[out[i][j]];
return oout;
}
//上述函数的重载版本 返回0:m的相应n组合
matrix combn(size_t m,size_t n)
{
vector<double>temp;
for(size_t i=0;i<m+1;i++)
temp.push_back(i);
return combn(temp,n);
}
//向量相应元素在前面是否出现
//出现显示1
c duplicated(const c&v)
{
//函数a返回向量cc的前n个元素构成的向量
auto a=[](const c&cc,size_t n)->c
{
c out(rep(0,n));
for(size_t i=0;i<n;i++)
out[i]=cc[i];
return out;
};
vector<c>vc0(v.re_dim());
for(size_t i=0;i<vc0.size();i++)
vc0[i]=a(v,i+1);
vector<vector<double>>vs0(vc0.size());
for(size_t i=0;i<vs0.size();i++)
vs0[i]=vector<double>(vc0[i].const_re_c_val());
c out(rep(0,v.re_dim()));
for(size_t i=1;i<v.re_dim();i++)
{
if(find(vs0[i-1].begin(),vs0[i-1].end(),v[i])!=vs0[i-1].end())
out[i]=1;
}
return out;
}
//实现unique的操作看似与factor的levels相同
//实际在序上是不同的
//levels保持double的序（由set有序决定）
//unique应保持原向量的序
//unique由duplicated定义
c unique(const c&v)
{
c d0=duplicated(v);
vector<pair<double,double>>vm0;
for(size_t i=0;i<v.re_dim();i++)
vm0.push_back(make_pair(v[i],d0[i]));
vector<double>out;
for(size_t i=0;i<vm0.size();i++)
{
if(vm0[i].second==0)
out.push_back(vm0[i].first);
}
return c(out);
}
//矩阵叉乘 T(a)%b (向量内积的矩阵版本)
matrix crossprod(const matrix&a,const matrix&b)
{
return T(a)%b;
}
//outer 生成外积 相应的可将外机作用于函数
matrix outer(const c&v1,const c&v2,function<double(const double&,const double&)>p=function<double(const double&,const double&)>())
{
if(!p)
{
matrix out(v1.re_dim(),v2.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=v1[i]*v2[j];
return out;
}
else
{
matrix out(v1.re_dim(),v2.re_dim(),0);
for(size_t i=0;i<out.re_dim().first;i++)
for(size_t j=0;j<out.re_dim().second;j++)
out[i][j]=p(v1[i],v2[j]);
return out;
}
}
//矩阵Kronecker积
matrix kronecker(const matrix&a,const matrix&b)
{
matrix out(a.re_dim().first*b.re_dim().first,a.re_dim().second*b.re_dim().second,0);
vector<vector<matrix>>vvm0(a.re_dim().first,vector<matrix>(a.re_dim().second));
for(size_t i=0;i<a.re_dim().first;i++)
for(size_t j=0;j<a.re_dim().second;j++)
{
vvm0[i][j]=a[i][j]*b;
}
for(size_t i=0;i<vvm0.size();i++)
for(size_t j=0;j<vvm0[0].size();j++)
{
matrix temp=vvm0[i][j];
for(size_t k=0;k<temp.re_dim().first;k++)
for(size_t l=0;l<temp.re_dim().second;l++)
out[i*b.re_dim().first+k][j*b.re_dim().second+l]=temp[k][l];
}
return out;
}
//由方阵上下三角部分组成的0-1阵 不包含对角线元素
matrix lower_tri(const matrix&m)
{
if(m.re_dim().first!=m.re_dim().second)
throw runtime_error("dim error!\n");
matrix out(m.re_dim().first,m.re_dim().second,0);
for(size_t i=1;i<out.re_dim().first;i++)
for(size_t j=0;j<i;j++)
out[i][j]=1;
return out;
}
matrix upper_tri(const matrix&m)
{
return T(lower_tri(m));
}
//矩阵拉直(按列)
c vec(const matrix&m)
{
c out=m.col(0);
for(size_t i=1;i<m.re_dim().second;i++)
out=connect({out,m.col(i)});
return out;
}
//矩阵半拉直 拉直运算忽略主对角线以上部分
c vech(const matrix&m)
{
c out=m.col(0);
for(size_t i=1;i<m.re_dim().second;i++)
{
vector<double> temp=m.col(i).const_re_c_val();
auto a=temp.begin();
for(size_t j=0;j<i;j++)
a++;
vector<double>temp0(a,temp.end());
out=connect({out,c(temp0)});
}
return out;
}
//在方阵满秩情况下 给出方阵A的QR分解
pair<matrix,matrix> qr_of_full_rank(const matrix&A)
{
if(A.re_dim().first!=A.re_dim().second)
throw runtime_error("dim error!\n");
// /*
size_t rk=rank(A);
cout<<"The rank: "<<rk<<endl;
if(rk!=A.re_dim().first)
throw runtime_error("singular!\n");
//*/
cout<<"A=QR"<<endl;
cout<<endl;
cout<<"Q of A:\n";
matrix Q=orth_normal(A);
Q.out();
cout<<"R of m0:\n";
matrix R=T(Q)%A;
R.val_to_double(lower_tri(R),0);
R.out();
// /*
cout<<"The F_norm gap between A and QR:\n";
cout<<F_norm(A-Q%R)<<endl;
//*/
return make_pair(Q,R);
}
//要求自变量矩阵为可逆方阵 利用qr分解 求解线性方程组
//这只是一个形式上的实现 因为求解R逆的方法是用伴随矩阵法
//如能改为初等行变换法应会提高效率
c qr_solve(const matrix&m,const c&v)
{
auto qr=qr_of_full_rank(m);
return solve(qr.second)%T(qr.first)%v;
}
int main()
{
//测试程序：
matrix m2(c({5,0,10,11,0,4,4,3,10,4,3,3,11,3,3,2,10,10,10,10,11,12,13,14,15}),5,5);
m2.out();
c v1({1,1,1,1,2});
c v2=solve(m2,v1);
c v3=qr_solve(m2,v1);
v2.out();
v3.out();
cout<<"The gap of F_norm:\n";
cout<<F_norm(v2-v3)<<endl;
return 0;
}
转载自斯温的博客 http://blog.csdn.net/sinat_30665603
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- 模仿R语言c++ 向量类c 矩阵类matrix等（持续更新 欢迎指点）
- iphone info.plist 各字段含义总结（持续更新，欢迎指点）
- js学习资料笔记合集，持续学习更新中，欢迎大牛们指点
- C++实现R语言向量化运算（向量类：c 矩阵类：matrix）2015.9.11
- 有意思的记录-C++(持续更新）
- C++面试题整理--持续更新
- C/C++文档注释神器——Doxygen常用知识整理（持续更新）
- C++ 备忘录 ( 持续更新中... )
- 帮朋友招开发工程师，欢迎大家关注，持续更新！
- (转）Eclipse 3.6 插件资源-持续更新ing-更新7次--欢迎有好用插件的朋友推荐
- wince系统运行中如何实现更新flash中的NK.BIN？持续探讨，欢迎关注
- 有意思的记录-C++(持续更新）
- C++使用小结（持续更新中）
- 持续集成C++ 实践记（一） 实现定时更新、编译、发送报告
- .NET面试题大全_持续更新，欢迎大家提意见、拍砖 (持续维护中...)
- 你必须了解的c++的特性(持续更新中)
- C++ 低级错误（持续更新中ing）
- C/C++易错点56条（持续更新中）
- C++实用小知识点---持续更新中
- AVL类模板C++(持续更新)