1 /* 2 * Interface to the statistical functions 3 * 4 * 2010 by Jian Yang <jian.yang@qimr.edu.au> 5 * 6 * This file is distributed under the GNU General Public 7 * License, Version 2. Please see the file COPYING for more 8 * details 9 */ 10 11 #ifndef _STATFUNC_H 12 #define _STATFUNC_H 13 14 #include <vector> 15 #include <algorithm> 16 #include <iostream> 17 #include "CommFunc.h" 18 #include "dcdflib.h" 19 #include <Eigen/Dense> 20 #include <Eigen/Sparse> 21 #include <unsupported/Eigen/SparseExtra> 22 23 using namespace Eigen; 24 using namespace std; 25 26 namespace StatFunc 27 { 28 ////////// P-value Calculatiion Functions Start //////////////// 29 // if two_tail=true, return two tail test probability; 30 // otherwise, return one tail test probability; 31 double t_prob(double df, double t_value, bool two_tail = true); 32 33 double F_prob(double df_1, double df_2, double F_value); 34 35 double betai(double a, double b, double x); 36 37 double gammln(double xx); 38 39 double betacf(double a, double b, double x); 40 41 double chi_prob(double df, double chi_sqr_val); 42 43 double gammp(const double a, const double x); 44 45 void gser(double &gamser, const double a, const double x, double &gln); 46 47 void gcf(double &gammcf, const double a, const double x, double &gln); 48 ////////// P-value Calculatiion Functions End //////////////// 49 50 ///////// Random Number Generation Functions Start //////// 51 // (0, 1) uniform distribution generator 52 double ran1(int &idum); 53 54 // (a, b) uniform distribution generator 55 double UniformDev(double a, double b, int &idum); 56 57 // normal distribution generator 58 double gasdev(int &idum); 59 60 // generate a sequence following normal distribution 61 void gasdev_seq(int &idum, vector<double> &vec, int size, double means, double var); 62 63 // generate a sequence following normal distribution with zero average 64 void gasdev_seq(int &idum, vector<double> &vec, int size, double var); 65 66 // generate a gamma random variable when alpha is larger than 1.0 67 double cheng_gamdev(int &idum, const double alpha); 68 69 // generate a random variable follow chi-square distribution 70 double chidev(int &idum, const double df); 71 72 // generate a integer between a and b 73 int RandAbs(int a, int b, int &seed); 74 75 // Function for get a Chi value of right tail when given a df and prob. 76 double chi_val(double df, double prob); 77 78 // Function for get a t value of right tail when given a df and prob. 79 double t_val(double df, double prob); 80 81 // Function for get a F value of right tail when given a df and prob. 82 double F_val(double df_1, double df_2, double prob); 83 84 // Control the experimental-wise type I error by FDR method 85 double ControlFDR(const vector<double> &P_Value, double alpha, bool Restrict); 86 double ControlFDR_Zou(const vector<double> &GenePValue, double FDR); 87 double ControlFDR_Storey(vector<double> &P_Value, vector<double> &Q_Value, double CrtQ, double &FDR); 88 double CalcuPi0(vector<double> &P_Value, vector<double> &Lambda); 89 void spline(vector<double> &x, vector<double> &y, const double yp1, const double ypn, vector<double> &y2); 90 void splint(vector<double> &xa, vector<double> &ya, vector<double> &y2a, const double x, double &y); 91 92 // normal distribution 93 double erf(double x); 94 double pnorm(double x); 95 double dnorm(double x); 96 double qnorm_sub(double x, double y); 97 double qnorm(double p, bool upper = true); 98 99 // chisq distribution 100 double pchisq(double x, double df); 101 double qchisq(double q, double df); 102 103 // sum of chisq distribution 104 double pchisqsum(double x, VectorXd lambda); 105 double psadd(double x, VectorXd lambda); 106 double psatt(double x, VectorXd lambda); 107 double K(double zeta, VectorXd &lambda); 108 double Kp(double zeta, VectorXd &lambda); 109 double Kpp(double zeta, VectorXd &lambda); 110 double Kp_min_x(double zeta, VectorXd &lambda, double x); 111 double Brents_Kp_min_x(VectorXd &lambda, double x, double lowerLimit, double upperLimit, double errorTol); 112 } 113 114 #endif 115