From: Ben Pfaff Date: Sun, 21 Nov 2021 21:16:19 +0000 (-0800) Subject: All the distribution functions, though they're untested. X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8539b9672ca634e0bedf7a531709e845a6b451d6;p=pspp All the distribution functions, though they're untested. --- diff --git a/src/language/expressions/operations.def b/src/language/expressions/operations.def index df2887044f..42728e46ea 100644 --- a/src/language/expressions/operations.def +++ b/src/language/expressions/operations.def @@ -872,14 +872,14 @@ function NPDF.T (x, df > 0, nc) = unimplemented; // Type-1 Gumbel distribution. extension function CDF.T1G (x, a, b) = gsl_cdf_gumbel1_P (x, a, b); extension function IDF.T1G (P >= 0 && P <= 1, a, b) - = gsl_cdf_gumbel1_P (P, a, b); + = gsl_cdf_gumbel1_Pinv (P, a, b); extension function PDF.T1G (x, a, b) = gsl_ran_gumbel1_pdf (x, a, b); no_opt extension function RV.T1G (a, b) = gsl_ran_gumbel1 (get_rng (), a, b); // Type-2 Gumbel distribution. extension function CDF.T2G (x, a, b) = gsl_cdf_gumbel2_P (x, a, b); extension function IDF.T2G (P >= 0 && P <= 1, a, b) - = gsl_cdf_gumbel2_P (P, a, b); + = gsl_cdf_gumbel2_Pinv (P, a, b); extension function PDF.T2G (x, a, b) = gsl_ran_gumbel2_pdf (x, a, b); no_opt extension function RV.T2G (a, b) = gsl_ran_gumbel2 (get_rng (), a, b); diff --git a/src/language/stats/matrix.c b/src/language/stats/matrix.c index b709d65e01..b507f022fa 100644 --- a/src/language/stats/matrix.c +++ b/src/language/stats/matrix.c @@ -181,31 +181,31 @@ matrix_var_set (struct matrix_var *var, gsl_matrix *value) } #define MATRIX_FUNCTIONS \ - F(ABS, "ABS", m_m_e, NULL) \ + F(ABS, "ABS", m_e, NULL) \ F(ALL, "ALL", d_m, NULL) \ F(ANY, "ANY", d_m, NULL) \ - F(ARSIN, "ARSIN", m_m_e, "a[-1,1]") \ - F(ARTAN, "ARTAN", m_m_e, NULL) \ + F(ARSIN, "ARSIN", m_e, "a[-1,1]") \ + F(ARTAN, "ARTAN", m_e, NULL) \ F(BLOCK, "BLOCK", m_any, NULL) \ F(CHOL, "CHOL", m_m, NULL) \ F(CMIN, "CMIN", m_m, NULL) \ F(CMAX, "CMAX", m_m, NULL) \ - F(COS, "COS", m_m_e, NULL) \ + F(COS, "COS", m_e, NULL) \ F(CSSQ, "CSSQ", m_m, NULL) \ F(CSUM, "CSUM", m_m, NULL) \ F(DESIGN, "DESIGN", m_m, NULL) \ F(DET, "DET", d_m, NULL) \ F(DIAG, "DIAG", m_m, NULL) \ F(EVAL, "EVAL", m_m, NULL) \ - F(EXP, "EXP", m_m_e, NULL) \ + F(EXP, "EXP", m_e, NULL) \ F(GINV, "GINV", m_m, NULL) \ F(GRADE, "GRADE", m_m, NULL) \ F(GSCH, "GSCH", m_m, NULL) \ F(IDENT, "IDENT", IDENT, NULL) \ F(INV, "INV", m_m, NULL) \ F(KRONEKER, "KRONEKER", m_mm, NULL) \ - F(LG10, "LG10", m_m_e, "a>0") \ - F(LN, "LN", m_m_e, "a>0") \ + F(LG10, "LG10", m_e, "a>0") \ + F(LN, "LN", m_e, "a>0") \ F(MAGIC, "MAGIC", m_d, NULL) \ F(MAKE, "MAKE", m_ddd, NULL) \ F(MDIAG, "MDIAG", m_v, NULL) \ @@ -220,11 +220,11 @@ matrix_var_set (struct matrix_var *var, gsl_matrix *value) F(RESHAPE, "RESHAPE", m_mdd, NULL) \ F(RMAX, "RMAX", m_m, NULL) \ F(RMIN, "RMIN", m_m, NULL) \ - F(RND, "RND", m_m_e, NULL) \ + F(RND, "RND", m_e, NULL) \ F(RNKORDER, "RNKORDER", m_m, NULL) \ F(RSSQ, "RSSQ", m_m, NULL) \ F(RSUM, "RSUM", m_m, NULL) \ - F(SIN, "SIN", m_m_e, NULL) \ + F(SIN, "SIN", m_e, NULL) \ F(SOLVE, "SOLVE", m_mm, NULL) \ F(SQRT, "SQRT", m_m, NULL) \ F(SSCP, "SSCP", m_m, NULL) \ @@ -233,36 +233,118 @@ matrix_var_set (struct matrix_var *var, gsl_matrix *value) F(T, "T", m_m, NULL) \ F(TRACE, "TRACE", d_m, NULL) \ F(TRANSPOS, "TRANSPOS", m_m, NULL) \ - F(TRUNC, "TRUNC", m_m_e, NULL) \ + F(TRUNC, "TRUNC", m_e, NULL) \ F(UNIFORM, "UNIFORM", m_dd, NULL) \ - F(PDF_BETA, "PDF.BETA", m_mdd_e, "a[0,1] b>0 c>0") \ - F(CDF_BETA, "CDF.BETA", m_mdd_e, "a[0,1] b>0 c>0") \ - F(IDF_BETA, "IDF.BETA", m_mdd_e, "a[0,1] b>0 c>0") \ + F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \ + F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \ + F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0") \ F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \ - F(NCDF_BETA, "NCDF.BETA", m_mddd_e, "a>=0 b>0 c>0 d>0") \ - F(NPDF_BETA, "NCDF.BETA", m_mddd_e, "a>=0 b>0 c>0 d>0") \ - /* XXX CDF.BVNOR */ \ - F(PDF_BVNOR, "PDF.BVNOR", m_mdd_e, "c[-1,1]") \ - F(CDF_CAUCHY, "CDF.CAUCHY", m_mdd_e, "c>0") \ - F(IDF_CAUCHY, "IDF.CAUCHY", m_mdd_e, "a(0,1) c>0") \ - F(PDF_CAUCHY, "PDF.CAUCHY", m_mdd_e, "c>0") \ + F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \ + F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \ + F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \ + F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \ + F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \ + F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \ + F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0") \ F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \ - F(CDF_CHISQ, "CDF.CHISQ", m_md_e, "a>=0 b>0") \ - F(IDF_CHISQ, "IDF.CHISQ", m_md_e, "a[0,1) b>0") \ - F(PDF_CHISQ, "PDF.CHISQ", m_md_e, "a>=0 b>0") \ + F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \ + F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \ + F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \ + F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0") \ F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \ - F(SIG_CHISQ, "SIG.CHISQ", m_md_e, "a>=0 b>0") \ - F(CDF_EXP, "CDF.EXP", m_md_e, "a>=0 b>=0") \ - F(IDF_EXP, "IDF.EXP", m_md_e, "a[0,1) b>0") \ - F(PDF_EXP, "PDF.EXP", m_md_e, "a>=0 b>0") \ + F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \ + F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \ + F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \ + F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0") \ F(RV_EXP, "RV.EXP", d_d, "a>0") \ - F(PDF_XPOWER, "PDF.XPOWER", m_mdd_e, "b>0 c>=0") \ + F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \ F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \ - F(CDF_F, "CDF.F", m_mdd_e, "a>=0 b>0 c>0") \ - F(IDF_F, "IDF.F", m_mdd_e, "a[0,1) b>0 c>0") \ - F(PDF_F, "PDF.F", m_mdd_e, "a>=0 b>0 c>0") \ + F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \ + F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \ + F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \ + F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \ F(RV_F, "RV.F", d_dd, "a>0 b>0") \ - F(SIG_F, "SIG.F", m_mdd_e, "a>=0 b>0 c>0") + F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \ + F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \ + F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \ + F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \ + F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \ + F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \ + F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \ + F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \ + F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \ + F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \ + F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \ + F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \ + F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \ + F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \ + F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \ + F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \ + F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \ + F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \ + F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \ + F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \ + F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \ + F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \ + F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \ + F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \ + F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \ + F(CDFNORM, "CDFNORM", m_e, NULL) \ + F(PROBIT, "PROBIT", m_e, "a(0,1)") \ + F(NORMAL, "NORMAL", m_e, "a>0") \ + F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \ + F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \ + F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \ + F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \ + F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \ + F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \ + F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \ + F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \ + F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \ + F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \ + F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \ + F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \ + F(CDF_T, "CDF.T", m_ed, "b>0") \ + F(TCDF, "TCDF", m_ed, "b>0") \ + F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \ + F(PDF_T, "PDF.T", m_ed, "b>0") \ + F(RV_T, "RV.T", d_d, "a>0") \ + F(CDF_T1G, "CDF.T1G", m_edd, NULL) \ + F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \ + F(PDF_T1G, "PDF.T1G", m_edd, NULL) \ + F(RV_T1G, "RV.T1G", d_dd, NULL) \ + F(CDF_T2G, "CDF.T2G", m_edd, NULL) \ + F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \ + F(PDF_T2G, "PDF.T2G", m_edd, NULL) \ + F(RV_T2G, "RV.T2G", d_dd, NULL) \ + F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \ + F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \ + F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \ + F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \ + F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \ + F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \ + F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \ + F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \ + F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \ + F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \ + F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \ + F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \ + F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \ + F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \ + F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \ + F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \ + F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \ + F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \ + F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \ + F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \ + F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \ + F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \ + F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \ + F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \ + F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \ + F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \ + F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \ + F(RV_POISSON, "RV.POISSON", d_d, "a>0") struct matrix_function_properties { @@ -270,94 +352,43 @@ struct matrix_function_properties const char *constraints; }; -/* -() -(a) -(a > 0) -(a >= 0) -(a > 0 && a < 1) -(a > 0 && a <= 1) -(a >= 0 && a <= 1) -(a > 0 && a < 1, b > 0) -(a >= 0 && a < 1, b > 0) -(a >= 0 && a <= 1, b > 0) -(a == 0 || a == 1, b >= 0 && b <= 1) -(a >= 0 && a < 1, b > 0, c > 0) -(a >= 0 && a <= 1, b > 0, c > 0) -(a >= 0 && a < 1, b >= 1, c >= 1) -(a >= 0 && a <= 1, b, c) -(a > 0 && a < 1, b, c > 0) -(a >= 0 && a <= 1, b <= c, c) -(a > 0 && a == floor (a), -(a >= 0 && a == floor (a) && a <= b, -(a >= 0 && a == floor (a) && a <= d, -(a >= 0 && a == floor (a), b > 0) -(a > 0 && a == floor (a), b >= 0 && b <= 1) -(a > 0, b > 0) -(a > 0, b >= 0) -(a >= 0, b > 0) -(a >= 0, b > 0, c) -(a > 0, b > 0, c > 0) -(a >= 0, b > 0, c > 0) -(a >= 0, b > 0, c > 0, d > 0) -(a >= 0, b > 0, c > 0, d >= 0) -(a > 0, b >= 1, c >= 1) -(a >= 1 && a == floor (a), b >= 0 && b <= 1) -(a >= 1, b > 0 && b <= 1) -(a >= 1, b == floor (b), c > 0 && c <= 1) -(a, b) -(a, b > 0) -(a, b > 0 && b <= 2) -(a, b > 0 && b <= 2, c >= -1 && c <= 1) -(a, b > 0 && b == floor (c), c >= 0 && c <= 1) -(a, b > 0, c) -(a, b > 0, c > 0) -(a, b > 0, c >= 0) -(a <= b, b) -(a >= b, b > 0, c > 0) -(a, b, c) -(a, b, c > 0) -(a, b, c >= -1 && c <= 1) -(a <= c, b <= a, c) -(a == floor (a), b > 0 && b <= 1) -b > 0 && b == floor (b), -b > 0 && b == floor (b) && b <= a, -c >= 0 && c <= 1) -c > 0 && c == floor (c) && c <= a) -c > 0 && c == floor (c) && c <= b, -d > 0 && d == floor (d) && d <= b) -*/ - +enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 }; enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 }; enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 }; +enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 }; enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 }; enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 }; enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 }; enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 }; -enum { m_m_e_MIN_ARGS = 1, m_m_e_MAX_ARGS = 1 }; +enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 }; enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 }; -enum { m_md_e_MIN_ARGS = 2, m_md_e_MAX_ARGS = 2 }; +enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 }; enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 }; -enum { m_mdd_e_MIN_ARGS = 3, m_mdd_e_MAX_ARGS = 3 }; -enum { m_mddd_e_MIN_ARGS = 4, m_mddd_e_MAX_ARGS = 4 }; +enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 }; +enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 }; +enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 }; enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 }; enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 }; enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 }; enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX }; enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 }; +typedef double matrix_proto_d_none (void); typedef double matrix_proto_d_d (double); typedef double matrix_proto_d_dd (double, double); +typedef double matrix_proto_d_dd (double, double); +typedef double matrix_proto_d_ddd (double, double, double); typedef gsl_matrix *matrix_proto_m_d (double); typedef gsl_matrix *matrix_proto_m_dd (double, double); typedef gsl_matrix *matrix_proto_m_ddd (double, double, double); typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *); -typedef double matrix_proto_m_m_e (double); +typedef double matrix_proto_m_e (double); typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double); -typedef double matrix_proto_m_md_e (double, double); +typedef double matrix_proto_m_ed (double, double); typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double); -typedef double matrix_proto_m_mdd_e (double, double, double); -typedef double matrix_proto_m_mddd_e (double, double, double, double); +typedef double matrix_proto_m_edd (double, double, double); +typedef double matrix_proto_m_eddd (double, double, double, double); +typedef double matrix_proto_m_eed (double, double, double); typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *); typedef gsl_matrix *matrix_proto_m_v (gsl_vector *); typedef double matrix_proto_d_m (gsl_matrix *); @@ -1781,6 +1812,12 @@ matrix_eval_NPDF_BETA (double x, double a, double b, double lambda) return npdf_beta (x, a, b, lambda); } +static double +matrix_eval_CDF_BVNOR (double x0, double x1, double r) +{ + return cdf_bvnor (x0, x1, r); +} + static double matrix_eval_PDF_BVNOR (double x0, double x1, double r) { @@ -1817,6 +1854,12 @@ matrix_eval_CDF_CHISQ (double x, double df) return gsl_cdf_chisq_P (x, df); } +static double +matrix_eval_CHICDF (double x, double df) +{ + return matrix_eval_CDF_CHISQ (x, df); +} + static double matrix_eval_IDF_CHISQ (double P, double df) { @@ -1883,6 +1926,12 @@ matrix_eval_CDF_F (double x, double df1, double df2) return gsl_cdf_fdist_P (x, df1, df2); } +static double +matrix_eval_FCDF (double x, double df1, double df2) +{ + return matrix_eval_CDF_F (x, df1, df2); +} + static double matrix_eval_IDF_F (double P, double df1, double df2) { @@ -1907,6 +1956,486 @@ matrix_eval_SIG_F (double x, double df1, double df2) return gsl_cdf_fdist_Q (x, df1, df2); } +static double +matrix_eval_CDF_GAMMA (double x, double a, double b) +{ + return gsl_cdf_gamma_P (x, a, 1. / b); +} + +static double +matrix_eval_IDF_GAMMA (double P, double a, double b) +{ + return gsl_cdf_gamma_Pinv (P, a, 1. / b); +} + +static double +matrix_eval_PDF_GAMMA (double x, double a, double b) +{ + return gsl_ran_gamma_pdf (x, a, 1. / b); +} + +static double +matrix_eval_RV_GAMMA (double a, double b) +{ + return gsl_ran_gamma (get_rng (), a, 1. / b); +} + +static double +matrix_eval_PDF_LANDAU (double x) +{ + return gsl_ran_landau_pdf (x); +} + +static double +matrix_eval_RV_LANDAU (void) +{ + return gsl_ran_landau (get_rng ()); +} + +static double +matrix_eval_CDF_LAPLACE (double x, double a, double b) +{ + return gsl_cdf_laplace_P ((x - a) / b, 1); +} + +static double +matrix_eval_IDF_LAPLACE (double P, double a, double b) +{ + return a + b * gsl_cdf_laplace_Pinv (P, 1); +} + +static double +matrix_eval_PDF_LAPLACE (double x, double a, double b) +{ + return gsl_ran_laplace_pdf ((x - a) / b, 1); +} + +static double +matrix_eval_RV_LAPLACE (double a, double b) +{ + return a + b * gsl_ran_laplace (get_rng (), 1); +} + +static double +matrix_eval_RV_LEVY (double c, double alpha) +{ + return gsl_ran_levy (get_rng (), c, alpha); +} + +static double +matrix_eval_RV_LVSKEW (double c, double alpha, double beta) +{ + return gsl_ran_levy_skew (get_rng (), c, alpha, beta); +} + +static double +matrix_eval_CDF_LOGISTIC (double x, double a, double b) +{ + return gsl_cdf_logistic_P ((x - a) / b, 1); +} + +static double +matrix_eval_IDF_LOGISTIC (double P, double a, double b) +{ + return a + b * gsl_cdf_logistic_Pinv (P, 1); +} + +static double +matrix_eval_PDF_LOGISTIC (double x, double a, double b) +{ + return gsl_ran_logistic_pdf ((x - a) / b, 1) / b; +} + +static double +matrix_eval_RV_LOGISTIC (double a, double b) +{ + return a + b * gsl_ran_logistic (get_rng (), 1); +} + +static double +matrix_eval_CDF_LNORMAL (double x, double m, double s) +{ + return gsl_cdf_lognormal_P (x, log (m), s); +} + +static double +matrix_eval_IDF_LNORMAL (double P, double m, double s) +{ + return gsl_cdf_lognormal_Pinv (P, log (m), s);; +} + +static double +matrix_eval_PDF_LNORMAL (double x, double m, double s) +{ + return gsl_ran_lognormal_pdf (x, log (m), s); +} + +static double +matrix_eval_RV_LNORMAL (double m, double s) +{ + return gsl_ran_lognormal (get_rng (), log (m), s); +} + +static double +matrix_eval_CDF_NORMAL (double x, double u, double s) +{ + return gsl_cdf_gaussian_P (x - u, s); +} + +static double +matrix_eval_IDF_NORMAL (double P, double u, double s) +{ + return u + gsl_cdf_gaussian_Pinv (P, s); +} + +static double +matrix_eval_PDF_NORMAL (double x, double u, double s) +{ + return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s; +} + +static double +matrix_eval_RV_NORMAL (double u, double s) +{ + return u + gsl_ran_gaussian (get_rng (), s); +} + +static double +matrix_eval_CDFNORM (double x) +{ + return gsl_cdf_ugaussian_P (x); +} + +static double +matrix_eval_PROBIT (double P) +{ + return gsl_cdf_ugaussian_Pinv (P); +} + +static double +matrix_eval_NORMAL (double s) +{ + return gsl_ran_gaussian (get_rng (), s); +} + +static double +matrix_eval_PDF_NTAIL (double x, double a, double sigma) +{ + return gsl_ran_gaussian_tail_pdf (x, a, sigma);; +} + +static double +matrix_eval_RV_NTAIL (double a, double sigma) +{ + return gsl_ran_gaussian_tail (get_rng (), a, sigma); +} + +static double +matrix_eval_CDF_PARETO (double x, double a, double b) +{ + return gsl_cdf_pareto_P (x, b, a); +} + +static double +matrix_eval_IDF_PARETO (double P, double a, double b) +{ + return gsl_cdf_pareto_Pinv (P, b, a); +} + +static double +matrix_eval_PDF_PARETO (double x, double a, double b) +{ + return gsl_ran_pareto_pdf (x, b, a); +} + +static double +matrix_eval_RV_PARETO (double a, double b) +{ + return gsl_ran_pareto (get_rng (), b, a); +} + +static double +matrix_eval_CDF_RAYLEIGH (double x, double sigma) +{ + return gsl_cdf_rayleigh_P (x, sigma); +} + +static double +matrix_eval_IDF_RAYLEIGH (double P, double sigma) +{ + return gsl_cdf_rayleigh_Pinv (P, sigma); +} + +static double +matrix_eval_PDF_RAYLEIGH (double x, double sigma) +{ + return gsl_ran_rayleigh_pdf (x, sigma); +} + +static double +matrix_eval_RV_RAYLEIGH (double sigma) +{ + return gsl_ran_rayleigh (get_rng (), sigma); +} + +static double +matrix_eval_PDF_RTAIL (double x, double a, double sigma) +{ + return gsl_ran_rayleigh_tail_pdf (x, a, sigma); +} + +static double +matrix_eval_RV_RTAIL (double a, double sigma) +{ + return gsl_ran_rayleigh_tail (get_rng (), a, sigma); +} + +static double +matrix_eval_CDF_T (double x, double df) +{ + return gsl_cdf_tdist_P (x, df); +} + +static double +matrix_eval_TCDF (double x, double df) +{ + return matrix_eval_CDF_T (x, df); +} + +static double +matrix_eval_IDF_T (double P, double df) +{ + return gsl_cdf_tdist_Pinv (P, df); +} + +static double +matrix_eval_PDF_T (double x, double df) +{ + return gsl_ran_tdist_pdf (x, df); +} + +static double +matrix_eval_RV_T (double df) +{ + return gsl_ran_tdist (get_rng (), df); +} + +static double +matrix_eval_CDF_T1G (double x, double a, double b) +{ + return gsl_cdf_gumbel1_P (x, a, b); +} + +static double +matrix_eval_IDF_T1G (double P, double a, double b) +{ + return gsl_cdf_gumbel1_Pinv (P, a, b); +} + +static double +matrix_eval_PDF_T1G (double x, double a, double b) +{ + return gsl_ran_gumbel1_pdf (x, a, b); +} + +static double +matrix_eval_RV_T1G (double a, double b) +{ + return gsl_ran_gumbel1 (get_rng (), a, b); +} + +static double +matrix_eval_CDF_T2G (double x, double a, double b) +{ + return gsl_cdf_gumbel1_P (x, a, b); +} + +static double +matrix_eval_IDF_T2G (double P, double a, double b) +{ + return gsl_cdf_gumbel1_Pinv (P, a, b); +} + +static double +matrix_eval_PDF_T2G (double x, double a, double b) +{ + return gsl_ran_gumbel1_pdf (x, a, b); +} + +static double +matrix_eval_RV_T2G (double a, double b) +{ + return gsl_ran_gumbel1 (get_rng (), a, b); +} + +static double +matrix_eval_CDF_UNIFORM (double x, double a, double b) +{ + return gsl_cdf_flat_P (x, a, b); +} + +static double +matrix_eval_IDF_UNIFORM (double P, double a, double b) +{ + return gsl_cdf_flat_Pinv (P, a, b); +} + +static double +matrix_eval_PDF_UNIFORM (double x, double a, double b) +{ + return gsl_ran_flat_pdf (x, a, b); +} + +static double +matrix_eval_RV_UNIFORM (double a, double b) +{ + return gsl_ran_flat (get_rng (), a, b); +} + +static double +matrix_eval_CDF_WEIBULL (double x, double a, double b) +{ + return gsl_cdf_weibull_P (x, a, b); +} + +static double +matrix_eval_IDF_WEIBULL (double P, double a, double b) +{ + return gsl_cdf_weibull_Pinv (P, a, b); +} + +static double +matrix_eval_PDF_WEIBULL (double x, double a, double b) +{ + return gsl_ran_weibull_pdf (x, a, b); +} + +static double +matrix_eval_RV_WEIBULL (double a, double b) +{ + return gsl_ran_weibull (get_rng (), a, b); +} + +static double +matrix_eval_CDF_BERNOULLI (double k, double p) +{ + return k ? 1 : 1 - p; +} + +static double +matrix_eval_PDF_BERNOULLI (double k, double p) +{ + return gsl_ran_bernoulli_pdf (k, p); +} + +static double +matrix_eval_RV_BERNOULLI (double p) +{ + return gsl_ran_bernoulli (get_rng (), p); +} + +static double +matrix_eval_CDF_BINOM (double k, double n, double p) +{ + return gsl_cdf_binomial_P (k, p, n); +} + +static double +matrix_eval_PDF_BINOM (double k, double n, double p) +{ + return gsl_ran_binomial_pdf (k, p, n); +} + +static double +matrix_eval_RV_BINOM (double n, double p) +{ + return gsl_ran_binomial (get_rng (), p, n); +} + +static double +matrix_eval_CDF_GEOM (double k, double p) +{ + return gsl_cdf_geometric_P (k, p); +} + +static double +matrix_eval_PDF_GEOM (double k, double p) +{ + return gsl_ran_geometric_pdf (k, p); +} + +static double +matrix_eval_RV_GEOM (double p) +{ + return gsl_ran_geometric (get_rng (), p); +} + +static double +matrix_eval_CDF_HYPER (double k, double a, double b, double c) +{ + return gsl_cdf_hypergeometric_P (k, c, a - c, b); +} + +static double +matrix_eval_PDF_HYPER (double k, double a, double b, double c) +{ + return gsl_ran_hypergeometric_pdf (k, c, a - c, b); +} + +static double +matrix_eval_RV_HYPER (double a, double b, double c) +{ + return gsl_ran_hypergeometric (get_rng (), c, a - c, b); +} + +static double +matrix_eval_PDF_LOG (double k, double p) +{ + return gsl_ran_logarithmic_pdf (k, p); +} + +static double +matrix_eval_RV_LOG (double p) +{ + return gsl_ran_logarithmic (get_rng (), p); +} + +static double +matrix_eval_CDF_NEGBIN (double k, double n, double p) +{ + return gsl_cdf_negative_binomial_P (k, p, n); +} + +static double +matrix_eval_PDF_NEGBIN (double k, double n, double p) +{ + return gsl_ran_negative_binomial_pdf (k, p, n); +} + +static double +matrix_eval_RV_NEGBIN (double n, double p) +{ + return gsl_ran_negative_binomial (get_rng (), p, n); +} + +static double +matrix_eval_CDF_POISSON (double k, double mu) +{ + return gsl_cdf_poisson_P (k, mu); +} + +static double +matrix_eval_PDF_POISSON (double k, double mu) +{ + return gsl_ran_poisson_pdf (k, mu); +} + +static double +matrix_eval_RV_POISSON (double mu) +{ + return gsl_ran_poisson (get_rng (), mu); +} + struct matrix_function { const char *name; @@ -3065,10 +3594,80 @@ parse_constraint_value (const char **constraintsp) { char *tail; long retval = strtol (*constraintsp, &tail, 10); + assert (tail > *constraintsp); *constraintsp = tail; return retval; } +static void +argument_invalid (const struct matrix_function_properties *props, + size_t arg_index, gsl_matrix *a, size_t y, size_t x, + struct string *out) +{ + if (is_scalar (a)) + ds_put_format (out, _("Argument %zu to matrix function %s " + "has invalid value %g."), + arg_index, props->name, gsl_matrix_get (a, y, x)); + else + ds_put_format (out, _("Row %zu, column %zu of argument %zu " + "to matrix function %s has " + "invalid value %g."), + y + 1, x + 1, arg_index, props->name, + gsl_matrix_get (a, y, x)); +} + +static void +argument_inequality_error (const struct matrix_function_properties *props, + size_t a_index, gsl_matrix *a, size_t y, size_t x, + size_t b_index, double b, + bool greater, bool equal) +{ + struct string s = DS_EMPTY_INITIALIZER; + + argument_invalid (props, a_index, a, y, x, &s); + ds_put_cstr (&s, " "); + if (greater && equal) + ds_put_format (&s, _("This argument must be greater than or " + "equal to argument %zu, but that has value %g."), + b_index, b); + else if (greater && !equal) + ds_put_format (&s, _("This argument must be greater than argument %zu, " + "but that has value %g."), + b_index, b); + else if (equal) + ds_put_format (&s, _("This argument must be less than or " + "equal to argument %zu, but that has value %g."), + b_index, b); + else + ds_put_format (&s, _("This argument must be less than argument %zu, " + "but that has value %g."), + b_index, b); + msg (ME, ds_cstr (&s)); + ds_destroy (&s); +} + +static void +argument_value_error (const struct matrix_function_properties *props, + size_t a_index, gsl_matrix *a, size_t y, size_t x, + double b, + bool greater, bool equal) +{ + struct string s = DS_EMPTY_INITIALIZER; + argument_invalid (props, a_index, a, y, x, &s); + ds_put_cstr (&s, " "); + if (greater && equal) + ds_put_format (&s, _("This argument must be greater than or equal to %g."), + b); + else if (greater && !equal) + ds_put_format (&s, _("This argument must be greater than %g."), b); + else if (equal) + ds_put_format (&s, _("This argument must be less than or equal to %g."), b); + else + ds_put_format (&s, _("This argument must be less than %g."), b); + msg (ME, ds_cstr (&s)); + ds_destroy (&s); +} + static bool check_constraints (const struct matrix_function_properties *props, gsl_matrix *args[], size_t n_args) @@ -3129,56 +3728,84 @@ check_constraints (const struct matrix_function_properties *props, } else equal = false; - int comparand = parse_constraint_value (&constraints); - assert (arg_index < n_args); - MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index]) - if (greater - ? (equal ? !(*d >= comparand) : !(*d > comparand)) - : (equal ? !(*d <= comparand) : !(*d < comparand))) - { - struct string s = DS_EMPTY_INITIALIZER; - if (!is_scalar (args[arg_index])) - ds_put_format (&s, _("Row %zu, column %zu of argument %zu " - "to matrix function %s has " - "invalid value %g."), - y + 1, x + 1, arg_index + 1, props->name, *d); - else - ds_put_format (&s, _("Argument %zu to matrix function %s " - "has invalid value %g."), - arg_index + 1, props->name, *d); - - ds_put_cstr (&s, " "); - if (greater && equal) - ds_put_format (&s, _("This argument must be greater than or " - "equal to %d."), comparand); - else if (greater && !equal) - ds_put_format (&s, _("This argument must be greater than %d."), - comparand); - else if (equal) - ds_put_format (&s, _("This argument must be less than or " - "equal to %d."), comparand); - else - ds_put_format (&s, _("This argument must be less than %d."), - comparand); - msg (ME, ds_cstr (&s)); - ds_destroy (&s); - return false; - } + if (*constraints >= 'a' && *constraints <= 'd') + { + size_t a_index = arg_index; + size_t b_index = *constraints - 'a'; + assert (a_index < n_args); + assert (b_index < n_args); + + /* We only support one of the two arguments being non-scalar. + It's easier to support only the first one being non-scalar, so + flip things around if it's the other way. */ + if (!is_scalar (args[b_index])) + { + assert (is_scalar (args[a_index])); + size_t tmp_index = a_index; + a_index = b_index; + b_index = tmp_index; + + greater = !greater; + } + + double b = to_scalar (args[b_index]); + MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index]) + if (greater + ? (equal ? !(*a >= b) : !(*a > b)) + : (equal ? !(*a <= b) : !(*a < b))) + { + argument_inequality_error (props, + a_index + 1, args[a_index], y, x, + b_index + 1, b, + greater, equal); + return false; + } + } + else + { + int comparand = parse_constraint_value (&constraints); + + MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index]) + if (greater + ? (equal ? !(*d >= comparand) : !(*d > comparand)) + : (equal ? !(*d <= comparand) : !(*d < comparand))) + { + argument_value_error (props, + arg_index + 1, args[arg_index], y, x, + comparand, + greater, equal); + return false; + } + } } else { assert (*constraints == ' '); constraints++; + arg_index = SIZE_MAX; } } return true; } +static gsl_matrix * +matrix_expr_evaluate_d_none ( + const struct matrix_function_properties *props UNUSED, + gsl_matrix *subs[] UNUSED, size_t n_subs, + matrix_proto_d_none *f) +{ + assert (n_subs == 0); + + gsl_matrix *m = gsl_matrix_alloc (1, 1); + gsl_matrix_set (m, 0, 0, f ()); + return m; +} + static gsl_matrix * matrix_expr_evaluate_d_d (const struct matrix_function_properties *props, - gsl_matrix *subs[], size_t n_subs, - matrix_proto_d_d *f) + gsl_matrix *subs[], size_t n_subs, + matrix_proto_d_d *f) { assert (n_subs == 1); @@ -3213,6 +3840,25 @@ matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props, return m; } +static gsl_matrix * +matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_d_ddd *f) +{ + assert (n_subs == 3); + + double d[3]; + if (!to_scalar_args (props->name, subs, n_subs, d)) + return NULL; + + if (!check_constraints (props, subs, n_subs)) + return NULL; + + gsl_matrix *m = gsl_matrix_alloc (1, 1); + gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2])); + return m; +} + static gsl_matrix * matrix_expr_evaluate_m_d (const struct matrix_function_properties *props, gsl_matrix *subs[], size_t n_subs, @@ -3256,9 +3902,9 @@ matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED, } static gsl_matrix * -matrix_expr_evaluate_m_m_e (const struct matrix_function_properties *props, +matrix_expr_evaluate_m_e (const struct matrix_function_properties *props, gsl_matrix *subs[], size_t n_subs, - matrix_proto_m_m_e *f) + matrix_proto_m_e *f) { assert (n_subs == 1); @@ -3282,9 +3928,9 @@ matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED } static gsl_matrix * -matrix_expr_evaluate_m_md_e (const struct matrix_function_properties *props, - gsl_matrix *subs[], size_t n_subs, - matrix_proto_m_md_e *f) +matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_ed *f) { assert (n_subs == 2); if (!check_scalar_arg (props->name, subs, 1)) @@ -3311,9 +3957,9 @@ matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSE } static gsl_matrix * -matrix_expr_evaluate_m_mdd_e (const struct matrix_function_properties *props, - gsl_matrix *subs[], size_t n_subs, - matrix_proto_m_mdd_e *f) +matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_edd *f) { assert (n_subs == 3); if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2)) @@ -3330,9 +3976,9 @@ matrix_expr_evaluate_m_mdd_e (const struct matrix_function_properties *props, } static gsl_matrix * -matrix_expr_evaluate_m_mddd_e (const struct matrix_function_properties *props, - gsl_matrix *subs[], size_t n_subs, - matrix_proto_m_mddd_e *f) +matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_eddd *f) { assert (n_subs == 4); for (size_t i = 1; i < 4; i++) @@ -3350,6 +3996,49 @@ matrix_expr_evaluate_m_mddd_e (const struct matrix_function_properties *props, return subs[0]; } +static gsl_matrix * +matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_eed *f) +{ + assert (n_subs == 3); + if (!check_scalar_arg (props->name, subs, 2)) + return NULL; + + if (!is_scalar (subs[0]) && !is_scalar (subs[1]) + && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2)) + { + msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and " + "%zu×%zu, but %s requires these arguments either to have " + "the same dimensions or for one of them to be a scalar."), + props->name, + subs[0]->size1, subs[0]->size2, + subs[1]->size1, subs[1]->size2, + props->name); + return NULL; + } + + if (!check_constraints (props, subs, n_subs)) + return NULL; + + double c = to_scalar (subs[2]); + + if (is_scalar (subs[0])) + { + double a = to_scalar (subs[0]); + MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1]) + *b = f (a, *b, c); + return subs[1]; + } + else + { + double b = to_scalar (subs[1]); + MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0]) + *a = f (*a, b, c); + return subs[0]; + } +} + static gsl_matrix * matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED, gsl_matrix *subs[], size_t n_subs,