Updates to cov and mean - In progress

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);

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))