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 <vector> 20 #include <utility> 21 #include <cmath> 22 #include "Jackknife.hpp" 23 24 namespace Jackknife { 25 26 using std::vector; 27 using std::pair; 28 using std::make_pair; 29 stddev(const vector<double> & x,int n)30 double stddev(const vector <double> &x, int n) { 31 for (int i = 0; i < n; i++) if (std::isnan(x[i])) return NAN; 32 for (int i = 0; i < n; i++) if (std::isinf(x[i])) return INFINITY; 33 double s = 0.0, s2 = 0.0; 34 for (int i = 0; i < n; i++) { 35 s += x[i]; 36 s2 += x[i]*x[i]; 37 } 38 return sqrt((s2 - s*s/n) * (n-1) / n); 39 } 40 41 // assumes last element of x is the leave-no-data-out estimator (previous are jackknife reps) mean_std(const vector<double> & x)42 pair <double, double> mean_std(const vector <double> &x) { 43 int n = x.size()-1; // number of jackknife reps 44 return make_pair(x[n], stddev(x, n)); 45 } 46 47 // assumes last element of x is the leave-no-data-out estimator (previous are jackknife reps) zscore(const vector<double> & x)48 double zscore(const vector <double> &x) { 49 pair <double, double> mu_sigma = mean_std(x); 50 return mu_sigma.first / mu_sigma.second; 51 } 52 53 // assumes last element of x is the leave-no-data-out estimator (previous are jackknife reps) diff_mean_std(vector<double> x,const vector<double> & x_ref)54 pair <double, double> diff_mean_std(vector <double> x, const vector <double> &x_ref) { 55 int n = x.size()-1; // number of jackknife reps 56 for (int i = 0; i <= n; i++) x[i] -= x_ref[i]; 57 return make_pair(x[n], stddev(x, n)); 58 } 59 ratioOfSumsMeanStd(const vector<double> & x,const vector<double> & y)60 pair <double, double> ratioOfSumsMeanStd(const vector <double> &x, const vector <double> &y) { 61 int n = x.size(); 62 vector <double> ratioOfSumsJacks(n+1); 63 for (int j = 0; j <= n; j++) { 64 double xSum = 0, ySum = 0; 65 for (int i = 0; i < n; i++) 66 if (i != j) { 67 xSum += x[i]; 68 ySum += y[i]; 69 } 70 ratioOfSumsJacks[j] = xSum / ySum; 71 } 72 return mean_std(ratioOfSumsJacks); 73 } 74 75 } 76