~ubuntu-branches/ubuntu/raring/plink/raring

« back to all changes in this revision

Viewing changes to stats.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Andreas Tille
  • Date: 2009-10-23 13:35:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20091023133502-rmnkrfi314mevnvt
Tags: 1.07-1
* New upstream version
* Removed debian/patches/20_plink-1.06-gcc4.4.patch because it
  is applied upstream
* Standards-Version: 3.8.3 (no changes needed)
* Debhelper 7
* Build-Depends: zlib1g-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
801
801
}
802
802
 
803
803
 
804
 
void svd(matrix_t & u, vector_t &w, matrix_t &v)
 
804
bool svd(matrix_t & u, vector_t &w, matrix_t &v)
805
805
{
806
806
 
807
807
// #ifdef WITH_LAPACK
820
820
  w.resize(c);
821
821
  sizeMatrix(v,c,c);
822
822
 
823
 
  svdcmp(u,w,v); 
 
823
  bool flag = svdcmp(u,w,v); 
 
824
 
 
825
  return flag;
824
826
  
825
827
  // Look for singular values
826
828
//   double wmax = 0;
836
838
// #endif
837
839
}
838
840
 
839
 
vector< vector<double> > svd_inverse(vector< vector<double> > & u )
 
841
vector< vector<double> > svd_inverse(vector< vector<double> > & u , bool & flag )
840
842
{
841
 
 
842
 
  //  const double eps = 1e-12;
 
843
  
843
844
  const double eps = 1e-24; 
844
845
  
845
846
  if (u.size() == 0) 
854
855
  for (int i=0; i<n; i++) 
855
856
    v[i].resize(n,0);
856
857
 
857
 
  svdcmp(u,w,v); 
 
858
  flag = svdcmp(u,w,v); 
858
859
  
859
860
  // Look for singular values
860
861
  double wmax = 0;
870
871
  // u w t(v)
871
872
 
872
873
  // row U * 1/w
 
874
  
 
875
  // results matrix
873
876
  vector<vector<double> > r(n);
874
877
  for (int i=0; i<n; i++)
875
878
    {
1024
1027
}
1025
1028
 
1026
1029
 
1027
 
void svdcmp(vector<vector<double> > & a, 
 
1030
bool svdcmp(vector<vector<double> > & a, 
1028
1031
            vector<double> & w, 
1029
1032
            vector<vector<double> > &v)
1030
1033
{
1160
1163
        break;
1161
1164
      }
1162
1165
      if (its == 29) 
1163
 
        error("SVD function cannot converge: multicollinearity issues?");
 
1166
        return false; // cannot converge: multi-collinearity?
1164
1167
      x=w[l];
1165
1168
      nm=k-1;
1166
1169
      y=w[nm];
1211
1214
      w[k]=x;
1212
1215
    }
1213
1216
  }
 
1217
  return true;
1214
1218
}
1215
1219
 
1216
1220
 
1528
1532
    d[i][i] = v[i];
1529
1533
  return d;
1530
1534
}
 
1535
 
 
1536
double rnorm()
 
1537
{
 
1538
 
 
1539
  double u1 = CRandom::rand();
 
1540
  double u2 = CRandom::rand();
 
1541
  return sqrt(-2*log(u1)) * cos(2*M_PI*u2);
 
1542
 
 
1543
  // z2 = sqrt(-2*log(u1)) * sin(2*M_PI*u2);
 
1544
}