#include <xmath.h>

Public Member Functions
void	GetXV (const int RecN, TFltV &VarV) const
double	GetY (const int RecN) const
double	GetSig (const int RecN) const
void	NR_covsrt (TFltVV &CovarVV, const int &Vars, const TIntV &ia, const int &mfit)
void	NR_gaussj (TFltVV &a, const int &n, TFltVV &b, const int &m)
void	NR_lfit ()
	TLinReg ()
	~TLinReg ()
	TLinReg (TSIn &)
void	Save (TSOut &)
TLinReg &	operator= (const TLinReg &)
int	GetRecs () const
int	GetVars () const
double	GetCf (const int &VarN) const
double	GetCfUncer (const int &VarN) const
double	GetCovar (const int &VarN1, const int &VarN2) const
double	GetChiSq () const
void	Wr () const
Static Public Member Functions
static PLinReg	New (const TFltVV &XVV, const TFltV &YV, const TFltV &SigV=TFltV())
static PLinReg	Load (TSIn &SIn)
static double	LinInterp (const double &x1, const double &y1, const double &x2, const double &y2, const double &AtX) _CMPWARN
Public Attributes
TFltVV	XVV
TFltV	YV
TFltV	SigV
int	Recs
int	Vars
TFltVV	CovarVV
TFltV	CfV
double	ChiSq
Private Attributes
TCRef	CRef
Friends
class	TPt< TLinReg >

Detailed Description

Definition at line 343 of file xmath.h.

Constructor & Destructor Documentation

TLinReg::TLinReg ( ) [inline]

Definition at line 362 of file xmath.h.

{}

TLinReg::~TLinReg ( ) [inline]

Definition at line 365 of file xmath.h.

{}

TLinReg::TLinReg ( TSIn & ) [inline]

Definition at line 366 of file xmath.h.

{Fail;}

Member Function Documentation

double TLinReg::GetCf ( const int & VarN ) const [inline]

Definition at line 375 of file xmath.h.

{return CfV[VarN+1];}

double TLinReg::GetCfUncer ( const int & VarN ) const [inline]

Definition at line 376 of file xmath.h.

                                           {
    return sqrt(double(CovarVV.At(VarN+1, VarN+1)));}

double TLinReg::GetChiSq ( ) const [inline]

Definition at line 381 of file xmath.h.

{return ChiSq;}

double TLinReg::GetCovar	(	const int &	VarN1,
		const int &	VarN2
	)		const `[inline]`

Definition at line 378 of file xmath.h.

                                                            {
    return CovarVV.At(VarN1+1, VarN2+1);}

int TLinReg::GetRecs ( ) const [inline]

Definition at line 372 of file xmath.h.

{return Recs;}

double TLinReg::GetSig ( const int RecN ) const [inline]

Definition at line 357 of file xmath.h.

{return SigV[RecN-1];}

int TLinReg::GetVars ( ) const [inline]

Definition at line 373 of file xmath.h.

{return Vars;}

void TLinReg::GetXV	(	const int	RecN,
		TFltV &	VarV
	)		const `[inline]`

Definition at line 352 of file xmath.h.

                                                {
    VarV.Gen(Vars+1);
    for (int VarN=0; VarN<Vars; VarN++){VarV[VarN+1]=XVV.At(RecN-1, VarN);}
  }

double TLinReg::GetY ( const int RecN ) const [inline]

Definition at line 356 of file xmath.h.

{return YV[RecN-1];}

static double TLinReg::LinInterp	(	const double &	x1,
		const double &	y1,
		const double &	x2,
		const double &	y2,
		const double &	AtX
	)		`[inline, static]`

Definition at line 383 of file xmath.h.

                                                                  {
    if (x1 == x2) return (y1+y2)/2.0;
    const double k = (y2 - y1) / (x2 - x1);
    return k*(AtX - x1) + y1;
  }

static PLinReg TLinReg::Load ( TSIn & SIn ) [inline, static]

Definition at line 367 of file xmath.h.

{return new TLinReg(SIn);}

PLinReg TLinReg::New	(	const TFltVV &	XVV,
		const TFltV &	YV,
		const TFltV &	SigV = `TFltV()`
	)		`[static]`

Definition at line 767 of file xmath.cpp.

                                                                            {
  PLinReg LinReg=PLinReg(new TLinReg());
  LinReg->XVV=_XVV;
  LinReg->YV=_YV;
  if (_SigV.Empty()){
    LinReg->SigV.Gen(LinReg->YV.Len());
    LinReg->SigV.PutAll(1);
  } else {
    LinReg->SigV=_SigV;
  }
  LinReg->Recs=LinReg->XVV.GetXDim();
  LinReg->Vars=LinReg->XVV.GetYDim();
  IAssert(LinReg->Recs>0);
  IAssert(LinReg->Vars>0);
  IAssert(LinReg->YV.Len()==LinReg->Recs);
  IAssert(LinReg->SigV.Len()==LinReg->Recs);
  LinReg->CovarVV.Gen(LinReg->Vars+1, LinReg->Vars+1);
  LinReg->CfV.Gen(LinReg->Vars+1);
  LinReg->NR_lfit();
  return LinReg;
}

void TLinReg::NR_covsrt	(	TFltVV &	CovarVV,
		const int &	Vars,
		const TIntV &	ia,
		const int &	mfit
	)

Definition at line 789 of file xmath.cpp.

                                                                    {
  for (int i=mfit+1; i<=Vars; i++){
    for (int j=1; j<=i; j++){
      CovarVV.At(i, j)=0; CovarVV.At(j, i)=0.0;}
  }
  int k=mfit;
  for (int j=Vars; j>=1; j--){
    if (ia[j]!=0){
      for (int i=1; i<=Vars; i++){Swap(CovarVV.At(i, k), CovarVV.At(i, j));}
      {for (int i=1; i<=Vars; i++){Swap(CovarVV.At(k, i), CovarVV.At(j, i));}}
      k--;
    }
  }
}

