*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
* to the Inverse Beta Distribution to Find Percentiles of the F-Distribution,"
* ACM Transactions on Mathematical Software, volume 19, number 4, December 1993,
* pages 474-480.
*
* to the Inverse Beta Distribution to Find Percentiles of the F-Distribution,"
* ACM Transactions on Mathematical Software, volume 19, number 4, December 1993,
* pages 474-480.
*
* Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8,
* August 1968, pages 1264-1273.
*/
* Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8,
* August 1968, pages 1264-1273.
*/
s_bisect (double x, double y)
{
double result = GSL_MIN(x,y) + fabs(x - y) / 2.0;
return result;
}
static double
s_bisect (double x, double y)
{
double result = GSL_MIN(x,y) + fabs(x - y) / 2.0;
return result;
}
static double
-new_guess_P ( double old_guess, double x, double y,
+new_guess_P ( double old_guess, double x, double y,
-new_guess_Q ( double old_guess, double x, double y,
+new_guess_Q ( double old_guess, double x, double y,
* three terms of the Cornish-Fisher expansion
* without recursion. The recursive functions
* make the code more legible when higher order coefficients
* three terms of the Cornish-Fisher expansion
* without recursion. The recursive functions
* make the code more legible when higher order coefficients
gam_ab = gsl_sf_lngamma(a + b);
gam_a = gsl_sf_lngamma (a);
gam_b = gsl_sf_lngamma (b);
gam_ab = gsl_sf_lngamma(a + b);
gam_a = gsl_sf_lngamma (a);
gam_b = gsl_sf_lngamma (b);
- * Theoretically, use of this term should give a better approximation,
- * but in practice inclusion of the cubic term worsens the
+ * Theoretically, use of this term should give a better approximation,
+ * but in practice inclusion of the cubic term worsens the
* iterative procedure in gsl_cdf_beta_Pinv and gsl_cdf_beta_Qinv
* for extreme values of p, a or b.
* iterative procedure in gsl_cdf_beta_Pinv and gsl_cdf_beta_Qinv
* for extreme values of p, a or b.
* starting with the nth derivative of s_psi = -f'(x)/f(x),
* where f is the beta density.
*
* starting with the nth derivative of s_psi = -f'(x)/f(x),
* where f is the beta density.
*
asgn = (n % 2) ? 1.0:-1.0;
bsgn = (n % 2) ? -1.0:1.0;
result = gsl_sf_gamma(np1) * ((bsgn * bm1 / (pow(mx, np1)))
asgn = (n % 2) ? 1.0:-1.0;
bsgn = (n % 2) ? -1.0:1.0;
result = gsl_sf_gamma(np1) * ((bsgn * bm1 / (pow(mx, np1)))
if (n == 2)
{
result = s_d_psi(x, a, b, k);
if (n == 2)
{
result = s_d_psi(x, a, b, k);
k_fac = gsl_sf_lngamma(k+1.0);
i_fac = gsl_sf_lngamma(i+1.0);
kmi_fac = gsl_sf_lngamma(k-i+1.0);
k_fac = gsl_sf_lngamma(k+1.0);
i_fac = gsl_sf_lngamma(i+1.0);
kmi_fac = gsl_sf_lngamma(k-i+1.0);
- * get_d_coeff( x, a, b, 2.0, i)
+ * get_d_coeff( x, a, b, 2.0, i)
* get_d_coeff( x, a, b, (n - 1.0), (k - i));
}
result += get_d_coeff ( x, a, b, (n-1.0), (k+1.0));
* get_d_coeff( x, a, b, (n - 1.0), (k - i));
}
result += get_d_coeff ( x, a, b, (n-1.0), (k+1.0));
* 0, an initial state value of a/(a+b) will
* cause the interpolating polynomial
* to overflow.
* 0, an initial state value of a/(a+b) will
* cause the interpolating polynomial
* to overflow.
- tmp = new_guess_P ( state, lower, upper,
+ tmp = new_guess_P ( state, lower, upper,
/*
* The cubic term does not help, and can can
* harm the approximation for extreme values of
/*
* The cubic term does not help, and can can
* harm the approximation for extreme values of
#if 0
c3 = get_corn_fish_cube (state, a, b);
state += err * (c1 + (err / 2.0 ) * (c2 + c3 * err / 3.0));
#if 0
c3 = get_corn_fish_cube (state, a, b);
state += err * (c1 + (err / 2.0 ) * (c2 + c3 * err / 3.0));
* The recursion makes coding higher-order terms
* easier, but did not improve the result beyond
* the use of three terms. Since explicitly coding
* those three terms in the get_corn_fish_* functions
* was not difficult, the recursion was abandoned.
*/
* The recursion makes coding higher-order terms
* easier, but did not improve the result beyond
* the use of three terms. Since explicitly coding
* those three terms in the get_corn_fish_* functions
* was not difficult, the recursion was abandoned.
*/
coeff = 1.0;
for(i = 1.0; i < BETADISTINV_N_TERMS; i += 1.0)
{
i_fac *= i;
coeff = get_corn_fish (coeff, prior_state, a, b, i);
coeff = 1.0;
for(i = 1.0; i < BETADISTINV_N_TERMS; i += 1.0)
{
i_fac *= i;
coeff = get_corn_fish (coeff, prior_state, a, b, i);
(i_fac * pow (gsl_ran_beta_pdf(prior_state,a,b), i));
}
#endif
(i_fac * pow (gsl_ran_beta_pdf(prior_state,a,b), i));
}
#endif
* When q is close to 0, the bisection
* and interpolation done in the rest of
* this routine will not give the correct
* When q is close to 0, the bisection
* and interpolation done in the rest of
* this routine will not give the correct
* gsl_cdf_beta_Qinv is called instead.
*/
state = gslextras_cdf_beta_Pinv ( q, a, b);
* gsl_cdf_beta_Qinv is called instead.
*/
state = gslextras_cdf_beta_Pinv ( q, a, b);
- tmp = new_guess_Q ( state, lower, upper,
+ tmp = new_guess_Q ( state, lower, upper,