You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@systemml.apache.org by de...@apache.org on 2017/07/14 17:35:47 UTC
systemml git commit: [SYSTEMML-1419] Cleanup nested if-elses in GLM
and GLM-predict
Repository: systemml
Updated Branches:
refs/heads/master a4ce06461 -> ccac6dd37
[SYSTEMML-1419] Cleanup nested if-elses in GLM and GLM-predict
Closes #562.
Project: http://git-wip-us.apache.org/repos/asf/systemml/repo
Commit: http://git-wip-us.apache.org/repos/asf/systemml/commit/ccac6dd3
Tree: http://git-wip-us.apache.org/repos/asf/systemml/tree/ccac6dd3
Diff: http://git-wip-us.apache.org/repos/asf/systemml/diff/ccac6dd3
Branch: refs/heads/master
Commit: ccac6dd37de3dba745069caeab9803aa5d69190f
Parents: a4ce064
Author: Janardhan <pu...@nitkkr.ac.in>
Authored: Fri Jul 14 10:29:55 2017 -0700
Committer: Deron Eriksson <de...@apache.org>
Committed: Fri Jul 14 10:29:55 2017 -0700
----------------------------------------------------------------------
scripts/algorithms/GLM-predict.dml | 40 +++----
scripts/algorithms/GLM.dml | 202 ++++++++++++++++----------------
2 files changed, 123 insertions(+), 119 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/systemml/blob/ccac6dd3/scripts/algorithms/GLM-predict.dml
----------------------------------------------------------------------
diff --git a/scripts/algorithms/GLM-predict.dml b/scripts/algorithms/GLM-predict.dml
index 251d85a..49a593c 100644
--- a/scripts/algorithms/GLM-predict.dml
+++ b/scripts/algorithms/GLM-predict.dml
@@ -105,13 +105,13 @@ dispersion = as.double (dispersion);
if (dist_type == 3) {
link_type = 2;
-} else { if (link_type == 0) { # Canonical Link
+} else if (link_type == 0) { # Canonical Link
if (dist_type == 1) {
link_type = 1;
link_power = 1.0 - var_power;
- } else { if (dist_type == 2) {
+ } else if (dist_type == 2) {
link_type = 2;
-}}} }
+} }
X = read (fileX);
num_records = nrow (X);
@@ -343,36 +343,36 @@ glm_means_and_vars =
# POWER DISTRIBUTION
if (link_power == 0) {
y_mean = exp (linear_terms);
- } else { if (link_power == 1.0) {
+ } else if (link_power == 1.0) {
y_mean = linear_terms;
- } else { if (link_power == -1.0) {
+ } else if (link_power == -1.0) {
y_mean = 1.0 / linear_terms;
} else {
y_mean = linear_terms ^ (1.0 / link_power);
- }}}
+ }
if (var_power == 0) {
var_function = matrix (1.0, rows = num_points, cols = 1);
- } else { if (var_power == 1.0) {
+ } else if (var_power == 1.0) {
var_function = y_mean;
} else {
var_function = y_mean ^ var_power;
- }}
+ }
means = y_mean;
vars = var_function;
- } else { if (dist_type == 2 & link_type >= 1 & link_type <= 5) {
+ } else if (dist_type == 2 & link_type >= 1 & link_type <= 5) {
# BINOMIAL/BERNOULLI DISTRIBUTION
y_prob = matrix (0.0, rows = num_points, cols = 2);
if (link_type == 1 & link_power == 0) { # Binomial.log
y_prob [, 1] = exp (linear_terms);
y_prob [, 2] = 1.0 - y_prob [, 1];
- } else { if (link_type == 1 & link_power != 0) { # Binomial.power_nonlog
+ } else if (link_type == 1 & link_power != 0) { # Binomial.power_nonlog
y_prob [, 1] = linear_terms ^ (1.0 / link_power);
y_prob [, 2] = 1.0 - y_prob [, 1];
- } else { if (link_type == 2) { # Binomial.logit
+ } else if (link_type == 2) { # Binomial.logit
elt = exp (linear_terms);
y_prob [, 1] = elt / (1.0 + elt);
y_prob [, 2] = 1.0 / (1.0 + elt);
- } else { if (link_type == 3) { # Binomial.probit
+ } else if (link_type == 3) { # Binomial.probit
sign_lt = 2 * (linear_terms >= 0) - 1;
t_gp = 1.0 / (1.0 + abs (linear_terms) * 0.231641888); # 0.231641888 = 0.3275911 / sqrt (2.0)
erf_corr =
@@ -384,20 +384,20 @@ glm_means_and_vars =
y_prob [, 1] = (1 + sign_lt) - erf_corr;
y_prob [, 2] = (1 - sign_lt) + erf_corr;
y_prob = y_prob / 2;
- } else { if (link_type == 4) { # Binomial.cloglog
+ } else if (link_type == 4) { # Binomial.cloglog
elt = exp (linear_terms);
is_too_small = ((10000000 + elt) == 10000000);
y_prob [, 2] = exp (- elt);
y_prob [, 1] = (1 - is_too_small) * (1.0 - y_prob [, 2]) + is_too_small * elt * (1.0 - elt / 2);
- } else { if (link_type == 5) { # Binomial.cauchit
+ } else if (link_type == 5) { # Binomial.cauchit
atan_linear_terms = atan (linear_terms);
y_prob [, 1] = 0.5 + atan_linear_terms / 3.1415926535897932384626433832795;
y_prob [, 2] = 0.5 - atan_linear_terms / 3.1415926535897932384626433832795;
- }}}}}}
+ }
means = y_prob;
ones_ctg = matrix (1, rows = 2, cols = 1);
vars = means * (means %*% (1 - diag (ones_ctg)));
- } else { if (dist_type == 3) {
+ } else if (dist_type == 3) {
# MULTINOMIAL LOGIT DISTRIBUTION
elt = exp (linear_terms);
ones_pts = matrix (1, rows = num_points, cols = 1);
@@ -408,7 +408,7 @@ glm_means_and_vars =
} else {
means = matrix (0.0, rows = num_points, cols = 1);
vars = matrix (0.0, rows = num_points, cols = 1);
-} }}}
+} }
glm_partial_loglikelihood_for_power_dist_and_link = # Assumes: dist_type == 1 & link_type == 1
function (Matrix[double] linear_terms, Matrix[double] Y, double var_power, double link_power)
@@ -431,10 +431,10 @@ glm_partial_loglikelihood_for_power_dist_and_link = # Assumes: dist_type == 1
if (var_power == 2.0 & link_power == 0) { # Gamma.log
natural_parameters = - exp (- linear_terms);
b_cumulant = linear_terms;
- } else { if (var_power == 2.0) { # Gamma.power_nonlog
+ } else if (var_power == 2.0) { # Gamma.power_nonlog
natural_parameters = - linear_terms ^ (- 1.0 / link_power);
b_cumulant = log (linear_terms) / link_power;
- } else { if (link_power == 0) { # PowerDist.log
+ } else if (link_power == 0) { # PowerDist.log
natural_parameters = exp (linear_terms * (1.0 - var_power)) / (1.0 - var_power);
b_cumulant = exp (linear_terms * (2.0 - var_power)) / (2.0 - var_power);
} else { # PowerDist.power_nonlog
@@ -442,6 +442,6 @@ glm_partial_loglikelihood_for_power_dist_and_link = # Assumes: dist_type == 1
natural_parameters = (linear_terms ^ power_np) / (1.0 - var_power);
power_cu = (2.0 - var_power) / link_power;
b_cumulant = (linear_terms ^ power_cu) / (2.0 - var_power);
- }}}
+ }
log_l_part = Y * natural_parameters - b_cumulant;
} }
http://git-wip-us.apache.org/repos/asf/systemml/blob/ccac6dd3/scripts/algorithms/GLM.dml
----------------------------------------------------------------------
diff --git a/scripts/algorithms/GLM.dml b/scripts/algorithms/GLM.dml
index 2524d71..db911d6 100644
--- a/scripts/algorithms/GLM.dml
+++ b/scripts/algorithms/GLM.dml
@@ -251,9 +251,9 @@ if (link_type == 0)
if (distribution_type == 1) {
link_type = 1;
link_as_power_of_the_mean = 1.0 - variance_as_power_of_the_mean;
- } else { if (distribution_type == 2) {
+ } else if (distribution_type == 2) {
link_type = 2;
-} } }
+} }
# For power distributions and/or links, we use two constants,
# "variance as power of the mean" and "link_as_power_of_the_mean",
@@ -519,30 +519,30 @@ check_if_supported =
if (ncol_y == 1 & dist_type == 1 & link_type == 1)
{ # POWER DISTRIBUTION
is_supported = 1;
- if (var_power == 0 & link_power == -1.0) {print ("Gaussian.inverse"); } else {
- if (var_power == 0 & link_power == 0) {print ("Gaussian.log"); } else {
- if (var_power == 0 & link_power == 0.5) {print ("Gaussian.sqrt"); } else {
- if (var_power == 0 & link_power == 1.0) {print ("Gaussian.id"); } else {
- if (var_power == 0 ) {print ("Gaussian.power_nonlog"); } else {
- if (var_power == 1.0 & link_power == -1.0) {print ("Poisson.inverse"); } else {
- if (var_power == 1.0 & link_power == 0) {print ("Poisson.log"); } else {
- if (var_power == 1.0 & link_power == 0.5) {print ("Poisson.sqrt"); } else {
- if (var_power == 1.0 & link_power == 1.0) {print ("Poisson.id"); } else {
- if (var_power == 1.0 ) {print ("Poisson.power_nonlog"); } else {
- if (var_power == 2.0 & link_power == -1.0) {print ("Gamma.inverse"); } else {
- if (var_power == 2.0 & link_power == 0) {print ("Gamma.log"); } else {
- if (var_power == 2.0 & link_power == 0.5) {print ("Gamma.sqrt"); } else {
- if (var_power == 2.0 & link_power == 1.0) {print ("Gamma.id"); } else {
- if (var_power == 2.0 ) {print ("Gamma.power_nonlog"); } else {
- if (var_power == 3.0 & link_power == -2.0) {print ("InvGaussian.1/mu^2"); } else {
- if (var_power == 3.0 & link_power == -1.0) {print ("InvGaussian.inverse"); } else {
- if (var_power == 3.0 & link_power == 0) {print ("InvGaussian.log"); } else {
- if (var_power == 3.0 & link_power == 0.5) {print ("InvGaussian.sqrt"); } else {
- if (var_power == 3.0 & link_power == 1.0) {print ("InvGaussian.id"); } else {
- if (var_power == 3.0 ) {print ("InvGaussian.power_nonlog");}else{
- if ( link_power == 0) {print ("PowerDist.log"); } else {
- print ("PowerDist.power_nonlog");
- } }}}}} }}}}} }}}}} }}}}} }}
+ if (var_power == 0 & link_power == -1.0) print ("Gaussian.inverse");
+ else if (var_power == 0 & link_power == 0) print ("Gaussian.log");
+ else if (var_power == 0 & link_power == 0.5) print ("Gaussian.sqrt");
+ else if (var_power == 0 & link_power == 1.0) print ("Gaussian.id");
+ else if (var_power == 0 ) print ("Gaussian.power_nonlog");
+ else if (var_power == 1.0 & link_power == -1.0) print ("Poisson.inverse");
+ else if (var_power == 1.0 & link_power == 0) print ("Poisson.log");
+ else if (var_power == 1.0 & link_power == 0.5) print ("Poisson.sqrt");
+ else if (var_power == 1.0 & link_power == 1.0) print ("Poisson.id");
+ else if (var_power == 1.0 ) print ("Poisson.power_nonlog");
+ else if (var_power == 2.0 & link_power == -1.0) print ("Gamma.inverse");
+ else if (var_power == 2.0 & link_power == 0) print ("Gamma.log");
+ else if (var_power == 2.0 & link_power == 0.5) print ("Gamma.sqrt");
+ else if (var_power == 2.0 & link_power == 1.0) print ("Gamma.id");
+ else if (var_power == 2.0 ) print ("Gamma.power_nonlog");
+ else if (var_power == 3.0 & link_power == -2.0) print ("InvGaussian.1/mu^2");
+ else if (var_power == 3.0 & link_power == -1.0) print ("InvGaussian.inverse");
+ else if (var_power == 3.0 & link_power == 0) print ("InvGaussian.log");
+ else if (var_power == 3.0 & link_power == 0.5) print ("InvGaussian.sqrt");
+ else if (var_power == 3.0 & link_power == 1.0) print ("InvGaussian.id");
+ else if (var_power == 3.0 ) print ("InvGaussian.power_nonlog");
+ else if ( link_power == 0) print ("PowerDist.log");
+ else print ("PowerDist.power_nonlog");
+ }
if (ncol_y == 1 & dist_type == 2)
{
print ("Error: Bernoulli response matrix has not been converted into two-column format.");
@@ -550,16 +550,16 @@ check_if_supported =
if (ncol_y == 2 & dist_type == 2 & link_type >= 1 & link_type <= 5)
{ # BINOMIAL/BERNOULLI DISTRIBUTION
is_supported = 1;
- if (link_type == 1 & link_power == -1.0) {print ("Binomial.inverse"); } else {
- if (link_type == 1 & link_power == 0) {print ("Binomial.log"); } else {
- if (link_type == 1 & link_power == 0.5) {print ("Binomial.sqrt"); } else {
- if (link_type == 1 & link_power == 1.0) {print ("Binomial.id"); } else {
- if (link_type == 1) {print ("Binomial.power_nonlog"); } else {
- if (link_type == 2) {print ("Binomial.logit"); } else {
- if (link_type == 3) {print ("Binomial.probit"); } else {
- if (link_type == 4) {print ("Binomial.cloglog"); } else {
- if (link_type == 5) {print ("Binomial.cauchit"); }
- } }}}}} }}}
+ if (link_type == 1 & link_power == -1.0) print ("Binomial.inverse");
+ else if (link_type == 1 & link_power == 0) print ("Binomial.log");
+ else if (link_type == 1 & link_power == 0.5) print ("Binomial.sqrt");
+ else if (link_type == 1 & link_power == 1.0) print ("Binomial.id");
+ else if (link_type == 1) print ("Binomial.power_nonlog");
+ else if (link_type == 2) print ("Binomial.logit");
+ else if (link_type == 3) print ("Binomial.probit");
+ else if (link_type == 4) print ("Binomial.cloglog");
+ else if (link_type == 5) print ("Binomial.cauchit");
+ }
if (is_supported == 0) {
print ("Response matrix with " + ncol_y + " columns, distribution family (" + dist_type + ", " + var_power
+ ") and link family (" + link_type + ", " + link_power + ") are NOT supported together.");
@@ -578,21 +578,22 @@ return (Matrix[double] beta, double saturated_log_l, int isNaN)
y_corr = Y [, 1] / (n_corr + is_n_zero) + (0.5 - Y [, 1]) * is_n_zero;
}
linear_terms = y_corr;
- if (dist_type == 1 & link_type == 1) { # POWER DISTRIBUTION
- if (link_power == 0) {
+ if (dist_type == 1 & link_type == 1)
+ { # POWER DISTRIBUTION
+ if (link_power == 0) {
if (sum (y_corr < 0) == 0) {
is_zero_y_corr = (y_corr == 0);
linear_terms = log (y_corr + is_zero_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
} else { isNaN = 1; }
- } else { if (link_power == 1.0) {
+ } else if (link_power == 1.0) {
linear_terms = y_corr;
- } else { if (link_power == -1.0) {
+ } else if (link_power == -1.0) {
linear_terms = 1.0 / y_corr;
- } else { if (link_power == 0.5) {
+ } else if (link_power == 0.5) {
if (sum (y_corr < 0) == 0) {
linear_terms = sqrt (y_corr);
} else { isNaN = 1; }
- } else { if (link_power > 0) {
+ } else if (link_power > 0) {
if (sum (y_corr < 0) == 0) {
is_zero_y_corr = (y_corr == 0);
linear_terms = (y_corr + is_zero_y_corr) ^ link_power - is_zero_y_corr;
@@ -601,21 +602,21 @@ return (Matrix[double] beta, double saturated_log_l, int isNaN)
if (sum (y_corr <= 0) == 0) {
linear_terms = y_corr ^ link_power;
} else { isNaN = 1; }
- }}}}}
+ }
}
if (dist_type == 2 & link_type >= 1 & link_type <= 5)
{ # BINOMIAL/BERNOULLI DISTRIBUTION
- if (link_type == 1 & link_power == 0) { # Binomial.log
+ if (link_type == 1 & link_power == 0) { # Binomial.log
if (sum (y_corr < 0) == 0) {
is_zero_y_corr = (y_corr == 0);
linear_terms = log (y_corr + is_zero_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
} else { isNaN = 1; }
- } else { if (link_type == 1 & link_power > 0) { # Binomial.power_nonlog pos
+ } else if (link_type == 1 & link_power > 0) { # Binomial.power_nonlog pos
if (sum (y_corr < 0) == 0) {
is_zero_y_corr = (y_corr == 0);
linear_terms = (y_corr + is_zero_y_corr) ^ link_power - is_zero_y_corr;
} else { isNaN = 1; }
- } else { if (link_type == 1) { # Binomial.power_nonlog neg
+ } else if (link_type == 1) { # Binomial.power_nonlog neg
if (sum (y_corr <= 0) == 0) {
linear_terms = y_corr ^ link_power;
} else { isNaN = 1; }
@@ -626,20 +627,20 @@ return (Matrix[double] beta, double saturated_log_l, int isNaN)
if (link_type == 2) { # Binomial.logit
linear_terms = log (y_corr / (1.0 - y_corr))
+ is_one_y_corr / (1.0 - is_one_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
- } else { if (link_type == 3) { # Binomial.probit
+ } else if (link_type == 3) { # Binomial.probit
y_below_half = y_corr + (1.0 - 2.0 * y_corr) * (y_corr > 0.5);
t = sqrt (- 2.0 * log (y_below_half));
approx_inv_Gauss_CDF = - t + (2.515517 + t * (0.802853 + t * 0.010328)) / (1.0 + t * (1.432788 + t * (0.189269 + t * 0.001308)));
linear_terms = approx_inv_Gauss_CDF * (1.0 - 2.0 * (y_corr > 0.5))
+ is_one_y_corr / (1.0 - is_one_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
- } else { if (link_type == 4) { # Binomial.cloglog
+ } else if (link_type == 4) { # Binomial.cloglog
linear_terms = log (- log (1.0 - y_corr))
- log (- log (0.5)) * (is_zero_y_corr + is_one_y_corr)
+ is_one_y_corr / (1.0 - is_one_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
- } else { if (link_type == 5) { # Binomial.cauchit
+ } else if (link_type == 5) { # Binomial.cauchit
linear_terms = tan ((y_corr - 0.5) * 3.1415926535897932384626433832795)
+ is_one_y_corr / (1.0 - is_one_y_corr) - is_zero_y_corr / (1.0 - is_zero_y_corr);
- }} }}}}}
+ } }
}
if (isNaN == 0) {
@@ -651,13 +652,13 @@ return (Matrix[double] beta, double saturated_log_l, int isNaN)
(dist_type == 2 & link_type >= 2))
{
desired_eta = 0.0;
- } else { if (link_type == 1 & link_power == 0) {
+ } else if (link_type == 1 & link_power == 0) {
desired_eta = log (0.5);
- } else { if (link_type == 1) {
+ } else if (link_type == 1) {
desired_eta = 0.5 ^ link_power;
} else {
desired_eta = 0.5;
- }}}
+ }
beta = matrix (0.0, rows = ncol(X), cols = 1);
@@ -710,14 +711,15 @@ glm_dist = function (Matrix[double] linear_terms, Matrix[double] Y,
flip_neg [1, 2] = -1;
flip_neg [2, 1] = 1;
- if (dist_type == 1 & link_type == 1) { # POWER DISTRIBUTION
+ if (dist_type == 1 & link_type == 1)
+ { # POWER DISTRIBUTION
y_mean = zeros_r;
- if (link_power == 0) {
+ if (link_power == 0) {
y_mean = exp (linear_terms);
y_mean_pow = y_mean ^ (1 - var_power);
w = y_mean_pow * y_mean;
g_Y = y_mean_pow * (Y - y_mean);
- } else { if (link_power == 1.0) {
+ } else if (link_power == 1.0) {
y_mean = linear_terms;
w = y_mean ^ (- var_power);
g_Y = w * (Y - y_mean);
@@ -727,7 +729,8 @@ glm_dist = function (Matrix[double] linear_terms, Matrix[double] Y,
c2 = (2 - var_power) / link_power - 2;
g_Y = (linear_terms ^ c1) * (Y - y_mean) / link_power;
w = (linear_terms ^ c2) / (link_power ^ 2);
- } }}
+ }
+ }
if (dist_type == 2 & link_type >= 1 & link_type <= 5)
{ # BINOMIAL/BERNOULLI DISTRIBUTION
if (link_type == 1) { # BINOMIAL.POWER LINKS
@@ -739,9 +742,9 @@ glm_dist = function (Matrix[double] linear_terms, Matrix[double] Y,
vec1 = zeros_r;
if (link_power == 0.5) {
vec1 = 1 / (1 - linear_terms ^ 2);
- } else { if (sum (linear_terms < 0) == 0) {
+ } else if (sum (linear_terms < 0) == 0) {
vec1 = linear_terms ^ (- 2 + 1 / link_power) / (1 - linear_terms ^ (1 / link_power));
- } else {isNaN = 1;}}
+ } else {isNaN = 1;}
# We want a "zero-protected" version of
# vec2 = Y [, 1] / linear_terms;
is_y_0 = (Y [, 1] == 0);
@@ -761,7 +764,7 @@ glm_dist = function (Matrix[double] linear_terms, Matrix[double] Y,
Y_prob = Y_prob * ((1.0 - rowSums (is_LT_infinite)) %*% ones_2) + is_LT_infinite;
g_Y = rowSums (Y * (Y_prob %*% flip_neg)); ### = y_residual;
w = rowSums (Y * (Y_prob %*% flip_pos) * Y_prob); ### = y_variance;
- } else { if (link_type == 3) { # Binomial.probit
+ } else if (link_type == 3) { # Binomial.probit
is_lt_pos = (linear_terms >= 0);
t_gp = 1.0 / (1.0 + abs (finite_linear_terms) * 0.231641888); # 0.231641888 = 0.3275911 / sqrt (2.0)
pt_gp = t_gp * ( 0.254829592
@@ -774,14 +777,14 @@ glm_dist = function (Matrix[double] linear_terms, Matrix[double] Y,
vec2 = Y [, 1] - rowSums (Y) * is_lt_pos + the_gauss_exp * pt_gp * rowSums (Y) * (is_lt_pos - 0.5);
w = the_gauss_exp * (one_over_sqrt_two_pi ^ 2) * rowSums (Y) / vec1;
g_Y = one_over_sqrt_two_pi * vec2 / vec1;
- } else { if (link_type == 4) { # Binomial.cloglog
+ } else if (link_type == 4) { # Binomial.cloglog
the_exp = exp (linear_terms)
the_exp_exp = exp (- the_exp);
is_too_small = ((10000000 + the_exp) == 10000000);
the_exp_ratio = (1 - is_too_small) * (1 - the_exp_exp) / (the_exp + is_too_small) + is_too_small * (1 - the_exp / 2);
g_Y = (rowSums (Y) * the_exp_exp - Y [, 2]) / the_exp_ratio;
w = the_exp_exp * the_exp * rowSums (Y) / the_exp_ratio;
- } else { if (link_type == 5) { # Binomial.cauchit
+ } else if (link_type == 5) { # Binomial.cauchit
Y_prob = 0.5 + (atan (finite_linear_terms) %*% p_one_m_one) / 3.1415926535897932384626433832795;
Y_prob = Y_prob * ((1.0 - rowSums (is_LT_infinite)) %*% ones_2) + is_LT_infinite;
y_residual = Y [, 1] * Y_prob [, 2] - Y [, 2] * Y_prob [, 1];
@@ -789,7 +792,7 @@ glm_dist = function (Matrix[double] linear_terms, Matrix[double] Y,
link_gradient_normalized = (1 + linear_terms ^ 2) * 3.1415926535897932384626433832795;
g_Y = rowSums (Y) * y_residual / (var_function * link_gradient_normalized);
w = (rowSums (Y) ^ 2) / (var_function * link_gradient_normalized ^ 2);
- }}}}
+ }
}
}
}
@@ -813,80 +816,81 @@ glm_log_likelihood_part = function (Matrix[double] linear_terms, Matrix[double]
b_cumulant = exp (linear_terms);
is_natural_parameter_log_zero = (linear_terms == -1.0/0.0);
natural_parameters = replace (target = linear_terms, pattern = -1.0/0.0, replacement = 0);
- } else { if (var_power == 1.0 & link_power == 1.0) { # Poisson.id
+ } else if (var_power == 1.0 & link_power == 1.0) { # Poisson.id
if (sum (linear_terms < 0) == 0) {
b_cumulant = linear_terms;
is_natural_parameter_log_zero = (linear_terms == 0);
natural_parameters = log (linear_terms + is_natural_parameter_log_zero);
} else {isNaN = 1;}
- } else { if (var_power == 1.0 & link_power == 0.5) { # Poisson.sqrt
+ } else if (var_power == 1.0 & link_power == 0.5) { # Poisson.sqrt
if (sum (linear_terms < 0) == 0) {
b_cumulant = linear_terms ^ 2;
is_natural_parameter_log_zero = (linear_terms == 0);
natural_parameters = 2.0 * log (linear_terms + is_natural_parameter_log_zero);
} else {isNaN = 1;}
- } else { if (var_power == 1.0 & link_power > 0) { # Poisson.power_nonlog, pos
+ } else if (var_power == 1.0 & link_power > 0) { # Poisson.power_nonlog, pos
if (sum (linear_terms < 0) == 0) {
is_natural_parameter_log_zero = (linear_terms == 0);
b_cumulant = (linear_terms + is_natural_parameter_log_zero) ^ (1.0 / link_power) - is_natural_parameter_log_zero;
natural_parameters = log (linear_terms + is_natural_parameter_log_zero) / link_power;
} else {isNaN = 1;}
- } else { if (var_power == 1.0) { # Poisson.power_nonlog, neg
+ } else if (var_power == 1.0) { # Poisson.power_nonlog, neg
if (sum (linear_terms <= 0) == 0) {
b_cumulant = linear_terms ^ (1.0 / link_power);
natural_parameters = log (linear_terms) / link_power;
} else {isNaN = 1;}
- } else { if (var_power == 2.0 & link_power == -1.0) { # Gamma.inverse
+ } else if (var_power == 2.0 & link_power == -1.0) { # Gamma.inverse
if (sum (linear_terms <= 0) == 0) {
b_cumulant = - log (linear_terms);
natural_parameters = - linear_terms;
} else {isNaN = 1;}
- } else { if (var_power == 2.0 & link_power == 1.0) { # Gamma.id
+ } else if (var_power == 2.0 & link_power == 1.0) { # Gamma.id
if (sum (linear_terms <= 0) == 0) {
b_cumulant = log (linear_terms);
natural_parameters = - 1.0 / linear_terms;
} else {isNaN = 1;}
- } else { if (var_power == 2.0 & link_power == 0) { # Gamma.log
+ } else if (var_power == 2.0 & link_power == 0) { # Gamma.log
b_cumulant = linear_terms;
natural_parameters = - exp (- linear_terms);
- } else { if (var_power == 2.0) { # Gamma.power_nonlog
+ } else if (var_power == 2.0) { # Gamma.power_nonlog
if (sum (linear_terms <= 0) == 0) {
b_cumulant = log (linear_terms) / link_power;
natural_parameters = - linear_terms ^ (- 1.0 / link_power);
} else {isNaN = 1;}
- } else { if (link_power == 0) { # PowerDist.log
+ } else if (link_power == 0) { # PowerDist.log
natural_parameters = exp (linear_terms * (1.0 - var_power)) / (1.0 - var_power);
b_cumulant = exp (linear_terms * (2.0 - var_power)) / (2.0 - var_power);
} else { # PowerDist.power_nonlog
if (-2 * link_power == 1.0 - var_power) {
natural_parameters = 1.0 / (linear_terms ^ 2) / (1.0 - var_power);
- } else { if (-1 * link_power == 1.0 - var_power) {
+ } else if (-1 * link_power == 1.0 - var_power) {
natural_parameters = 1.0 / linear_terms / (1.0 - var_power);
- } else { if ( link_power == 1.0 - var_power) {
+ } else if ( link_power == 1.0 - var_power) {
natural_parameters = linear_terms / (1.0 - var_power);
- } else { if ( 2 * link_power == 1.0 - var_power) {
+ } else if ( 2 * link_power == 1.0 - var_power) {
natural_parameters = linear_terms ^ 2 / (1.0 - var_power);
+ } else if (sum (linear_terms <= 0) == 0) {
+ power = (1.0 - var_power) / link_power;
+ natural_parameters = (linear_terms ^ power) / (1.0 - var_power);
} else {
- if (sum (linear_terms <= 0) == 0) {
- power = (1.0 - var_power) / link_power;
- natural_parameters = (linear_terms ^ power) / (1.0 - var_power);
- } else {isNaN = 1;}
- }}}}
+ isNaN = 1;
+ }
+
if (-2 * link_power == 2.0 - var_power) {
b_cumulant = 1.0 / (linear_terms ^ 2) / (2.0 - var_power);
- } else { if (-1 * link_power == 2.0 - var_power) {
+ } else if (-1 * link_power == 2.0 - var_power) {
b_cumulant = 1.0 / linear_terms / (2.0 - var_power);
- } else { if ( link_power == 2.0 - var_power) {
+ } else if ( link_power == 2.0 - var_power) {
b_cumulant = linear_terms / (2.0 - var_power);
- } else { if ( 2 * link_power == 2.0 - var_power) {
+ } else if ( 2 * link_power == 2.0 - var_power) {
b_cumulant = linear_terms ^ 2 / (2.0 - var_power);
+ } else if (sum (linear_terms <= 0) == 0) {
+ power = (2.0 - var_power) / link_power;
+ b_cumulant = (linear_terms ^ power) / (2.0 - var_power);
} else {
- if (sum (linear_terms <= 0) == 0) {
- power = (2.0 - var_power) / link_power;
- b_cumulant = (linear_terms ^ power) / (2.0 - var_power);
- } else {isNaN = 1;}
- }}}}
- }}}}} }}}}}
+ isNaN = 1;
+ }
+ }
if (sum (is_natural_parameter_log_zero * abs (Y)) > 0) {
log_l = -1.0 / 0.0;
isNaN = 1;
@@ -951,13 +955,13 @@ binomial_probability_two_column =
if (link_type == 1) { # Binomial.power
if (link_power == 0) { # Binomial.log
Y_prob = exp (linear_terms) %*% p_one_m_one + ones_r %*% zero_one;
- } else { if (link_power == 0.5) { # Binomial.sqrt
+ } else if (link_power == 0.5) { # Binomial.sqrt
Y_prob = (linear_terms ^ 2) %*% p_one_m_one + ones_r %*% zero_one;
- } else { # Binomial.power_nonlog
- if (sum (linear_terms < 0) == 0) {
- Y_prob = (linear_terms ^ (1.0 / link_power)) %*% p_one_m_one + ones_r %*% zero_one;
- } else {isNaN = 1;}
- }}
+ } else if (sum (linear_terms < 0) == 0) { # Binomial.power_nonlog
+ Y_prob = (linear_terms ^ (1.0 / link_power)) %*% p_one_m_one + ones_r %*% zero_one;
+ } else {
+ isNaN = 1;
+ }
} else { # Binomial.non_power
is_LT_pos_infinite = (linear_terms == 1.0/0.0);
is_LT_neg_infinite = (linear_terms == -1.0/0.0);
@@ -967,7 +971,7 @@ binomial_probability_two_column =
if (link_type == 2) { # Binomial.logit
Y_prob = exp (finite_linear_terms) %*% one_zero + ones_r %*% zero_one;
Y_prob = Y_prob / (rowSums (Y_prob) %*% ones_2);
- } else { if (link_type == 3) { # Binomial.probit
+ } else if (link_type == 3) { # Binomial.probit
lt_pos_neg = (finite_linear_terms >= 0) %*% p_one_m_one + ones_r %*% zero_one;
t_gp = 1.0 / (1.0 + abs (finite_linear_terms) * 0.231641888); # 0.231641888 = 0.3275911 / sqrt (2.0)
pt_gp = t_gp * ( 0.254829592
@@ -977,17 +981,17 @@ binomial_probability_two_column =
+ t_gp * 1.061405429))));
the_gauss_exp = exp (- (finite_linear_terms ^ 2) / 2.0);
Y_prob = lt_pos_neg + ((the_gauss_exp * pt_gp) %*% ones_2) * (0.5 - lt_pos_neg);
- } else { if (link_type == 4) { # Binomial.cloglog
+ } else if (link_type == 4) { # Binomial.cloglog
the_exp = exp (finite_linear_terms);
the_exp_exp = exp (- the_exp);
is_too_small = ((10000000 + the_exp) == 10000000);
Y_prob [, 1] = (1 - is_too_small) * (1 - the_exp_exp) + is_too_small * the_exp * (1 - the_exp / 2);
Y_prob [, 2] = the_exp_exp;
- } else { if (link_type == 5) { # Binomial.cauchit
+ } else if (link_type == 5) { # Binomial.cauchit
Y_prob = 0.5 + (atan (finite_linear_terms) %*% p_one_m_one) / 3.1415926535897932384626433832795;
} else {
isNaN = 1;
- }}}}
+ }
Y_prob = Y_prob * ((1.0 - rowSums (is_LT_infinite)) %*% ones_2) + is_LT_infinite;
} }