docs: add method documentation

Signed-off-by: Athan <kgryte@gmail.com>

Author: kgryte

lib/node_modules/@stdlib/stats/strided/dpcorrwd/src/main.c

@@ -40,6 +40,109 @@ double API_SUFFIX(stdlib_strided_dpcorrwd)( const CBLAS_INT N, const double *X,
 /**
 * Computes the sample Pearson product-moment correlation coefficient of two double-precision floating-point strided arrays using Welford's algorithm and alternative indexing semantics.
 *
+* ## Method
+*
+* -   We begin by defining the co-moment \\(C_n\\)
+*
+*     ```tex
+*     C_n = \sum_{i=1}^{n} ( x_i - \bar{x}_n ) ( y_i - \bar{y}_n )
+*     ```
+*
+*     where \\(\bar{x}_n\\) and \\(\bar{y}_n\\) are the sample means for \\(x\\) and \\(y\\), respectively.
+*
+* -   Based on Welford's method, we know the update formulas for the sample means are given by
+*
+*     ```tex
+*     \bar{x}_n = \bar{x}_{n-1} + \frac{x_n - \bar{x}_{n-1}}{n}
+*     ```
+*
+*     and
+*
+*     ```tex
+*     \bar{y}_n = \bar{y}_{n-1} + \frac{y_n - \bar{y}_{n-1}}{n}
+*     ```
+*
+* -   Substituting into the equation for \\(C_n\\) and rearranging terms
+*
+*     ```tex
+*     C_n = C_{n-1} + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})
+*     ```
+*
+*     where the apparent asymmetry arises from
+*
+*     ```tex
+*     x_n - \bar{x}_n = \frac{n-1}{n} (x_n - \bar{x}_{n-1})
+*     ```
+*
+*     and, hence, the update term can be equivalently expressed
+*
+*     ```tex
+*     \frac{n-1}{n} (x_n - \bar{x}_{n-1}) (y_n - \bar{y}_{n-1})
+*     ```
+*
+* -   The covariance can be defined
+*
+*     ```tex
+*     \begin{align*}
+*     \operatorname{cov}_n(x,y) &= \frac{C_n}{n} \\
+*     &= \frac{C_{n-1} + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n} \\
+*     &= \frac{(n-1)\operatorname{cov}_{n-1}(x,y) + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n}
+*     \end{align*}
+*     ```
+*
+* -   Applying Bessel's correction, we arrive at an update formula for calculating an unbiased sample covariance
+*
+*     ```tex
+*     \begin{align*}
+*     \operatorname{cov}_n(x,y) &= \frac{n}{n-1}\cdot\frac{(n-1)\operatorname{cov}_{n-1}(x,y) + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n} \\
+*     &= \operatorname{cov}_{n-1}(x,y) + \frac{(x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n-1} \\
+*     &= \frac{C_{n-1} + (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})}{n-1}
+*     &= \frac{C_{n-1} + (x_n - \bar{x}_{n-1}) (y_n - \bar{y}_n)}{n-1}
+*     \end{align*}
+*     ```
+*
+* -   To calculate the corrected sample standard deviation, we can use Welford's method, which can be derived as follows. We can express the variance as
+*
+*     ```tex
+*     \begin{align*}
+*     S_n &= n \sigma_n^2 \\
+*         &= \sum_{i=1}^{n} (x_i - \mu_n)^2 \\
+*         &= \biggl(\sum_{i=1}^{n} x_i^2 \biggr) - n\mu_n^2
+*     \end{align*}
+*     ```
+*
+*     Accordingly,
+*
+*     ```tex
+*     \begin{align*}
+*     S_n - S_{n-1} &= \sum_{i=1}^{n} x_i^2 - n\mu_n^2 - \sum_{i=1}^{n-1} x_i^2 + (n-1)\mu_{n-1}^2 \\
+*                   &= x_n^2 - n\mu_n^2 + (n-1)\mu_{n-1}^2 \\
+*                   &= x_n^2 - \mu_{n-1}^2 + n(\mu_{n-1}^2 - \mu_n^2) \\
+*                   &= x_n^2 - \mu_{n-1}^2 + n(\mu_{n-1} - \mu_n)(\mu_{n-1} + \mu_n) \\
+*                   &= x_n^2 - \mu_{n-1}^2 + (\mu_{n-1} - x_n)(\mu_{n-1} + \mu_n) \\
+*                   &= x_n^2 - \mu_{n-1}^2 + \mu_{n-1}^2 - x_n\mu_n - x_n\mu_{n-1} + \mu_n\mu_{n-1} \\
+*                   &= x_n^2 - x_n\mu_n - x_n\mu_{n-1} + \mu_n\mu_{n-1} \\
+*                   &= (x_n - \mu_{n-1})(x_n - \mu_n) \\
+*                   &= S_{n-1} + (x_n - \mu_{n-1})(x_n - \mu_n)
+*     \end{align*}
+*     ```
+*
+*     where we use the identity
+*
+*     ```tex
+*     x_n - \mu_{n-1} = n (\mu_n - \mu_{n-1})
+*     ```
+*
+* -   To compute the corrected sample standard deviation, we apply Bessel's correction and take the square root.
+*
+* -   The sample Pearson product-moment correlation coefficient can thus be calculated as
+*
+*     ```tex
+*     r = \frac{\operatorname{cov}_n(x,y)}{\sigma_x \sigma_y}
+*     ```
+*
+*     where \\(\sigma_x\\) and \\(\sigma_y\\) are the corrected sample standard deviations for \\(x\\) and \\(y\\), respectively.
+*
 * @param N            number of indexed elements
 * @param X            first input array
 * @param strideX      stride length of `X`