---
title: "Robustness"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Robustness}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- markdownlint-disable MD024 -->

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

## How Robustness Works

Standard LOWESS can be biased by outliers. Robustness iterations downweight
points with large residuals:

1. Fit initial LOWESS
2. Compute residuals
3. Assign robustness weights (large residuals → low weight)
4. Refit using combined distance × robustness weights
5. Repeat steps 2–4

![Robustness Methods](../man/figures/robust_method_comparison.svg)

![Robustness Iterations](../man/figures/robust_iter_comparison.svg)

---

## Robustness Methods

### Bisquare (Default)

Smooth downweighting. Points transition gradually from full weight to zero.

$$w(u) = \begin{cases} (1 - u^2)^2 & |u| < 1 \\ 0 & |u| \geq 1 \end{cases}$$

**Use when**: General purpose, balanced approach.

```{r robustness_1}
library(rfastlowess)
set.seed(42)
x <- seq(0, 2 * pi, length.out = 100)
y <- sin(x) + rnorm(100, sd = 0.3)

model <- Lowess(iterations = 3, robustness_method = "bisquare")
result <- fit(model, x, y)
cat("Smoothed y[0]:", result$y[1], "\n")
```

---

### Huber

Linear penalty beyond threshold. Less aggressive than Bisquare.

$$w(u) = \begin{cases} 1 & |u| \leq k \\ k/|u| & |u| > k \end{cases}$$

**Use when**: Moderate outliers, want to retain some influence.

```{r robustness_2}
model <- Lowess(iterations = 3, robustness_method = "huber")
result <- fit(model, x, y)
cat("Smoothed y[0]:", result$y[1], "\n")
```

---

### Talwar

Hard threshold. Points are either fully weighted or completely excluded.

$$w(u) = \begin{cases} 1 & |u| \leq k \\ 0 & |u| > k \end{cases}$$

**Use when**: Extreme outliers, want binary exclusion.

```{r robustness_3}
model <- Lowess(iterations = 3, robustness_method = "talwar")
result <- fit(model, x, y)
cat("Smoothed y[0]:", result$y[1], "\n")
```

---

## Comparison

| Method | Transition | Aggressiveness | Use Case |
| --- | --- | --- | --- |
| **Bisquare** | Smooth | Moderate | General purpose |
| **Huber** | Gradual | Mild | Preserve influence |
| **Talwar** | Hard | Strong | Extreme contamination |

---

## Detecting Outliers

Use robustness weights to identify potential outliers:

```{r robustness_5}
library(rfastlowess)
set.seed(42)
x <- seq(0, 2 * pi, length.out = 100)
y <- sin(x) + rnorm(100, sd = 0.3)
y[c(20, 50, 80)] <- y[c(20, 50, 80)] + 5  # inject outliers

model <- Lowess(iterations = 5, outputs = "weights")
result <- fit(model, x, y)

shown <- 0
for (i in seq_along(result$robustness_weights)) {
    if (result$robustness_weights[i] < 0.5 && shown < 5) {
        cat(sprintf("Potential outlier at index %d: weight = %.3f\n",
                    i, result$robustness_weights[i]))
        shown <- shown + 1
    }
}
```

---

## Scale Estimation

Residuals are scaled before computing robustness weights. Two methods:

| Method | Formula | Robustness |
| --- | --- | --- |
| **MAD** | `median(\|r - median(r)\|)` | Very robust (default) |
| **MAR** | `median(\|r\|)` | Robust, uncentered |
| **Mean** | `mean(\|r\|)` | Less robust, fastest |

![Scaling Methods Comparison](../man/figures/scaling_comparison.svg)

```{r robustness_6}
library(rfastlowess)
set.seed(42)
x <- seq(0, 2 * pi, length.out = 100)
y <- sin(x) + rnorm(100, sd = 0.3)

model <- Lowess(iterations = 3, scaling_method = "mad")
result <- fit(model, x, y)
cat("Smoothed y[0]:", result$y[1], "\n")
```

---

## Auto-Convergence

Stop iterations early when weights stabilize:

> **Performance:** Auto-convergence can significantly reduce computation
> when weights stabilize before reaching max iterations.

```{r robustness_7}
library(rfastlowess)
set.seed(42)
x <- seq(0, 2 * pi, length.out = 100)
y <- sin(x) + rnorm(100, sd = 0.3)

model <- Lowess(iterations = 10, auto_converge = 1e-6)
result <- fit(model, x, y)
cat("Smoothed y[0]:", result$y[1], "\n")
```

```{r session-info}
sessionInfo()
```