From ab799b27b983ea4f20825fc37669c2e39f9b9770 Mon Sep 17 00:00:00 2001 From: Mario Rodriguez Date: Fri, 3 Mar 2017 23:02:42 +0100 Subject: [PATCH] Make use of unit_step(x) instead of conditional sentences on x in some continuous distributions. --- share/distrib/distrib.mac | 87 ++++++++++++++++------------------------------- 1 file changed, 30 insertions(+), 57 deletions(-) diff --git a/share/distrib/distrib.mac b/share/distrib/distrib.mac index 99c715f49..bfd4554a6 100644 --- a/share/distrib/distrib.mac +++ b/share/distrib/distrib.mac @@ -223,7 +223,7 @@ random_student_t(n,[num]) := pdf_noncentral_student_t(x,n,ncp) := if maybe(n > 0) = false then error("pdf_noncentral_student_t: degrees of freedom must be greater than 0") - elseif float(ncp)=0.0 + elseif sign(ncp) = 'zero then pdf_student_t(x,n) else n^(n/2) * factorial(n) * exp(-ncp^2/2) / (2^n * (n+x^2)^(n/2) * gamma(n/2)) * @@ -234,11 +234,13 @@ pdf_noncentral_student_t(x,n,ncp) := cdf_noncentral_student_t(x,n,ncp) := if maybe(n > 0) = false then error("cdf_noncentral_student_t: degrees of freedom must be greater than 0") - elseif float(ncp)=0.0 + elseif sign(ncp) = 'zero then cdf_student_t(x,n) - elseif numberp(x) and numberp(n) and numberp(ncp) - then ?cdfnt(float(x),float(n),float(ncp)) - else error("cdf_noncentral_student_t: need numeric arguments for approximate procedure") $ + else /* numerical approximation */ + block([fx: float(x), fn: float(n), fncp: float(ncp)], + if numberp(fx) and numberp(fn) and numberp(fncp) + then ?cdfnt(fx,fn,fncp) + else error("cdf_noncentral_student_t: need numeric arguments for approximate procedure")) $ /* R: qt(q,n,ncp) */ quantile_noncentral_student_t(q,n,ncp) := @@ -330,11 +332,9 @@ random_chi2(n,[num]) := pdf_noncentral_chi2(x,n,ncp) := if maybe(n > 0) = false then error("pdf_noncentral_chi2: degrees of freedom must be greater than 0") - elseif float(ncp)=0.0 + elseif sign(ncp) = 'zero then pdf_chi2(x,n) - elseif maybe(x >= 0) = false - then 0 - else 1/2 * exp(-(x+ncp)/2) * (x/ncp)^(n/4-1/2) * bessel_i(n/2-1, sqrt(x*ncp)) $ + else 1/2 * exp(-(x+ncp)/2) * (x/ncp)^(n/4-1/2) * bessel_i(n/2-1, sqrt(x*ncp)) * unit_step(x) $ /* R: pchisq(x,n,ncp) */ cdf_noncentral_chi2(x,n,ncp) := @@ -406,17 +406,13 @@ random_noncentral_chi2(n,ncp,[num]) := pdf_f(x,m,n) := if maybe(n > 0 and m > 0) = false then error("pdf_f: degrees of freedom must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) - then 0 - else gamma((m+n)/2)*(m/n)^(m/2)*x^(m/2-1)*(1+m*x/n)^(-(m+n)/2) / (gamma(m/2)*gamma(n/2)) $ + else gamma((m+n)/2)*(m/n)^(m/2)*x^(m/2-1)*(1+m*x/n)^(-(m+n)/2) / (gamma(m/2)*gamma(n/2)) * unit_step(x) $ /* R: pf(x,m,n) */ cdf_f(x,m,n) := if maybe(n > 0 and m > 0) = false then error("cdf_f: degrees of freedom must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) - then 0 - else 1 - beta_incomplete_regularized(n/2, m/2, n/(n+m*x)) $ + else (1 - beta_incomplete_regularized(n/2, m/2, n/(n+m*x))) * unit_step(x) $ /* R: qf(q,m,n) */ quantile_f(q,m,n) := @@ -501,14 +497,14 @@ random_exp(m,[num]) := pdf_lognormal(x,m,s) := if maybe(s > 0) = false then error("pdf_lognormal: parameter s must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) + elseif maybe(x > 0) = false then 0 else exp(-(log(x)-m)^2/(2*s^2))/(sqrt(2*%pi)*s*x) $ cdf_lognormal(x,m,s) := if maybe(s > 0) = false then error("cdf_lognormal: parameter s must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) + elseif maybe(x>0) = false then 0 else 1/2+erf((log(x)-m)/(s*sqrt(2)))/2 $ @@ -561,17 +557,13 @@ random_lognormal(m,s,[num]) := pdf_gamma(x,a,b) := if maybe(a > 0 and b > 0) = false then error("pdf_gamma: parameters a and b must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) - then 0 - else x^(a-1)*exp(-x/b)/(b^a*gamma(a)) $ + else x^(a-1)*exp(-x/b)/(b^a*gamma(a)) * unit_step(x) $ /* R: pgamma(x,a,1/b) */ cdf_gamma(x,a,b) := if maybe(a > 0 and b > 0) = false then error("cdf_gamma: parameters a and b must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) - then 0 - else 1 - gamma_incomplete_regularized(a,x/b) $ + else (1 - gamma_incomplete_regularized(a,x/b)) * unit_step(x) $ /* R: qgamma(q,a,1/b) */ quantile_gamma(q,a,b) := @@ -581,12 +573,13 @@ quantile_gamma(q,a,b) := then 0 elseif sign(q-1)='zero then 'inf - else block([fq: float(q), fa: float(a)], - if numberp(fq) and numberp(fa) - then(if fq = 0.0 then return(0), - if fq = 1.0 then return('inf), - b * ?iigamma(fq, fa)) - else error("quantile_gamma: need numeric arguments for approximate procedure")) $ + else /* approximate procedure */ + block([fq: float(q), fa: float(a)], + if numberp(fq) and numberp(fa) + then(if fq = 0.0 then return(0), + if fq = 1.0 then return('inf), + b * ?iigamma(fq, fa)) + else error("quantile_gamma: need numeric arguments for approximate procedure")) $ mean_gamma(a,b) := if maybe(a > 0 and b > 0) = false @@ -628,19 +621,13 @@ random_gamma(a,b,[num]) := pdf_beta(x,a,b) := if maybe(a > 0 and b > 0) = false then error("pdf_beta: parameters a and b must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) or member(sign(x-1), ['pos,'pz,'zero]) - then 0 - else x^(a-1)*(1-x)^(b-1)/beta(a,b) $ + else x^(a-1)*(1-x)^(b-1)/beta(a,b)*(unit_step(x)-unit_step(x-1)) $ /* R: pbeta(x,a,b) */ cdf_beta(x,a,b) := if maybe(a > 0 and b > 0) = false then error("cdf_beta: parameters a and b must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) - then 0 - elseif member(sign(x-1), ['pos,'pz,'zero]) - then 1 - else beta_incomplete_regularized(a,b,x) $ + else beta_incomplete_regularized(a,b,x) * (unit_step(x)-unit_step(x-1)) + unit_step(x-1)$ /* R: qbeta(q,a,b) */ quantile_beta(q,a,b) := @@ -697,18 +684,12 @@ random_beta(a,b,[num]) := pdf_continuous_uniform(x,a,b) := if maybe(b - a > 0) = false then error("pdf_continuous_uniform: parameter b must be greater than a") - elseif maybe(x >=a and x <= b) = false - then 0 - else 1/(b-a) $ + else (unit_step(x-a)-unit_step(x-b))/(b-a) $ cdf_continuous_uniform(x,a,b) := if maybe(b - a > 0) = false then error("cdf_continuous_uniform: parameter b must be greater than a") - elseif maybe(x >= a) = false - then 0 - elseif maybe(x <= b) = false - then 1 - else (x-a)/(b-a) $ + else (x-a)/(b-a)*(unit_step(x-a)-unit_step(x-b)) + unit_step(x-b) $ quantile_continuous_uniform(q,a,b) := if maybe(b - a > 0 and q >= 0 and q <= 1) = false @@ -822,16 +803,12 @@ random_logistic(a,b,[num]) := pdf_pareto(x,a,b) := if maybe(a > 0 and b > 0) = false then error("pdf_pareto: parameters a and b must be greater than 0") - elseif member(sign(x-b), ['nz,'neg,'zero]) - then 0 - else a*b^a*x^(-a-1) $ + else a*b^a*x^(-a-1) * unit_step(x-b) $ cdf_pareto(x,a,b) := if maybe(a > 0 and b > 0) = false then error("pdf_pareto: parameters a and b must be greater than 0") - elseif member(sign(x-b), ['nz,'neg,'zero]) - then 0 - else 1-(b/x)^a $ + else (1-(b/x)^a) * unit_step(x-b) $ quantile_pareto(q,a,b) := if maybe(a > 0 and b > 0 and q >= 0 and q <= 1) = false @@ -898,16 +875,12 @@ random_pareto(a,b,[num]) := pdf_weibull(x,a,b) := if maybe(a > 0 and b > 0) = false then error("pdf_weibull: parameters a and b must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) - then 0 - else a/b*(x/b)^(a-1)*exp(-(x/b)^a) $ + else a/b*(x/b)^(a-1)*exp(-(x/b)^a) * unit_step(x) $ cdf_weibull(x,a,b) := if maybe(a > 0 and b > 0) = false then error("cdf_weibull: parameters a and b must be greater than 0") - elseif member(sign(x), ['nz,'neg,'zero]) - then 0 - else 1-exp(-(x/b)^a) $ + else (1-exp(-(x/b)^a)) * unit_step(x) $ /* R: qweibull(q,a,b) */ quantile_weibull(q,a,b) := -- 2.11.4.GIT