stats/src/bandwidths.c

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  bandwidth.c by W. N. Venables and B. D. Ripley  Copyright (C) 1994-2001
 *  Copyright (C) 2012-2020  The R Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.R-project.org/Licenses/
 */

#include <stdlib.h> //abs
#include <math.h>
#include <Rmath.h> // M_* constants
#include <Rinternals.h>

// or include "stats.h"
#ifdef ENABLE_NLS
#include <libintl.h>
#define _(String) dgettext ("stats", String)
#else
#define _(String) (String)
#endif

#define DELMAX 1000
/* Avoid slow and possibly error-producing underflows by cutting off at
   plus/minus sqrt(DELMAX) std deviations */
/* Formulae (6.67) and (6.69) of Scott (1992), the latter corrected. */

SEXP bw_ucv(SEXP sn, SEXP sd, SEXP cnt, SEXP sh)
{
    double h = asReal(sh), d = asReal(sd), sum = 0.0, term, u;
    int n = asInteger(sn), nbin = LENGTH(cnt);
    double *x = REAL(cnt);
    for (int i = 0; i < nbin; i++) {
	double delta = i * d / h;
	delta *= delta;
	if (delta >= DELMAX) break;
	term = exp(-delta / 4.) - sqrt(8.0) * exp(-delta / 2.);
	sum += term * x[i];
    }
    u = (0.5 + sum/n) / (n * h * M_SQRT_PI);
    // = 1 / (2 * n * h * sqrt(PI)) + sum / (n * n * h * sqrt(PI));
    return ScalarReal(u);
}

SEXP bw_bcv(SEXP sn, SEXP sd, SEXP cnt, SEXP sh)
{
    double h = asReal(sh), d = asReal(sd), sum = 0.0, term, u;
    int n = asInteger(sn), nbin = LENGTH(cnt);
    double *x = REAL(cnt);

    sum = 0.0;
    for (int i = 0; i < nbin; i++) {
	double delta = i * d / h; delta *= delta;
	if (delta >= DELMAX) break;
	term = exp(-delta / 4) * (delta * delta - 12 * delta + 12);
	sum += term * x[i];
    }
    u = (1 + sum/(32.0*n)) / (2.0 * n * h * M_SQRT_PI);
    // = 1 / (2 * n * h * sqrt(PI)) + sum / (64 * n * n * h * sqrt(PI));
    return ScalarReal(u);
}

SEXP bw_phi4(SEXP sn, SEXP sd, SEXP cnt, SEXP sh)
{
    double h = asReal(sh), d = asReal(sd), sum = 0.0, term, u;
    int n = asInteger(sn), nbin = LENGTH(cnt);
    double *x = REAL(cnt);

    for (int i = 0; i < nbin; i++) {
	double delta = i * d / h; delta *= delta;
	if (delta >= DELMAX) break;
	term = exp(-delta / 2.) * (delta * delta - 6. * delta + 3.);
	sum += term * x[i];
    }
    sum = 2. * sum + n * 3.;	/* add in diagonal */
    u = sum / ((double)n * (n - 1) * pow(h, 5.0)) * M_1_SQRT_2PI;
    // = sum / (n * (n - 1) * pow(h, 5.0) * sqrt(2 * PI));
    return ScalarReal(u);
}

SEXP bw_phi6(SEXP sn, SEXP sd, SEXP cnt, SEXP sh)
{
    double h = asReal(sh), d = asReal(sd), sum = 0.0, term, u;
    int n = asInteger(sn), nbin = LENGTH(cnt);
    double *x = REAL(cnt);

    for (int i = 0; i < nbin; i++) {
	double delta = i * d / h; delta *= delta;
	if (delta >= DELMAX) break;
	term = exp(-delta / 2) *
	    (delta * delta * delta - 15 * delta * delta + 45 * delta - 15);
	sum += term * x[i];
    }
    sum = 2. * sum - 15. * n;	/* add in diagonal */
    u = sum / ((double)n * (n - 1) * pow(h, 7.0)) * M_1_SQRT_2PI;
    // = sum / (n * (n - 1) * pow(h, 7.0) * sqrt(2 * PI));
    return ScalarReal(u);
}

/*
   Use double cnt as from R 3.4.0, as counts can exceed INT_MAX for
   large n (65537 in the worse case but typically not at n = 1 million
   for a smooth distribution -- and this is by default no longer used
   for n > 500).
*/

SEXP bw_den(SEXP nbin, SEXP sx)
{
    int nb = asInteger(nbin), n = LENGTH(sx);
    double xmin, xmax, rang, dd, *x = REAL(sx);

    xmin = R_PosInf; xmax = R_NegInf;
    for (int i = 0; i < n; i++) {
	if(!R_FINITE(x[i]))
	    error(_("non-finite x[%d] in bandwidth calculation"), i+1);
	if(x[i] < xmin) xmin = x[i];
	if(x[i] > xmax) xmax = x[i];
    }
    rang = (xmax - xmin) * 1.01;
    dd = rang / nb;

    SEXP ans = PROTECT(allocVector(VECSXP, 2)),
	sc = SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, nb));
    SET_VECTOR_ELT(ans, 0, ScalarReal(dd));
    double *cnt = REAL(sc);
    for (int i = 0; i < nb; i++) cnt[i] = 0.0;

    for (int i = 1; i < n; i++) {
	int ii = (int)(x[i] / dd);
	for (int j = 0; j < i; j++) {
	    int jj = (int)(x[j] / dd);
	    cnt[abs(ii - jj)] += 1.0;
	}
    }

    UNPROTECT(1);
    return ans;
}

/* Input: counts for nb bins */
SEXP bw_den_binned(SEXP sx)
{
    int nb = LENGTH(sx);
    int *x = INTEGER(sx);

    SEXP ans = PROTECT(allocVector(REALSXP, nb));
    double *cnt = REAL(ans);
    for (int ib = 0; ib < nb; ib++) cnt[ib] = 0.0;

    for (int ii = 0; ii < nb; ii++) {
	double w = x[ii]; // avoid int overflows below
	cnt[0] += w*(w-1.); // don't count distances to self
	for (int jj = 0; jj < ii; jj++)
	    cnt[ii - jj] += w * x[jj];
    }
    cnt[0] *= 0.5; // counts in the same bin got double-counted

    UNPROTECT(1);
    return ans;
}