1 /* 2 This file is part of the BOLT-LMM linear mixed model software package 3 developed by Po-Ru Loh. Copyright (C) 2014-2019 Harvard University. 4 5 This program is free software: you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation, either version 3 of the License, or 8 (at your option) any later version. 9 10 This program is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU General Public License for more details. 14 15 You should have received a copy of the GNU General Public License 16 along with this program. If not, see <http://www.gnu.org/licenses/>. 17 */ 18 19 #include <cmath> 20 #include <cstdlib> 21 #include <utility> 22 #include <algorithm> 23 24 #include "NumericUtils.hpp" 25 26 namespace NumericUtils { sum(const double x[],uint64 N)27 double sum(const double x[], uint64 N) { 28 double ans = 0; 29 for (uint64 n = 0; n < N; n++) 30 ans += x[n]; 31 return ans; 32 } mean(const std::vector<double> & x)33 double mean(const std::vector <double> &x) { 34 uint64 N = x.size(); return sum(&x[0], N) / N; 35 } median(std::vector<double> x)36 double median(std::vector <double> x) { 37 std::sort(x.begin(), x.end()); 38 int n = x.size(); 39 return (x[n/2] + x[(n-1)/2]) / 2; 40 } 41 // takes into account that some 0 values may indicate missing/ignored: divide out by Nused, not N mean(const double x[],uint64 N,uint64 Nused)42 double mean(const double x[], uint64 N, uint64 Nused) { 43 return sum(x, N) / Nused; 44 } 45 // regress y on x, assuming both have been 0-centered (so 0-filled missing values ok) regCoeff(const double y[],const double x[],uint64 N)46 double regCoeff(const double y[], const double x[], uint64 N) { 47 /* WRONG! if not mean-centered already, need to mask missing indivs in loop 48 double xbar = mean(x, N, Nused); 49 double ybar = mean(y, N, Nused); 50 cout << "xbar: " << xbar << " ybar: " << ybar << endl; 51 double numer = 0, denom = 0; 52 for (uint64 n = 0; n < N; n++) { 53 numer += (x[n]-xbar) * (y[n]-ybar); 54 denom += sq(x[n]-xbar); 55 } 56 */ 57 double numer = 0, denom = 0; 58 for (uint64 n = 0; n < N; n++) { 59 numer += x[n] * y[n]; 60 denom += sq(x[n]); 61 } 62 return numer / denom; 63 } dot(const double x[],const double y[],uint64 N)64 double dot(const double x[], const double y[], uint64 N) { 65 double ans = 0; 66 for (uint64 n = 0; n < N; n++) 67 ans += x[n] * y[n]; 68 return ans; 69 } norm2(const double x[],uint64 N)70 double norm2(const double x[], uint64 N) { 71 double ans = 0; 72 for (uint64 n = 0; n < N; n++) 73 ans += sq(x[n]); 74 return ans; 75 } normalize(double x[],uint64 N)76 void normalize(double x[], uint64 N) { 77 double scale = 1.0 / sqrt(norm2(x, N)); 78 for (uint64 n = 0; n < N; n++) 79 x[n] *= scale; 80 } 81 meanStdDev(const double x[],uint64 N)82 std::pair <double, double> meanStdDev(const double x[], uint64 N) { 83 double mu = 0, s2 = 0; 84 for (uint64 n = 0; n < N; n++) mu += x[n]; 85 mu /= N; 86 for (uint64 n = 0; n < N; n++) s2 += sq(x[n]-mu); 87 s2 /= (N-1); 88 return std::make_pair(mu, sqrt(s2)); 89 } meanStdErr(const double x[],uint64 N)90 std::pair <double, double> meanStdErr(const double x[], uint64 N) { 91 std::pair <double, double> ret = meanStdDev(x, N); 92 ret.second /= sqrt((double) N); 93 return ret; 94 } meanStdDev(const std::vector<double> & x)95 std::pair <double, double> meanStdDev(const std::vector <double> &x) { 96 return meanStdDev(&x[0], x.size()); 97 } meanStdErr(const std::vector<double> & x)98 std::pair <double, double> meanStdErr(const std::vector <double> &x) { 99 return meanStdErr(&x[0], x.size()); 100 } 101 } 102