... ... @@ -35,32 +35,36 @@ my $header = qq!#pragma once namespace fracdist { constexpr size_t /// The number of q values (from 1 to$DataParser::NUM_Q) in fracdist::q_const and fracdist::q_noconst. q_length = $DataParser::NUM_Q, /// The number of b values in fracdist::bvalues. b_length == bvalues.size() will be true. b_length =$bvalues, /** The number of probability values and associated quantiles in fracdist::pvalues, * fracdist::q_const[i][j], and fracdist::q_noconst[i][j] (for any admissable i and j). */ p_length = $pvalues; /** The bvalues: bvalues[j]' is the b value corresponding to the quantiles contained in * q_const[i][j]' /** The bvalues: bvalues[j] is the b value corresponding to the quantiles contained in * fracdist::q_const[i][j] */ extern const std::array bvalues; /** The pvalues; pvalues[k]' is the pvalue associated with quantile * q_const[i][j][k]' /** The pvalues; pvalues[k] is the pvalue associated with quantile * fracdist::q_const[i][j][k] */ extern const std::array pvalues; /** A double[][][] (wrapped in std::array's) where double[x][y][z] corresponds to the quantile with q=z+1', b=bvalues[y]', and pval=pvalues[z]'. For example, if bvalues == 0.75' and pvalues == 0.05' then quantiles' is the 0.05 quantile * for q=4, b=0.75. This variable is for models estimated with a constant. /** A double[][][] (wrapped in nested std::array) where q_const[x][y][z] corresponds to the quantile with \\f\$q = z+1\\f\$, \\f\$b={}\\f\$fracdist::bvalues[y], and \\f\$p={}\\f\$fracdist::pvalues[z]. For example, if bvalues == 0.75 and pvalues == 0.05 then q_const is the 0.05 quantile for \\f\$q=4, b=0.75\\f\$for a model with a constant. This variable is for models estimated *with* a constant. */ extern const std::array, b_length>, q_length> q_const; /** A double[][][] (wrapped in std::array's) where double[x][y][z] corresponds to the quantile with q=z+1', b=bvalues[y]', and pval=pvalues[z]'. For example, if bvalues == 0.75' and pvalues == 0.05' then quantiles' is the 0.05 quantile for q=4, b=0.75. This variable is for models estimated without a constant. /** A double[][][] (wrapped in nested std::array) where q_noconst[x][y][z] corresponds to the quantile with \\f\$q = z+1\\f\$, \\f\$b={}\\f\$fracdist::bvalues[y], and \\f\$p={}\\f\$fracdist::pvalues[z]. For example, if bvalues == 0.75 and pvalues == 0.05 then q_noconst is the 0.05 quantile for \\f\$q=4, b=0.75\\f\$for a model without a constant. This variable is for models estimated *without* a constant. */ extern const std::array, b_length>, q_length> q_noconst; ... ...  /** @file fdcrit.c /** @file fdcrit.cpp * @brief Takes q value, b value, constant (1 or 0), and test levels, outputs critical values. * * If any invalid arguments are provided, a help message is written to stderr and the program exits ... ...  /** @file fdpval.c /** @file fdpval.cpp * @brief Takes q value, b value, constant (1 or 0), and test stats, outputs p-values. * * If any invalid arguments are provided, a help message is written to stderr and the program exits ... ...  ... ... @@ -175,7 +175,7 @@ std::pair find_bracket(const size_t ¢er, const size_t &max, unsigned int chisq_inv_cache_q = 0; std::array chisq_inv_cache; boost::math::chi_squared_distribution chisq_dist(1); /* Returns the inverse chi squared cdf at pvalues[pval_index]' with q^2 degrees of freedom. The /** Returns the inverse chi squared cdf at pvalues[pval_index] with q^2 degrees of freedom. The * value is cached (so long as q doesn't change) so that subsequent calls for the same value are * very fast. */ ... ...  ... ... @@ -4,44 +4,56 @@ #include #include /** @file fracdist/common.hpp * @brief Header file for various common fracdist functionality. */ /** Namespace for all fracdist library code. */ namespace fracdist { /** An enum giving the different quantile interpolation types supported by pvalue_advanced and * critical_advanced. /** An enum giving the different quantile interpolation types supported by pvalue_advanced() and * critical_advanced(). * * \todo add a spline interpolation mode */ enum class interpolation { /** Quadratic fitting of nearby points as described in MacKinnon and Nielsen (2014). This * always uses quadratic approximation across nearby b values, even when the requested b value * is one of the ones in the data file. This interpolation method gives smoother curves across * b' values than the other two methods, but is slightly less accurate at known b' * values (0.51, 0.55, 0.6, 0.65, ..., 1.95, 2.0). * b values than the other two methods, but is slightly less accurate at known b * values \f$(0.51, 0.55, 0.6, 0.65, \ldots, 1.95, 2.0)\f$. */ JGMMON14, /** Like JGMMON14, but when a b value is requested that exactly matches a b value in the /** Like interpolation::JGMMON14, but when a b value is requested that exactly matches a b value in the * quantile data, the exact data quantiles are used. Otherwise, interpolation occurs as in * JGMMON14. This has the advantage of offering more precise values for known b' values, but * the disadvantage that there are discontinuities in the calculated quantiles at the known b' * interpolation::JGMMON14. This has the advantage of offering more precise values for known b values, but * the disadvantage that there are discontinuities in the calculated quantiles at the known b * values. */ exact_or_JGMMON14, /** Linear interpolation between bracketing quantiles. If, for example, b=0.69 is provided * but the data only has b=0.65' and b=0.7', the resulting quantiles will be the weighted sum * of \f$0.2q_{b=0.65} + 0.8q_{b=0.7}\f$of the two quantiles. Like exact_or_JGMMON14, this * returns exactly the data's quantiles for an exact match of b' value. Unlike * exact_or_JGMMON14, this method has no discontinuities for changes in b' (but doesn't have * kinks at each known b' value). * but the data only has quantiles \f$q_{0.65}\f$and \f$q_{0.7}\f$(for \f$b=0.65\f$and * \f$b=0.7\f$), the resulting quantiles will be the weighted sum of \f$0.2q_{0.65} + * 0.8q_{0.7}\f$of the two quantiles. Like interpolation::exact_or_JGMMON14, this returns exactly the data's * quantiles for an exact match of b value. Unlike interpolation::exact_or_JGMMON14, this method has no * discontinuities for changes in b (but does have kinks at each known b value). */ linear }; /** Takes q, b, constant, and interpolation mode values and calculates the quantiles for the given * values. If anything is invalid, throws an exception. /** Takes \f$q\f$, \f$b\f$, constant, and interpolation mode values and calculates the quantiles for * the given set of values. If any of the values is invalid, throws an exception. * * The result of the previous call is cached so that calling quantiles() a second time with the same * q, b, constant, and interp values will not re-perform the necessary calculations. * * This function is mainly used for internal use by the other functions in this file, but may be * useful for other purposes. * * \throws std::out_of_range for an invalid b, q value * \throws std::runtime_error if approx_points is there are no enough data points to estimate a * quadratic approximation. This shouldn't happen, normally, as the data and weights ensure that * there will always be enough data points. */ const std::array quantiles(const unsigned int &q, const double &b, const bool &constant, const interpolation &interp); ... ... @@ -63,32 +75,36 @@ size_t find_closest(const double &value, const Container &array) { return min_at; } /** Finds a bracket of size at most size' of indices centered (if possible) on the given index. If the * given index is too close to 0 or max', the first and last values are truncated to the end points * (and a bracket smaller than n' results). /** Finds a bracket of size at most size of indices centered (if possible) on the given index. If the * given index is too close to 0 or max, the first and last values are truncated to the end points * (and a bracket smaller than n results). * * Returns an std::pair with .first being the first bracket element, .second being * the last bracket element. last-first+1 <= size is guaranteed. * Returns an std::pair with pair.first being the first bracket element, * pair.second being the last bracket element. \f$pair_{second}-pair_{first}+1 \leq size\f$is * guaranteed; the weak inequality results from end-point truncation. */ std::pair find_bracket(const size_t ¢er, const size_t &max, const size_t &size); /** Returns the inverse chi squared cdf at pvalues[pval_index]' with q^2 degrees of freedom. The * value is cached (so long as q doesn't change) so that subsequent calls for the same value are * very fast. /** Returns the inverse chi squared cdf at pvalues[pval_index] with \f$q^2\f$degrees of freedom. * The value is cached (so long as the same q is used) so that subsequent calls for the same value * are very fast. */ double chisq_inv_p_i(const size_t &pval_index, const unsigned int &q); /** Wrapper object around std::ostringstream that overrides << and is castable to a std::string so * that a construction like: /** Wrapper object around std::ostringstream that overrides the << operator and is castable to a * std::string. Thie is primarily intended for quickly building strings via a construction such as: * * somethingrequiringastring(fdstringstream() << "a" << "b") * somethingrequiringastring(fdstringstream() << "value: " << var); * * is valid without requiring the caller to store an intermediate ostringstream object. */ class ostringstream : public std::ostringstream { public: /// Constructs a new empty std::ostringstream wrapper ostringstream() : std::ostringstream() {} /// Forwards anything shifted onto this object to std::ostringstream template ostringstream& operator<<(const T &v) { std::ostringstream::operator<<(v); return *this; } /// Can be cast implicitly to a std::string whenever required. operator std::string() const { return str(); } }; ... ...  #pragma once #include /** @file fracdist/critical.hpp * @brief Header file for fracdist's interface to finding a critical test statistic from a test * level. */ namespace fracdist { /** Calculates a critical value for a given level of the test. Takes the level, q value, b value, * and whether the model contains a constant (0 for no constant, anything else for constant). /** Calculates a critical value for a given level of the test. Takes the test level, \f$q\f$value, * \f$b\f$value, and whether the model contains a constant (true for a constant, false for no * constant). * * \throws std::out_of_range for an invalid b or q value, or for a test level outside [0, 1]. * \throws std::out_of_range for an invalid b or q value, or for a test level outside \f$[0, 1]\f\$. */ double critical(const double &test_level, const unsigned int &q, const double &b, const bool &constant); /** Like critical(), but also takes an interpolation mode and number of P-value approximation * points. approx_points' must be at least 3 (and depending on the test_stat and parameters, might * points. approx_points must be at least 3 (and depending on the test_stat and parameters, might * need to be at least 5). * * Note that for values near the limit of the data (i.e. with pvalues close to 0 or 1), fewer points ... ... @@ -21,9 +26,9 @@ double critical(const double &test_level, const unsigned int &q, const double &b * * \throws std::out_of_range for an invalid b or q value, or for a test level outside [0, 1]. * \throws std::runtime_error if approx_points is too small to perform the required quadratic * approximation (in other words, fewer than 3 points). This will happen with approx_points < 5' * approximation (in other words, fewer than 3 points). This will happen with approx_points < 5 * for test_stats closest to those associated with limit p-values (0.0001 and 0.9999). Thus, while * approx_points' of 3 or 4 may work for some test_stat' values, 5 is the minimum value that never * approx_points of 3 or 4 may work for some test_stat values, 5 is the minimum value that never * results in this error. */ double critical_advanced(double test_level, const unsigned int &q, const double &b, const bool &constant, ... ...
 #pragma once #include /** @file fracdist/pvalue.hpp * @brief Header file for fracdist's interface to finding a pvalue from a test statistic. */ namespace fracdist { /** Calculates a p-value for a given test statistic, q value, b value, and whether the model ... ... @@ -17,7 +21,7 @@ namespace fracdist { double pvalue(const double &test_stat, const unsigned int &q, const double &b, const bool &constant); /** Like pvalue(), but requires an interpolation mode and number of P-value approximation * points. approx_points' must be at least 3 (and depending on the test_stat and parameters, might * points. approx_points must be at least 3 (and depending on the test_stat and parameters, might * need to be at least 5). * * Note that for values near the limit of the data (i.e. with pvalues close to 0 or 1), fewer points ... ... @@ -27,9 +31,9 @@ double pvalue(const double &test_stat, const unsigned int &q, const double &b, c * * \throws std::out_of_range for an invalid b, q value * \throws std::runtime_error if approx_points is too small to perform the required quadratic * approximation (in other words, fewer than 3 points). This will happen with approx_points < 5' * approximation (in other words, fewer than 3 points). This will happen with approx_points < 5 * for test_stats closest to those associated with limit p-values (0.0001 and 0.9999). Thus, while * approx_points' of 3 or 4 may work for some test_stat' values, 5 is the minimum value that never * approx_points of 3 or 4 may work for some test_stat values, 5 is the minimum value that never * results in this error. */ double pvalue_advanced(const double &test_stat, const unsigned int &q, const double &b, const bool &constant, ... ...