-
Christian Beckmann authoredChristian Beckmann authored
melhlprfns.cc 24.78 KiB
/* MELODIC - Multivariate exploratory linear optimized decomposition into
independent components
melhlprfns.cc - misc functions
Christian F. Beckmann, FMRIB Image Analysis Group
Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
#include "melhlprfns.h"
#include "libprob.h"
#include "miscmaths/miscprob.h"
#include "miscmaths/t2z.h"
#include "miscmaths/f2z.h"
namespace Melodic{
void update_mask(volume<float>& mask, Matrix& Data)
{
Matrix DStDev=stdev(Data);
volume4D<float> tmpMask, RawData;
tmpMask.setmatrix(DStDev,mask);
RawData.setmatrix(Data,mask);
float tMmax;
volume<float> tmpMask2;
tmpMask2 = tmpMask[0];
tMmax = tmpMask2.max();
double st_mean = DStDev.Sum()/DStDev.Ncols();
double st_std = stdev(DStDev.t()).AsScalar();
mask = binarise(tmpMask2,(float) max((float) st_mean-3*st_std,(float) 0.01*st_mean),tMmax);
Data = RawData.matrix(mask);
}
void del_vols(volume4D<float>& in, int howmany)
{
for(int ctr=1; ctr<=howmany; ctr++){
in.deletevolume(ctr);
}
}
Matrix calc_FFT(const Matrix& Mat, const bool logpwr)
{
Matrix res;
for(int ctr=1; ctr <= Mat.Ncols(); ctr++){
ColumnVector tmpCol;
tmpCol=Mat.Column(ctr);
ColumnVector FtmpCol_real;
ColumnVector FtmpCol_imag;
ColumnVector tmpPow;
if(tmpCol.Nrows()%2 != 0){
Matrix empty(1,1); empty=0;
tmpCol &= empty;
}
RealFFT(tmpCol,FtmpCol_real,FtmpCol_imag);
tmpPow = pow(FtmpCol_real,2)+pow(FtmpCol_imag,2);
tmpPow = tmpPow.Rows(2,tmpPow.Nrows());
if(logpwr) tmpPow = log(tmpPow);
if(res.Storage()==0){res= tmpPow;}else{res|=tmpPow;}
}
return res;
} //Matrix calc_FFT()
Matrix smoothColumns(const Matrix& inp)
{ Matrix temp(inp);
int ctr1 = temp.Nrows();
Matrix temp2(temp);
temp2=0;
temp = temp.Row(4) & temp.Row(3) & temp.Row(2) & temp & temp.Row(ctr1-1)
& temp.Row(ctr1-2) &temp.Row(ctr1-3);
double kern[] ={0.0045 , 0.055, 0.25, 0.4, 0.25, 0.055, 0.0045};
double fac = 0.9090909;
for(int cc=1;cc<=temp2.Ncols();cc++){
for(int cr=1;cr<=temp2.Nrows();cr++){
temp2(cr,cc) = fac*( kern[0] * temp(cr,cc) + kern[1] * temp(cr+1,cc) +
kern[2] * temp(cr+2,cc) + kern[3] * temp(cr+3,cc) +
kern[4] * temp(cr+4,cc) + kern[5] * temp(cr+5,cc) +
kern[6] * temp(cr+6,cc));
}
}
return temp2;
} //Matrix smoothColumns()
Matrix convert_to_pbsc(Matrix& inp)
{
Matrix meanimg;
meanimg = mean(inp);
float eps = 0.00001;
for(int ctr=1; ctr <= inp.Ncols(); ctr++){
if(meanimg(1,ctr) < eps)
meanimg(1,ctr) = eps;
}
for(int ctr=1; ctr <= inp.Nrows(); ctr++){
Matrix tmp;
tmp << inp.Row(ctr);
inp.Row(ctr) << 100 * SP((tmp - meanimg),pow(meanimg,-1));
}
inp = remmean(inp);
return meanimg;
} //void convert_to_pbsc
RowVector varnorm(Matrix& in, int dim, float level)
{
Matrix Corr;
Corr = calc_corr(in);
RowVector out;
out = varnorm(in,Corr,dim,level);
return out;
} //RowVector varnorm
void varnorm(Matrix& in, const RowVector& vars)
{
Matrix tmp = vars;
in = SP(in,pow(ones(in.Nrows(),1) * tmp,-1));
}
RowVector varnorm(Matrix& in, Matrix& Corr, int dim, float level)
{
Matrix tmpE, white, dewhite;
RowVector tmpD, tmpD2;
std_pca(remmean(in,2), Corr, tmpE, tmpD);
calc_white(tmpE,tmpD, dim, white, dewhite);
Matrix ws = white * in;
for(int ctr1 = 1; ctr1<=ws.Ncols(); ctr1++)
for(int ctr2 = 1; ctr2<=ws.Nrows(); ctr2++) if(std::abs(ws(ctr2,ctr1)) < level)
ws(ctr2,ctr1)=0;
tmpD = stdev(in - dewhite*ws);
for(int ctr1 = 1; ctr1<=tmpD.Ncols(); ctr1++)
if(tmpD(ctr1) < 0.01){
tmpD(ctr1) = 1.0;
in.Column(ctr1) = 0.0*in.Column(ctr1);
}
varnorm(in,tmpD);
return tmpD;
} //RowVector varnorm
Matrix SP2(const Matrix& in, const Matrix& weights, bool econ)
{
Matrix Res;
Res = in;
if(econ){
ColumnVector tmp;
for(int ctr=1; ctr <= in.Ncols(); ctr++){
tmp = in.Column(ctr);
tmp = tmp * weights(1,ctr);
Res.Column(ctr) = tmp;
}
}
else{
Res = ones(in.Nrows(),1)*weights.Row(1);
Res = SP(in,Res);
}
return Res;
} //Matrix SP
Matrix corrcoef(const Matrix& in1, const Matrix& in2){
Matrix tmp = in1;
tmp |= in2;
Matrix out;
out=MISCMATHS::corrcoef(tmp,0);
return out.SubMatrix(1,in1.Ncols(),in1.Ncols()+1,out.Ncols());
}
Matrix calc_corr(const Matrix& in, bool econ)
{
Matrix Res;
Res = zeros(in.Nrows(),in.Nrows());
if(econ){
ColumnVector tmp;
for(int ctr=1; ctr <= in.Ncols(); ctr++){
tmp = in.Column(ctr);
tmp = tmp - mean(tmp).AsScalar();
Res += (tmp * tmp.t()) / in.Ncols();
}
}
else
Res = cov(in.t());
return Res;
} //Matrix calc_corr
Matrix calc_corr(const Matrix& in, const Matrix& weights, bool econ)
{
Matrix Res;
Res = zeros(in.Nrows(),in.Nrows());
Matrix localweights;
if(weights.Storage() == 0)
localweights = ones(1,in.Ncols());
else
localweights = weights;
if(econ){
ColumnVector tmp;
for(int ctr=1; ctr <= in.Ncols(); ctr++){
tmp = in.Column(ctr);
tmp = tmp - mean(tmp).AsScalar(); tmp = tmp * localweights(1,ctr);
Res += (tmp * tmp.t()) / in.Ncols();
}
}
else{
Res = SP2(in,localweights);
Res = calc_corr(Res, 0);
}
return Res;
} //Matrix calc_corr
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite)
{
// tmpD2= tmpD | tmpPD.AsRow().Columns(tmpPE.Ncols()-param.Ncols()+1,tmpPE.Ncols());
// cerr << tmpPD.AsRow().Columns(tmpPE.Ncols()-param.Ncols()+1,tmpPE.Ncols()) << endl;
//
// Matrix tmpPE;
// tmpPE = SP(param,ones(param.Nrows(),1)*pow(stdev(param,1)*std::sqrt((float)param.Nrows()),-1));
// RE |= tmpPE;
// RowVector tmpD2;
// tmpD2 = tmpD | stdev(param,1).AsRow()*std::sqrt((float)param.Nrows());
// RD << abs(diag(tmpD2.t()));
// RD << RD.SymSubMatrix(N-dim+1,N);
Matrix RE;
DiagonalMatrix RD;
int N = tmpE.Ncols();
dim = std::min(dim,N);
// cerr << stdev(param) << endl;
RE = tmpE.Columns(std::min(N-dim+1+param.Ncols(),N-2),N);
RE |= param;
// cerr << paramS.Nrows() << " x " << paramS.Ncols() << endl;
// cerr << paramS << endl;
RowVector tmpD2;
tmpD2 = tmpD | pow(paramS,2).AsRow();
RD << abs(diag(tmpD2.t()));
// cerr << " " <<tmpD2.Ncols() << " " << N << " " << dim << endl;
RD << RD.SymSubMatrix(N-dim+1+param.Ncols(),N+param.Ncols());
float res = 1.0;
white = sqrt(abs(RD.i()))*RE.t();
dewhite = RE*sqrt(RD);
if(dim < N)
res = PercEV(dim);
return res;
} //Matrix calc_white
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& white, Matrix& dewhite)
{
Matrix RE;
DiagonalMatrix RD;
int N = tmpE.Ncols();
dim = std::min(dim,N);
RE = tmpE.Columns(N-dim+1,N);
RD << abs(diag(tmpD.t()));
RD << RD.SymSubMatrix(N-dim+1,N);
float res = 1.0;
white = sqrt(abs(RD.i()))*RE.t();
dewhite = RE*sqrt(RD);
if(dim < N) res = PercEV(dim);
return res;
} //Matrix calc_white
void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& white, Matrix& dewhite)
{
RowVector tmp(tmpE.Ncols());
float tmp2;
tmp2 = calc_white(tmpE,tmpD, tmp, dim, white, dewhite);
} //Matrix calc_white
void calc_white(const Matrix& tmpE, const RowVector& tmpD, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite)
{
RowVector tmp(tmpE.Ncols());
float tmp2;
tmp2 = calc_white(tmpE,tmpD, tmp, dim, param, paramS, white, dewhite);
} //Matrix calc_white
void calc_white(const Matrix& Corr, int dim, Matrix& white, Matrix& dewhite)
{
Matrix RE;
DiagonalMatrix RD;
SymmetricMatrix tmp;
RowVector tmp2;
tmp << Corr;
EigenValues(tmp,RD,RE);
tmp2 = diag(RD).t();
calc_white(RE,tmp2, dim, white, dewhite);
} //Matrix calc_white
void std_pca(const Matrix& Mat, const Matrix& weights, Matrix& Corr, Matrix& evecs, RowVector& evals)
{
if(weights.Storage()>0)
Corr << calc_corr(Mat, weights);
else
Corr << calc_corr(Mat);
SymmetricMatrix tmp;
tmp << Corr;
DiagonalMatrix tmpD;
EigenValues(tmp,tmpD,evecs);
evals = tmpD.AsRow();
} //void std_pca
void std_pca(const Matrix& Mat, Matrix& Corr, Matrix& evecs, RowVector& evals)
{
Matrix weights;
std_pca(Mat,weights,Corr,evecs,evals);
} //void std_pca
void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc, int iter)
{
Matrix guess;
guess = normrnd(Mat.Nrows(),num_pc);
em_pca(Mat,guess,evecs,evals,num_pc,iter);
} //void em_pca
void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc, int iter)
{
Matrix C;
if(guess.Ncols() < num_pc){
C=normrnd(Mat.Nrows(),num_pc);
C.Columns(1,guess.Ncols()) = guess;
}
else
C = guess;
Matrix tmp, tmp2;
for(int ctr=1; ctr <= iter; ctr++){
// E-Step tmp = C.t()*C;
tmp = tmp.i();
tmp = tmp * C.t();
tmp = tmp * Mat;
// M-Step
tmp2 = tmp * tmp.t();
tmp2 = tmp2.i();
tmp2 = Mat*tmp.t()*tmp2;
C = tmp2;
}
symm_orth(C);
Matrix Evc, tmpC;
RowVector Evl;
tmp = C.t() * Mat;
std_pca(tmp,tmpC,Evc,Evl);
evals = Evl;
evecs = C * Evc;
} //void em_pca
float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim)
{
Matrix Corr, Evecs, tmpWM, tmpDWM, tmp;
RowVector Evals;
std_pca(Mat.t(), Corr, Evecs, Evals);
calc_white(Corr, dim, tmpWM, tmpDWM);
tmp = tmpWM * Mat.t();
cols = tmp.t();
rows << tmpDWM;
float res;
Evals=fliplr(Evals);
res = sum(Evals.Columns(1,dim),2).AsScalar()/sum(Evals,2).AsScalar()*100;
return res;
} // rankapprox
RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows)
{
Matrix out; RowVector res(Mat.Ncols());
for(int ctr1 = 1; ctr1 <= Mat.Ncols(); ctr1++){
Matrix tmpVals(cols.Nrows(),rows.Nrows());
for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() *
(ctr2 - 1) + 1,cols.Nrows()*ctr2 ,ctr1,ctr1);
Matrix tmpcols, tmprows;
res(ctr1) =rankapprox(tmpVals, tmpcols, tmprows);
cols.Column(ctr1) = tmpcols;
rows.Column(ctr1) = tmprows;
}
return res;
} // krfact
RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows)
{
RowVector res;
cols = zeros(colnum,Mat.Ncols());
rows = zeros(int(Mat.Nrows() / colnum),Mat.Ncols());
res = krfact(Mat,cols,rows);
return res;
} // krfact
Matrix krprod(const Matrix& cols, const Matrix& rows)
{
Matrix out;
out = zeros(cols.Nrows()*rows.Nrows(),cols.Ncols());
for(int ctr1 = 1; ctr1 <= cols.Ncols(); ctr1++)
for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
{
out.SubMatrix(cols.Nrows()*(ctr2-1)+1,cols.Nrows()*ctr2,ctr1,ctr1) << cols.Column(ctr1) * rows(ctr2,ctr1);
} return out;
} // krprod
Matrix krapprox(const Matrix& Mat, int size_cols, int dim)
{
Matrix out, cols, rows;
out = zeros(Mat.Nrows(), Mat.Ncols());
cols = zeros(size_cols,Mat.Ncols());
rows = zeros(int(Mat.Nrows() / size_cols), Mat.Ncols());
krfact(Mat,cols,rows);
out = krprod(cols, rows);
return out;
} // krapprox
void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2, RowVector& out3, int& out4, int num_vox, float resels)
{
RowVector AdjEV;
AdjEV << in.AsRow();
AdjEV = AdjEV.Columns(3,AdjEV.Ncols());
AdjEV = AdjEV.Reverse();
RowVector CircleLaw;
int NumVox = (int) floor(num_vox/(2.5*resels));
CircleLaw = Feta(int(AdjEV.Ncols()), NumVox);
for(int ctr=1;ctr<=CircleLaw.Ncols(); ctr++){
if(CircleLaw(ctr)<5*10e-10){CircleLaw(ctr) = 5*10e-10;}
}
/* float slope;
slope = CircleLaw.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
int(AdjEV.Ncols()/4)).Sum() /
AdjEV.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
int(AdjEV.Ncols()/4)).Sum();*/
RowVector PercEV(AdjEV);
PercEV = cumsum(AdjEV / sum(AdjEV,2).AsScalar());
AdjEV << SP(AdjEV,pow(CircleLaw.Columns(1,AdjEV.Ncols()),-1));
SortDescending(AdjEV);
int maxEV = 1;
float threshold = 0.98;
for(int ctr_i = 1; ctr_i < PercEV.Ncols(); ctr_i++){
if((PercEV(ctr_i)<threshold)&&(PercEV(ctr_i+1)>=threshold)){maxEV=ctr_i;}
}
if(maxEV<3){maxEV=PercEV.Ncols()/2;}
RowVector NewEV;
Matrix temp1;
temp1 = abs(AdjEV);
NewEV << temp1.SubMatrix(1,1,1,maxEV);
AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar());
out1 = AdjEV;
out2 = PercEV;
out3 = NewEV;
out4 = maxEV;
} //adj_eigspec
void adj_eigspec(const RowVector& in, RowVector& out1, RowVector& out2)
{
RowVector AdjEV, PercEV;
AdjEV = in.Reverse();
SortDescending(AdjEV);
PercEV = cumsum(AdjEV / sum(AdjEV,2).AsScalar());
AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar()); out1 = AdjEV;
out2 = PercEV;
} //adj_eigspec
RowVector Feta(int n1, int n2)
{
float nu = (float) n1/n2;
float bm = pow((1-sqrt(nu)),2.0);
float bp = pow((1+sqrt(nu)),2.0);
float lrange = 0.9*bm;
float urange = 1.1*bp;
RowVector eta(30*n1);
float rangestepsize = (urange - lrange) / eta.Ncols();
for(int ctr_i = 0; ctr_i < eta.Ncols(); ctr_i++){
eta(ctr_i+1) = lrange + rangestepsize * (ctr_i);
}
RowVector teta(10*n1);
teta = 0;
float stepsize = (bp - bm) / teta.Ncols();
for(int ctr_i = 0; ctr_i < teta.Ncols(); ctr_i++){
teta(ctr_i+1) = stepsize*(ctr_i);
}
RowVector feta(teta);
feta = SP(pow(2*M_PI*nu*(teta + bm),-1), pow(SP(teta, bp-bm-teta),0.5));
teta = teta + bm;
RowVector claw(eta);
claw = 0;
for(int ctr_i = 1; ctr_i <= eta.Ncols(); ctr_i++){
double tmpval = 0.0;
for(int ctr_j = 1; ctr_j <= teta.Ncols(); ctr_j++){
if(( double(teta(ctr_j))/double(eta(ctr_i)) )<1)
tmpval += feta(ctr_j);
}
claw(ctr_i) = n1*(1-stepsize*tmpval);
}
RowVector Res(n1); //invert the CDF
for(int ctr_i = 1; ctr_i < eta.Ncols(); ctr_i++){
if(floor(claw(ctr_i))>floor(claw(ctr_i+1))){
Res(int(floor(claw(ctr_i)))) = eta(ctr_i);
}
}
return Res;
} //RowVector Feta
RowVector cumsum(const RowVector& Inp)
{
UpperTriangularMatrix UT(Inp.Ncols());
UT=1.0;
RowVector Res;
Res = Inp * UT;
return Res;
} //RowVector cumsum
int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, Matrix& Corr, Matrix& tmpE, RowVector &tmpD, float resels, string which)
{
std_pca(in,weights,Corr,tmpE,tmpD);
int maxEV = 1;
RowVector NewEV;
adj_eigspec(tmpD.AsRow(),AdjEV,PercEV,NewEV,maxEV,in.Ncols(),resels);
int res; PPCA = ppca_est(NewEV, in.Ncols(),resels);
ColumnVector tmp = ppca_select(PPCA, res, maxEV, which);
PPCA = tmp | PPCA;
return res;
} //int ppca_dim
int ppca_dim(const Matrix& in, const Matrix& weights, Matrix& PPCA, RowVector& AdjEV, RowVector& PercEV, float resels, string which)
{
RowVector tmpD;
Matrix tmpE;
Matrix Corr;
int res = ppca_dim(in, weights, PPCA, AdjEV, PercEV, Corr, tmpE, tmpD, resels, which);
return res;
} //int ppca_dim
int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which)
{
ColumnVector PPCA;
RowVector AdjEV, PercEV;
int res = ppca_dim(in,weights,PPCA,AdjEV,PercEV,resels,which);
return res;
} //int ppca_dim
ColumnVector ppca_select(Matrix& PPCAest, int& dim, int maxEV, string which)
{
RowVector estimators(5);
estimators = 1.0;
for(int ctr=1; ctr<=PPCAest.Ncols(); ctr++){
PPCAest.Column(ctr) = (PPCAest.Column(ctr) -
min(PPCAest.Column(ctr)).AsScalar()) /
( max(PPCAest.Column(ctr)).AsScalar() -
min(PPCAest.Column(ctr)).AsScalar());
}
int ctr_i = 1;
while((ctr_i< PPCAest.Nrows()-1)&&
(PPCAest(ctr_i,2) < PPCAest(ctr_i+1,2))&&(ctr_i<maxEV))
{estimators(1)=ctr_i+1;ctr_i++;}
ctr_i = 1;
while((ctr_i< PPCAest.Nrows()-1)&&
(PPCAest(ctr_i,3) < PPCAest(ctr_i+1,3))&&(ctr_i<maxEV))
{estimators(2)=ctr_i+1;ctr_i++;}
ctr_i = 1;
while((ctr_i< PPCAest.Nrows()-1)&&
(PPCAest(ctr_i,4) < PPCAest(ctr_i+1,4))&&(ctr_i<maxEV))
{estimators(3)=ctr_i+1;ctr_i++;}
ctr_i = 1;
while((ctr_i< PPCAest.Nrows()-1)&&
(PPCAest(ctr_i,5) < PPCAest(ctr_i+1,5))&&(ctr_i<maxEV))
{estimators(4)=ctr_i+1;ctr_i++;}
ctr_i = 1;
while((ctr_i< PPCAest.Nrows()-1)&&
(PPCAest(ctr_i,6) < PPCAest(ctr_i+1,6))&&(ctr_i<maxEV))
{estimators(5)=ctr_i+1;ctr_i++;}
int res = 0;
ColumnVector PPCA;
RowVector PercEV(PPCAest.Column(1).t());
PercEV = cumsum(PercEV / sum(PercEV,2).AsScalar());
if(which == string("aut"))
if(int(estimators(2)) < int(estimators(1)) &&
float(PercEV(int(estimators(2))))>0.8){
res=int(estimators(2));
PPCA << PPCAest.Column(3);
}else{
res = int(estimators(1)); PPCA << PPCAest.Column(2);
}
if(which == string("lap")){
res = int(estimators(1));
PPCA << PPCAest.Column(2);
}
if(which == string("bic")){
res = int(estimators(2));
PPCA << PPCAest.Column(3);
}
if(which == string("mdl")){
res = int(estimators(3));
PPCA << PPCAest.Column(4);
}
if(which == string("rrn")){
res = int(estimators(4));
PPCA << PPCAest.Column(5);
}
if(which == string("aic")){
res = int(estimators(5));
PPCA << PPCAest.Column(6);
}
if(which == string("median")){
RowVector unsorted(estimators);
SortAscending(unsorted);
ctr_i=1;
res=int(unsorted(3));
while(res != int(estimators(ctr_i)))
ctr_i++;
PPCA << PPCAest.Column(ctr_i);
}
if(res==0 || which == string("mean")){//median estimator
PPCA = mean(PPCAest.Columns(2,6),2);
res=int(mean(estimators,2).AsScalar());
}
dim = res;
return PPCA;
} //RowVector ppca_select
Matrix ppca_est(const RowVector& eigenvalues, const int N1, const float N2)
{
Matrix Res;
Res = ppca_est(eigenvalues, (int) floor(N1/(2.5*N2)));
return Res;
} //Matrix ppca_est
Matrix ppca_est(const RowVector& eigenvalues, const int N)
{
RowVector logLambda(eigenvalues);
logLambda = log(eigenvalues);
int d = logLambda.Ncols();
RowVector k(d);
for(int ctr_i = 1; ctr_i <=d; ctr_i++){
k(ctr_i)=ctr_i;
}
RowVector m(d);
m=d*k-0.5*SP(k,k+1);
RowVector loggam(d);
loggam = 0.5*k.Reverse();
for(int ctr_i = 1; ctr_i <=d; ctr_i++){
loggam(ctr_i)=lgam(loggam(ctr_i));
}
loggam = cumsum(loggam);
RowVector l_probU(d); l_probU = -log(2)*k + loggam - cumsum(0.5*log(M_PI)*k.Reverse());
RowVector tmp1;
tmp1 = -cumsum(eigenvalues).Reverse()+sum(eigenvalues,2).AsScalar();
tmp1(1) = 0.95*tmp1(2);
tmp1=tmp1.Reverse();
RowVector tmp2;
tmp2 = -cumsum(logLambda).Reverse()+sum(logLambda,2).AsScalar();
tmp2(1)=tmp2(2);
tmp2=tmp2.Reverse();
RowVector tmp3;
tmp3 = d-k;
tmp3(d) = 1.0;
RowVector tmp4;
tmp4 = SP(tmp1,pow(tmp3,-1));
for(int ctr_i = 1; ctr_i <=d; ctr_i++){
if(tmp4(ctr_i)<0.01){tmp4(ctr_i)=0.01;}
if(tmp3(ctr_i)<0.01){tmp3(ctr_i)=0.01;}
if(tmp1(ctr_i)<0.01){tmp1(ctr_i)=0.01;}
}
RowVector l_nu;
l_nu = SP(-N/2*(d-k),log(tmp4));
l_nu(d) = 0;
RowVector l_lam;
l_lam = -(N/2)*cumsum(logLambda);
RowVector l_lhood;
l_lhood = SP(pow(tmp3,-1),tmp2)-log(SP(pow(tmp3,-1),tmp1));
Matrix t1,t2, t3;
UpperTriangularMatrix triu(d);
triu = 1.0;
for(int ctr_i = 1; ctr_i <= triu.Ncols(); ctr_i++){
triu(ctr_i,ctr_i)=0.0;
}
t1 = (ones(d,1) * eigenvalues);
t1 = SP(triu,t1.t() - t1);
t2 = pow(tmp4,-1).t()*ones(1,d);
t3 = ones(d,1)*pow(eigenvalues,-1);
t2 = SP(triu, t2.t()-t3.t());
for(int ctr_i = 1; ctr_i <= t1.Ncols(); ctr_i++){
for(int ctr_j = 1; ctr_j <= t1.Nrows(); ctr_j++){
if(t1(ctr_j,ctr_i)<=0){t1(ctr_j,ctr_i)=1;}
}
}
for(int ctr_i = 1; ctr_i <= t2.Ncols(); ctr_i++){
for(int ctr_j = 1; ctr_j <= t2.Nrows(); ctr_j++){
if(t2(ctr_j,ctr_i)<=0){t2(ctr_j,ctr_i)=1;}
}
}
t1 = cumsum(sum(log(t1),2).AsRow());
t2 = cumsum(sum(log(t2),2).AsRow());
RowVector l_Az(d);
l_Az << (t1+t2);
RowVector l_lap;
l_lap = l_probU + l_nu +l_Az + l_lam + 0.5*log(2*M_PI)*(m+k)-0.5*log(N)*k;
RowVector l_BIC;
l_BIC = l_lam + l_nu - 0.5*log(N)*(m+k);
RowVector l_RRN;
l_RRN = -0.5*N*SP(k,log(SP(cumsum(eigenvalues),pow(k,-1))))+l_nu;
RowVector l_AIC;
l_AIC = -2*N*SP(tmp3,l_lhood)+ 2*(1+d*k+0.5*(k-1));
l_AIC = -l_AIC;
RowVector l_MDL;
l_MDL = -N*SP(tmp3,l_lhood)+ 0.5*(1+d*k+0.5*(k-1))*log(N);
l_MDL = -l_MDL;
Matrix Res;
Res = eigenvalues.t();
Res |= l_lap.t();
Res |= l_BIC.t();
Res |= l_MDL.t();
Res |= l_RRN.t();
Res |= l_AIC.t();
return Res;
} //Matrix ppca_est
ColumnVector acf(const ColumnVector& in, int order)
{
double meanval;
meanval = mean(in).AsScalar();
int tpoints = in.Nrows();
ColumnVector y, res;
Matrix X, tmp;
y = in.Rows(order+1,tpoints) - meanval;
X = zeros(order + 1, order);
for(int ctr1 = 1; ctr1 <= order; ctr1++)
X.Column(ctr1) = in.Rows(order + 1 - ctr1, tpoints - ctr1) - meanval;
tmp = X.t()*X;
tmp = tmp.i();
tmp = tmp * X.t();
res << tmp * y;
return res;
} //ColumnVector acf
ColumnVector pacf(const ColumnVector& in, int maxorder)
{
int tpoint = in.Nrows();
ColumnVector res;
res = acf(in, maxorder);
for(int ctr1 = 1; ctr1 <= maxorder; ctr1++)
if ( res.Column(ctr1).AsScalar() < (1/tpoint) + 2/(float)std::pow(tpoint,0.5))
res.Column(ctr1) = 0;
return res;
} //ColumnVector pacf
Matrix est_ar(const Matrix& Mat, int maxorder)
{
Matrix res;
res = zeros(maxorder, Mat.Ncols());
ColumnVector tmp;
for (int ctr = 1; ctr <= Mat.Ncols(); ctr++){
tmp = pacf(Mat.Column(ctr));
res.Column(ctr) = tmp;
}
return res;
} //Matrix est_ar
ColumnVector gen_ar(const ColumnVector& in, int maxorder)
{
float sdev;
sdev = stdev(in).AsScalar();
ColumnVector x, arcoeff, scaled;
scaled = in / sdev; arcoeff = pacf( scaled, maxorder);
x = normrnd(in.Nrows(),1).AsColumn() * sdev;
for(int ctr1=2; ctr1 <= in.Nrows(); ctr1++)
for(int ctr2 = 1; ctr2 <= maxorder; ctr2++)
x(ctr1) = arcoeff(ctr2) * x(std::max(1,int(ctr1-ctr2))) + x(ctr1);
return x;
} //ColumnVector gen_ar
Matrix gen_ar(const Matrix& in, int maxorder)
{
Matrix res;
res = in;
ColumnVector tmp;
for(int ctr=1; ctr <= in.Ncols(); ctr++){
tmp = in.Column(ctr);
res.Column(ctr) = gen_ar(tmp, maxorder);
}
return res;
} //Matrix gen_ar
Matrix gen_arCorr(const Matrix& in, int maxorder)
{
Matrix res;
res = zeros(in.Nrows(), in.Nrows());
ColumnVector tmp;
for(int ctr=1; ctr<= in.Ncols(); ctr++){
tmp = in.Column(ctr);
tmp = gen_ar(tmp, maxorder);
res += tmp * tmp.t();
}
return res;
} //Matrix gen_arCorr
void basicGLM::olsfit(const Matrix& data, const Matrix& design,
const Matrix& contrasts, int DOFadjust)
{
beta = zeros(design.Ncols(),1);
residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1);
t = zeros(1); z = zeros(1); p=-1.0*ones(1);
dof = (int)-1; cbeta = -1.0*ones(1);
if(data.Nrows()==design.Nrows()){
Matrix dat = data;
Matrix tmp = design.t()*design;
Matrix pinvdes = tmp.i()*design.t();
beta = pinvdes * dat;
residu = dat - design*beta;
dof = design.Nrows() - design.Ncols()-1;
sigsq = sum(SP(residu,residu))/dof;
float fact = float(dof) / design.Ncols();
f_fmf = SP(sum(SP(design*beta,design*beta)),pow(sum(SP(residu,residu)),-1)) * fact;
pf_fmf = f_fmf.Row(1);
for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),
int(design.Nrows() -1 -design.Ncols()),f_fmf.Column(ctr1).AsScalar());
if(contrasts.Storage()>0 && contrasts.Ncols()==beta.Nrows()){
cbeta = contrasts*beta;
Matrix tmp = contrasts*pinvdes*pinvdes.t()*contrasts.t();
varcb = diag(tmp)*sigsq;
t = SP(cbeta,pow(varcb,-0.5));
z = t; p=t;
for(int ctr1=1;ctr1<=t.Ncols();ctr1++){
ColumnVector tmp = t.Column(ctr1);
T2z::ComputeZStats(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
z.Column(ctr1) << tmp; T2z::ComputePs(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
p.Column(ctr1) << exp(tmp);
}
}
}
}
}