Add files via upload

faosorios · Jan 22, 2024 · 2393c30 · 2393c30
1 parent b9bccad
commit 2393c30
Show file tree

Hide file tree

Showing 3 changed files with 49 additions and 6 deletions.
diff --git a/pkg/src/R_init_fastmatrix.c b/pkg/src/R_init_fastmatrix.c
@@ -1,4 +1,4 @@
-/* $ID: R_init_fastmatrix.c, last updated 2023-12-06, F.Osorio */
+/* $ID: R_init_fastmatrix.c, last updated 2024-01-03, F.Osorio */
 
 #include "fastmatrix.h"
 #include <R_ext/Rdynload.h>
@@ -62,6 +62,7 @@ static const R_CMethodDef CEntries[]  = {
   CALLDEF(symmetrizer_prod,       6),
   CALLDEF(urzua_ALM,              5),
   CALLDEF(wilson_hilferty_chisq,  4),
+  CALLDEF(wilson_hilferty_gamma,  5),
   CALLDEF(whitening_chol,         4),
   {NULL, NULL, 0}
 };
@@ -187,8 +188,11 @@ void R_init_fastmatrix(DllInfo *info) {
   /* distances */
   FM_REGDEF(FM_pythag);
   FM_REGDEF(FM_mahalanobis);
+  /* Wilson-Hilferty transformation */
   FM_REGDEF(FM_WH_chisq);
+  FM_REGDEF(FM_WH_gamma);
   FM_REGDEF(FM_WH_F);
+  FM_REGDEF(FM_WH_Laplace);
   /* products */
   FM_REGDEF(FM_compensated_product);
   FM_REGDEF(FM_two_product_FMA);

diff --git a/pkg/src/fastmatrix.h b/pkg/src/fastmatrix.h
@@ -1,4 +1,4 @@
-/* ID: fastmatrix.h, last updated 2023-12-06, F.Osorio */
+/* ID: fastmatrix.h, last updated 2024-01-03, F.Osorio */
 
 #ifndef FASTMATRIX_H
 #define FASTMATRIX_H
@@ -122,7 +122,10 @@ void F77_NAME(median_center)(double *, int *, int *, int *, double *, int *, int
 void geometric_mean(double *, int *, double *);
 void mahal_distances(double *, int *, int *, double *, double *, int *, double *);
 void skewness_and_kurtosis(double *, int *, int *, double *, double *, double *, int *);
+
+/* Wilson-Hilferty transformation */
 void wilson_hilferty_chisq(double *, int *, int *, double *);
+void wilson_hilferty_gamma(double *, int *, double *, double *, double *);
 
 /* tests for normality based on the standardized third and fourth moments */
 void doornik_hansen(double *, int *, double *, double *, double *);
@@ -272,8 +275,12 @@ void FM_gls_GQR(double *, int, int, int, double *, double *, double *, int *);
 double FM_pythag(double, double);
 double FM_mahalanobis(double *, int, double *, double *);
 void FM_mahal_distances(double *, int, int, double *, double *, int, double *);
+
+/* Wilson-Hilferty transformation */
 void FM_WH_chisq(double *, int, int, double *);
+void FM_WH_gamma(double *, int, double, double, double *);
 void FM_WH_F(double *, int, int, double, double *);
+void FM_WH_Laplace(double *, int, int, double *);
 
 /* products */
 void FM_two_product_FMA(double, double, double *, double *);

diff --git a/pkg/src/wilson_hilferty.c b/pkg/src/wilson_hilferty.c
@@ -1,4 +1,4 @@
-/* $ID: wilson_hilferty.c, last updated 2020-11-24, F.Osorio */
+/* $ID: wilson_hilferty.c, last updated 2024-01-03, F.Osorio */
 
 #include "fastmatrix.h"
 
@@ -10,10 +10,29 @@ wilson_hilferty_chisq(double *distances, int *n, int *p, double *z)
   FM_WH_chisq(distances, *n, *p, z);
 }
 
+void
+wilson_hilferty_gamma(double *y, int *n, double *shape, double *scale, double *z)
+{ /* Wilson-Hilferty transformation for gamma variables (to be called by R) */
+  FM_WH_gamma(y, *n, *shape, *scale, z);
+}
+
+void
+FM_WH_gamma(double *y, int n, double shape, double scale, double *z)
+{ /* Wilson-Hilferty transformation for gamma variables */
+  double t, mean, sd, q = 1./3.;
+
+  for (int i = 0; i < n; i++) {
+    t = *y++ / (shape / scale);
+    mean = 1. - 1. / (9. * shape);
+    sd = 1. / sqrt(9. * shape);
+    *z++ = (R_pow(t, q) - mean) / sd;
+  }
+}
+
 void
 FM_WH_chisq(double *distances, int n, int p, double *z)
 { /* Wilson-Hilferty transformation for chi-squared variables */
-  double f, q = 1. / 3., s = 2. / 9.;
+  double f, q = 1./3., s = 2./9.;
 
   for (int i = 0; i < n; i++) {
     f = *distances++ / p;
@@ -24,11 +43,24 @@ FM_WH_chisq(double *distances, int n, int p, double *z)
 void
 FM_WH_F(double *distances, int n, int p, double eta, double *z)
 { /* Wilson-Hilferty transformation for F variables */
-  double f, q = 1. / 3., r = 2. / 3., s = 2. / 9.;
+  double f, q = 1./3., r = 2./3., s = 2./9.;
 
   for (int i = 0; i < n; i++) {
-    f = *distances++ / p;
+    f  = *distances++ / p;
     f /= 1. - 2. * eta;
     *z++ = ((1. - s * eta) * R_pow(f, q) - (1. - s / p)) / sqrt(s * eta * R_pow(f, r) + s / p);
   }
 }
+
+void
+FM_WH_Laplace(double *distances, int n, int p, double *z)
+{ /* Wilson-Hilferty transformation for Laplace variables */
+  double f, mean, sd, q = 1./3.;
+
+  for (int i = 0; i < n; i++) {
+    f = sqrt(*distances++) / (2. * p);
+    mean = 1. - 1. / (9. * p);
+    sd = 1. / sqrt(9. * p);
+    *z++ = (R_pow(f, q) - mean) / sd;
+  }
+}