| 96 | } |
| 97 | |
| 98 | void LinearRegression::computeGoodness_(const std::vector<double>& X, const std::vector<double>& Y, double confidence_interval_P) |
| 99 | { |
| 100 | OPENMS_PRECONDITION(static_cast<unsigned>(X.size() == Y.size()), |
| 101 | "Fitted X and Y have different lengths."); |
| 102 | OPENMS_PRECONDITION(static_cast<unsigned>(X.size()) > 2, |
| 103 | "Cannot compute goodness of fit for regression with less than 3 data points"); |
| 104 | // specifically, boost throws an exception for a t-distribution with zero df |
| 105 | |
| 106 | Size N = X.size(); |
| 107 | |
| 108 | // Mean of abscissa and ordinate values |
| 109 | double x_mean = Math::mean(X.begin(), X.end()); |
| 110 | double y_mean = Math::mean(Y.begin(), Y.end()); |
| 111 | |
| 112 | // Variance and Covariances |
| 113 | double var_X = Math::variance(X.begin(), X.end(), x_mean); |
| 114 | double var_Y = Math::variance(Y.begin(), Y.end(), y_mean); |
| 115 | double cov_XY = Math::covariance(X.begin(), X.end(), Y.begin(), Y.end()); |
| 116 | |
| 117 | // S_xx |
| 118 | double s_XX = var_X * (N-1); |
| 119 | /*for (unsigned i = 0; i < N; ++i) |
| 120 | { |
| 121 | double d = (X[i] - x_mean); |
| 122 | s_XX += d * d; |
| 123 | }*/ |
| 124 | |
| 125 | // Compute the squared Pearson coefficient |
| 126 | r_squared_ = (cov_XY * cov_XY) / (var_X * var_Y); |
| 127 | |
| 128 | // The standard deviation of the residuals |
| 129 | double sum = 0; |
| 130 | for (unsigned i = 0; i < N; ++i) |
| 131 | { |
| 132 | double x_i = fabs(Y[i] - (intercept_ + slope_ * X[i])); |
| 133 | sum += x_i; |
| 134 | } |
| 135 | mean_residuals_ = sum / N; |
| 136 | stand_dev_residuals_ = sqrt((chi_squared_ - (sum * sum) / N) / (N - 1)); |
| 137 | |
| 138 | // The Standard error of the slope |
| 139 | stand_error_slope_ = stand_dev_residuals_ / sqrt(s_XX); |
| 140 | |
| 141 | // and the intersection of Y_hat with the x-axis |
| 142 | x_intercept_ = -(intercept_ / slope_); |
| 143 | |
| 144 | double P = 1 - (1 - confidence_interval_P) / 2; |
| 145 | boost::math::students_t tdist(N - 2); |
| 146 | t_star_ = boost::math::quantile(tdist, P); |
| 147 | |
| 148 | //Compute the asymmetric 95% confidence interval of around the X-intercept |
| 149 | double g = (t_star_ / (slope_ / stand_error_slope_)); |
| 150 | g *= g; |
| 151 | double left = (x_intercept_ - x_mean) * g; |
| 152 | double bottom = 1 - g; |
| 153 | double d = (x_intercept_ - x_mean); |
| 154 | double right = t_star_ * (stand_dev_residuals_ / slope_) * sqrt((d * d) / s_XX + (bottom / N)); |
| 155 | |