From 3b2fb6609c87e9c4b57e5c433a3bc2c8395c3067 Mon Sep 17 00:00:00 2001
From: Stephen Smith <steve@fmrib.ox.ac.uk>
Date: Fri, 27 Feb 2004 08:24:21 +0000
Subject: [PATCH] Added ols and SD
---
miscmaths.cc | 52 ++++++++++++++++++++++++++++++++++++++++++++++++++++
miscmaths.h | 12 +++++++++++-
2 files changed, 63 insertions(+), 1 deletion(-)
diff --git a/miscmaths.cc b/miscmaths.cc
index df4ac77..e9f6363 100644
--- a/miscmaths.cc
+++ b/miscmaths.cc
@@ -1407,6 +1407,27 @@ ReturnMatrix neq(const Matrix& mat1,const Matrix& mat2)
return res;
}
+ReturnMatrix SD(const Matrix& mat1,const Matrix& mat2)
+{
+ if((mat1.Nrows() != mat2.Nrows()) ||
+ (mat1.Ncols() != mat2.Ncols()) ){
+ cerr <<"MISCMATHS::SD - matrices are of different dimensions"<<endl;
+ exit(-1);
+ }
+ Matrix ret(mat1.Nrows(),mat1.Ncols());
+ for (int r = 1; r <= mat1.Nrows(); r++) {
+ for (int c =1; c <= mat1.Ncols(); c++) {
+ if( mat2(r,c)==0)
+ ret(r,c)=0;
+ else
+ ret(r,c) = mat1(r,c)/mat2(r,c);
+ }
+ }
+
+ ret.Release();
+ return ret;
+}
+
@@ -1648,6 +1669,37 @@ ReturnMatrix read_vest(string p_fname)
return p_mat;
}
+void ols(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Matrix& varcope){
+ // ols
+ // data is t x v
+ // des is t x ev (design matrix)
+ // tc is cons x ev (contrast matrix)
+ // cope and varcope will be cons x v
+ // but will be resized if they are wrong
+ // hence may be passed in uninitialised
+ // TB 2004
+ if(data.Nrows() != des.Nrows()){
+ cerr <<"MISCMATHS::ols - data and design have different number of time points"<<endl;
+ exit(-1);
+ }
+ if(des.Ncols() != tc.Ncols()){
+ cerr <<"MISCMATHS::ols - design and contrast matrix have different number of EVs"<<endl;
+ exit(-1);
+ }
+
+ Matrix pdes = pinv(des);
+ Matrix prevar=diag(tc*pdes*pdes.t()*tc.t());
+ Matrix R=Identity(des.Nrows())-des*pdes;
+ float tR=R.Trace();
+ Matrix pe=pdes*data;
+ cope=tc*pe;
+ Matrix res=data-des*pe;
+ Matrix sigsq=sum(SP(res,res))/tR;
+ varcope=prevar*sigsq;
+
+
+}
+
int conjgrad(ColumnVector& x, const Matrix& A, const ColumnVector& b, int maxit,
float reltol)
{
diff --git a/miscmaths.h b/miscmaths.h
index 557ce9d..5407049 100644
--- a/miscmaths.h
+++ b/miscmaths.h
@@ -186,7 +186,8 @@ namespace MISCMATHS {
ReturnMatrix leqt(const Matrix& mat1,const Matrix& mat2);
ReturnMatrix eq(const Matrix& mat1,const Matrix& mat2);
ReturnMatrix neq(const Matrix& mat1,const Matrix& mat2);
-
+ ReturnMatrix SD(const Matrix& mat1,const Matrix& mat2); // Schur (element-wise) divide
+
void remmean(const Matrix& mat, Matrix& demeanedmat, Matrix& Mean, const int dim = 1);
ReturnMatrix remmean(const Matrix& mat, const int dim = 1);
ReturnMatrix stdev(const Matrix& mat, const int dim = 1);
@@ -194,6 +195,15 @@ namespace MISCMATHS {
ReturnMatrix corrcoef(const Matrix& mat, const int norm = 0);
void symm_orth(Matrix &Mat);
void powerspectrum(const Matrix &Mat1, Matrix &Result, bool useLog);
+
+ // ols
+ // data is t x v
+ // des is t x ev (design matrix)
+ // tc is cons x ev (contrast matrix)
+ // cope and varcope will be cons x v
+ // but will be resized if they are wrong
+ void ols(const Matrix& data,const Matrix& des,const Matrix& tc, Matrix& cope,Matrix& varcope);
+
// Conjugate Gradient methods to solve for x in: A * x = b
// A must be symmetric and positive definite
--
GitLab