void TLinReg::NR_gaussj	(	TFltVV &	a,
		const int &	n,
		TFltVV &	b,
		const int &	m
	)

Definition at line 805 of file xmath.cpp.

                                                                       {
  int i, icol, irow=0, j, k, l, ll;
  double big, dum, pivinv;

  TIntV indxc(n+1);
  TIntV indxr(n+1);
  TIntV ipiv(n+1);
  for (j=1; j<=n; j++){ipiv[j]=0;}
  for (i=1; i<=n; i++){
    big=0.0;
    for (j=1; j<=n; j++){
      if (ipiv[j]!=1){
        for (k=1; k<=n; k++){
          if (ipiv[k]==0){
            if (fabs(double(a.At(j, k))) >= big){
              big=fabs(double(a.At(j, k)));
              irow=j;
              icol=k;
            }
          } else
          if (ipiv[k]>1){
            TExcept::Throw("Singular Matrix(1) in Gauss");}
        }
      }
    }
    ipiv[icol]++;
    if (irow != icol){
      for (l=1; l<=n; l++){Swap(a.At(irow, l), a.At(icol, l));}
      for (l=1; l<=m; l++){Swap(b.At(irow, l), b.At(icol, l));}
    }
    indxr[i]=irow;
    indxc[i]=icol;
    if (a.At(icol, icol)==0.0){
      TExcept::Throw("Singular Matrix(1) in Gauss");}
    pivinv=1.0/a.At(icol, icol);
    a.At(icol, icol)=1.0;
    for (l=1; l<=n; l++){a.At(icol, l)=a.At(icol, l)*pivinv;}
    for (l=1; l<=m; l++){b.At(icol, l)=b.At(icol, l)*pivinv;}
    for (ll=1; ll<=n; ll++){
      if (ll != icol){
        dum=a.At(ll, icol);
        a.At(ll, icol)=0.0;
        for (l=1;l<=n;l++){a.At(ll, l)-=a.At(icol, l)*dum;}
        for (l=1;l<=m;l++){b.At(ll, l)-=b.At(icol, l)*dum;}
      }
    }
  }
  for (l=n; l>=1; l--){
    if (indxr[l]!=indxc[l]){
      for (k=1; k<=n; k++){
        Swap(a.At(k, indxr[l]), a.At(k, indxc[l]));}
    }
  }
}

void TLinReg::NR_lfit ( )

Definition at line 860 of file xmath.cpp.

                     {
  int i,j,k,l,m,mfit=0;
  double ym,wt,sum,sig2i;

  TIntV ia(Vars+1); for (i=1; i<=Vars; i++){ia[i]=1;}
  TFltVV beta(Vars+1, 1+1);
  TFltV afunc(Vars+1);
  for (j=1;j<=Vars;j++){
    if (ia[j]!=0){mfit++;}}
  if (mfit==0){TExcept::Throw("No parameters to be fitted in LFit");}
  for (j=1; j<=mfit; j++){
    for (k=1; k<=mfit; k++){CovarVV.At(j, k)=0.0;}
    beta.At(j, 1)=0.0;
  }
  for (i=1; i<=Recs; i++){
    GetXV(i, afunc); // funcs(XVV[i],afunc,Vars);
    ym=GetY(i);
    if (mfit<Vars){
      for (j=1;j<=Vars;j++){
        if (ia[j]==0){ym-=CfV[j]*afunc[j];}}
    }
    sig2i=1.0/TMath::Sqr(GetSig(i));
    for (j=0, l=1; l<=Vars; l++){
      if (ia[l]!=0){
        wt=afunc[l]*sig2i;
        for (j++, k=0, m=1; m<=l; m++){
          if (ia[m]!=0){CovarVV.At(j, ++k)+=wt*afunc[m];}
        }
        beta.At(j, 1)+=ym*wt;
      }
    }
  }
  for (j=2; j<=mfit; j++){
    for (k=1; k<j; k++){CovarVV.At(k, j)=CovarVV.At(j, k);}
  }
  NR_gaussj(CovarVV, mfit, beta, 1);
  for (j=0, l=1; l<=Vars; l++){
    if (ia[l]!=0){CfV[l]=beta.At(++j, 1);}
  }
  ChiSq=0.0;
  for (i=1; i<=Recs; i++){
    GetXV(i, afunc); // funcs(XVV[i],afunc,Vars);
    for (sum=0.0, j=1; j<=Vars; j++){sum+=CfV[j]*afunc[j];}
    ChiSq+=TMath::Sqr((GetY(i)-sum)/GetSig(i));
  }
  NR_covsrt(CovarVV, Vars, ia, mfit);
}

TLinReg& TLinReg::operator= ( const TLinReg & ) [inline]

Definition at line 370 of file xmath.h.

{Fail; return *this;}

void TLinReg::Save ( TSOut & ) [inline]

Definition at line 368 of file xmath.h.

{Fail;}

void TLinReg::Wr ( ) const

Definition at line 908 of file xmath.cpp.

                       {
  printf("\n%11s %21s\n","parameter","uncertainty");
  for (int i=0; i<Vars;i++){
    printf("  a[%1d] = %8.6f %12.6f\n",
     i+1, GetCf(i), GetCfUncer(i));
  }
  printf("chi-squared = %12f\n", GetChiSq());

  printf("full covariance matrix\n");
  {for (int i=0;i<Vars;i++){
    for (int j=0;j<Vars;j++){
      printf("%12f", GetCovar(i, j));}
    printf("\n");
  }}
}