diff --git a/miscmaths.cc b/miscmaths.cc
index e48ae7e9c3f155bbf2e57ab955b39b1d51c9e340..2913b5f04e718a2f1b1f408e7a032f3abea2305a 100644
--- a/miscmaths.cc
+++ b/miscmaths.cc
@@ -1732,6 +1732,29 @@ ReturnMatrix sum(const Matrix& mat, const int dim)
return res;
}
+ReturnMatrix mean(const Matrix& mat, const RowVector& weights, const int dim) //weights are considered to be in the "direction" of dim and normalised to sum 1
+{
+ Matrix res;
+ if (dim == 1){
+ res = zeros(1,mat.Ncols());
+ for (int mc=1; mc<=mat.Ncols(); mc++) {
+ for (int mr=1; mr<=mat.Nrows(); mr++) {
+ res(1,mc) += weights(mr)*mat(mr,mc);
+ }
+ }
+ }
+ else{
+ res = zeros(mat.Nrows(),1);
+ for (int mr=1; mr<=mat.Nrows(); mr++) {
+ for (int mc=1; mc<=mat.Ncols(); mc++) {
+ res(mr,1) += weights(mc)*mat(mr,mc);
+ }
+ }
+ }
+ res.Release();
+ return res;
+}
+
ReturnMatrix mean(const Matrix& mat, const int dim)
{
Matrix res;
@@ -1758,6 +1781,7 @@ ReturnMatrix mean(const Matrix& mat, const int dim)
return res;
}
+
ReturnMatrix var(const Matrix& mat, const int dim)
{
Matrix res, matmean;
@@ -2038,7 +2062,7 @@ ReturnMatrix remmean(const Matrix& mat, const int dim)
return res;
}
-/* ReturnMatrix cov(const Matrix& mat, const int norm)
+ReturnMatrix oldcov(const Matrix& mat, const int norm)
{
SymmetricMatrix res;
Matrix tmp;
@@ -2052,27 +2076,95 @@ ReturnMatrix remmean(const Matrix& mat, const int dim)
res.Release();
return res;
}
-*/
-ReturnMatrix cov(const Matrix& mat, const int norm)
-{
+ReturnMatrix cov(const Matrix& data, const bool sampleCovariance, int econ)
+{
+ //This assumes vectors are stored using column order in data
SymmetricMatrix res;
- res << ones(mat.Ncols(),mat.Ncols());
+ res << zeros(data.Ncols(),data.Ncols());
+ Matrix meanM(mean(data));
+ int N=data.Nrows();
+ if (sampleCovariance && N>1)
+ N--;
+ if ( econ < 1 )
+ econ=data.Nrows();
+ for(int startRow=1; startRow <= data.Nrows(); startRow+=econ) {
+ Matrix suba=data.SubMatrix(startRow,Min(startRow+econ-1,data.Nrows()),1,data.Ncols());
+ for (int row=1; row <= suba.Nrows(); row++)
+ suba.Row(row)-=meanM;
+ res << res + suba.t()*suba/N;
+ }
+ res.Release();
+ return res;
+}
- Matrix meanM;
- int N;
- meanM = mean(mat);
- if (norm == 1) {N = mat.Nrows();}
- else {N = mat.Nrows()-1;}
- for (int ctr=1; ctr <= mat.Nrows(); ctr++)
- res << res + (mat.Row(ctr) - meanM ).t() * (mat.Row(ctr) - meanM);
- res = res/N;
+ReturnMatrix cov_r(const Matrix& data, const bool sampleCovariance, int econ)
+{
+ //This assumes vectors are stored using row order in data
+ SymmetricMatrix res;
+ res << zeros(data.Nrows(),data.Nrows());
+ Matrix meanM(mean(data,2));
+ int N=data.Ncols();
+ if (sampleCovariance && N>1)
+ N--;
+ if ( econ < 1 )
+ econ=data.Ncols();
+ for(int startCol=1; startCol <= data.Ncols(); startCol+=econ) {
+ Matrix suba=data.SubMatrix(1,data.Nrows(),startCol,Min(startCol+econ-1,data.Ncols()));
+ for (int col=1; col <= suba.Ncols(); col++)
+ suba.Column(col)-=meanM;
+ res << res + suba*suba.t()/N;
+ }
res.Release();
return res;
}
-ReturnMatrix corrcoef(const Matrix& mat, const int norm)
+ReturnMatrix cov_r(const Matrix& data, const Matrix& weights2, int econ)
+{
+ //This assumes vectors are stored using row order in data, weights are a single "row". No bool flag as biased vs unbiased isn't relevant here
+ RowVector weights=((weights2/weights2.Sum()).AsRow());
+ SymmetricMatrix res;
+ res << zeros(data.Nrows(),data.Nrows());
+ Matrix meanM(mean(data,weights,2));
+
+ int N=1-weights.SumSquare();//As weights.Sum() is equal to 1
+
+ if ( econ < 1 )
+ econ=data.Ncols();
+ for(int startCol=1; startCol <= data.Ncols(); startCol+=econ) {
+ Matrix suba=data.SubMatrix(1,data.Nrows(),startCol,Min(startCol+econ-1,data.Ncols()));
+ for (int col=1; col <= suba.Ncols(); col++) {
+ suba.Column(col)-=meanM;
+ suba.Column(col)*=sqrt(weights(startCol+col-1));
+ }
+ res << res + suba*suba.t()/N;
+ }
+ write_ascii_matrix("data.mat",data);
+ write_ascii_matrix("weights.mat",weights);
+ write_ascii_matrix("nonorm",cov_r(data,false));
+ write_ascii_matrix("old.mat",res);
+ //res.Release();
+ Matrix Data2=data;
+ for(int ctr=1; ctr <= data.Ncols(); ctr++)
+ Data2.Column(ctr)*=weights(ctr);
+
+
+
+ Matrix res2 = oldcov(Data2.t(),1);
+ res << res2;
+ write_ascii_matrix("new.mat",res);
+ exit(1);
+
+
+
+ return res;
+}
+
+
+
+
+ReturnMatrix corrcoef(const Matrix& mat, const bool norm)
{
SymmetricMatrix res;
res = cov(mat,norm);
diff --git a/miscmaths.h b/miscmaths.h
index 02fa63fa363a1159996dc1907ecba64e034efdc3..23ebfc3f1d4618d7fd7486a36db734785950e266 100644
--- a/miscmaths.h
+++ b/miscmaths.h
@@ -232,6 +232,7 @@ namespace MISCMATHS {
void pow_econ(Matrix& mat, const double exp);
ReturnMatrix sum(const Matrix& mat, const int dim = 1);
ReturnMatrix mean(const Matrix& mat, const int dim = 1);
+ ReturnMatrix mean(const Matrix& mat, const RowVector& weights, const int dim=1);
ReturnMatrix var(const Matrix& mat, const int dim = 1);
ReturnMatrix max(const Matrix& mat);
ReturnMatrix max(const Matrix& mat,ColumnVector& index);
@@ -256,8 +257,14 @@ namespace MISCMATHS {
ReturnMatrix remmean(const Matrix& mat, const int dim = 1);
ReturnMatrix stdev(const Matrix& mat, const int dim = 1);
- ReturnMatrix cov(const Matrix& mat, const int norm = 0);
- ReturnMatrix corrcoef(const Matrix& mat, const int norm = 0);
+ ReturnMatrix cov(const Matrix& mat, const bool sampleCovariance = false, const int econ=20000);
+ ReturnMatrix cov_r(const Matrix& mat, const bool sampleCovariance = false, const int econ=20000);
+ ReturnMatrix cov_r(const Matrix& data, const Matrix& weights, int econ=20000);
+
+
+
+ ReturnMatrix oldcov(const Matrix& mat, const bool norm = false);
+ ReturnMatrix corrcoef(const Matrix& mat, const bool norm = false);
void symm_orth(Matrix &Mat);
void powerspectrum(const Matrix &Mat1, Matrix &Result, bool useLog);
void element_mod_n(Matrix& Mat,double n); //represent each element in modulo n (useful for wrapping phases (n=2*M_PI))