Some statistics function impl in R

Co-authored-by: Copilot <copilot@github.com>

Author: XuhuaHuang

R/README.md

@@ -255,3 +255,134 @@ If the end of the function body is not an explicit `return` statement, the value
 - `warning("warning message")` to throw a warning and continue execution
 - `try(expr)` to execute an expression and catch any errors without stopping execution
 - `tryCatch(expr, error = function(e) { ... })` to execute an expression and handle errors with a custom function
+
+## Wichmann-Hill Pseudo-Random Number Generator
+
+```r
+wh = function(seed, start = 0, end = 1) {
+    seed[1] = (171 * seed[1]) %% 30269
+    seed[2] = (172 * seed[2]) %% 30307
+    seed[3] = (170 * seed[3]) %% 30323
+    x = (seed[1]/30269 + seed[2]/30307 + seed[3]/30323) %% 1
+    start + (end - start) * x
+}
+whv = function(seed, start = 0, end = 1) {
+    y = c(171, 172, 170)
+    z = c(30269, 30307, 30323)
+    seed = (y * seed) %% z
+    x = sum(seed / z) %% 1
+    start + (end - start) * x
+}
+```
+
+## Programming Style
+
+- Use meaningful variable and function names
+- Modularize code into functions
+- Document code with comments and function documentation
+- Use consistent indentation and spacing
+- Use existing functions and libraries when possible
+- Use parentheses to make grouping and operator precedence clear
+- Avoid unnecessary loops and use vectorized operations when possible
+- Avoid using recursive functions when iterative solutions are more efficient
+
+## IQR
+
+```r
+IQR(x, na.rm = FALSE)
+q <- quantile(x, probs = c(0.25, 0.75))
+iqr <- q[2] - q[1]
+iqr
+```
+
+## Monte Carlo Demonstration of `CLT` in `R`
+
+```r
+set.seed(42)
+
+# Population: Exponential distribution (skewed)
+n <- 30              # sample size
+num_sim <- 10000     # number of Monte Carlo repetitions
+
+# True parameters
+mu <- 1              # mean of Exp(1)
+sigma <- 1           # std dev of Exp(1)
+
+z_values <- replicate(num_sim, {
+  x <- rexp(n, rate = 1)         # draw sample
+  x_bar <- mean(x)               # sample mean
+  sqrt(n) * (x_bar - mu)         # standardized value
+})
+
+hist(z_values, probability = TRUE, breaks = 50,
+     main = "CLT Demonstration (Monte Carlo)",
+     xlab = expression(sqrt(n) * (bar(X) - mu)))
+
+# Overlay theoretical normal curve
+curve(dnorm(x, mean = 0, sd = sigma), 
+      col = "red", lwd = 2, add = TRUE)
+```
+
+## Alpha-trimmed Mean
+
+```r
+alpha_trimmed_mean_safe <- function(x, alpha) {
+  if (!is.numeric(x)) {
+    stop("x must be numeric")
+  }
+
+  if (alpha < 0 || alpha >= 0.5) {
+    stop("alpha must be in [0, 0.5)")
+  }
+
+  n <- length(x)
+  k <- floor(alpha * n)
+
+  x_sorted <- sort(x)
+
+  if (2 * k >= n) {
+    stop("Trimming removes all data")
+  }
+
+  mean(x_sorted[(k + 1):(n - k)])
+}
+```
+
+## Sample Correlation
+
+```r
+sample_correlation <- function(x, y) {
+  if (length(x) != length(y)) {
+    stop("Vectors must have the same length")
+  }
+
+  x_bar <- mean(x)
+  y_bar <- mean(y)
+
+  numerator <- sum((x - x_bar) * (y - y_bar))
+  denominator <- sqrt(sum((x - x_bar)^2)) * sqrt(sum((y - y_bar)^2))
+
+  numerator / denominator
+}
+
+sample_correlation_stable <- function(x, y) {
+  if (length(x) != length(y)) {
+    stop("Vectors must have the same length")
+  }
+
+  n <- length(x)
+
+  cov_xy <- sum(x * y) / n - mean(x) * mean(y)
+  sd_x <- sqrt(sum(x^2) / n - mean(x)^2)
+  sd_y <- sqrt(sum(y^2) / n - mean(y)^2)
+
+  cov_xy / (sd_x * sd_y)
+}
+```
+
+## Debugging and Maintenance
+
+- Error messages generated by `R` functions can be helpful for debugging.
+- Use `browser()` function to set breakpoints.
+- Use `print()` or `cat()` to display intermediate values for debugging.
+- Do not use `print()` in return statements of functions, just use `return(return_object)`.