Transforming Shape

"""The transform that makes the histogram symmetric is not the transform you want.

One skewed input column, three different truths about how the target depends on
it, five candidate transforms. `skew_after` is the usual criterion — how close
to symmetric the column looks once transformed — and the three columns beside
it are what a straight line through the transformed column actually achieves.

Read `skew_after` against the three, and the usual criterion falls apart:
yeo_johnson and quantile_normal symmetrize best and are the *worst* option when
the target is linear in x. The transform to want is the one that straightens
your relationship, and a histogram cannot tell you which that is.

Printed as a markdown table, in identifiers only, so one artifact serves both
the English and the Portuguese page.
"""

import numpy as np
from scipy.stats import skew
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import (FunctionTransformer, PowerTransformer,
                                   QuantileTransformer, StandardScaler)

N = 800

rng = np.random.default_rng(23)
x = np.exp(rng.normal(3.0, 1.0, N))            # a lognormal column: heavily right-skewed
X = x.reshape(-1, 1)

truths = {
    "y~log(x)": 2.0 * np.log(x),
    "y~x": 0.05 * x,
    "y~sqrt(x)": 0.6 * np.sqrt(x),
}
transforms = {
    "identity": FunctionTransformer(),
    "log1p": FunctionTransformer(np.log1p),
    "sqrt": FunctionTransformer(np.sqrt),
    "yeo_johnson": PowerTransformer(method="yeo-johnson"),
    "quantile_normal": QuantileTransformer(
        output_distribution="normal", n_quantiles=400, random_state=0),
}

targets = {name: base + 0.25 * np.std(base) * rng.normal(size=N) for name, base in truths.items()}
cv = KFold(5, shuffle=True, random_state=0)

print("| `transform` | `skew_after` | `" + "` | `".join(targets) + "` |")
print("|---|---:|" + "---:|" * len(targets))
for name, transform in transforms.items():
    after = skew(transform.fit_transform(X).ravel())
    scores = []
    for y in targets.values():
        model = make_pipeline(transform, StandardScaler(), LinearRegression())
        scores.append(cross_val_score(model, X, y, cv=cv, scoring="r2").mean())
    best = max(scores)
    cells = " | ".join(f"**{s:.3f}**" if s == best else f"{s:.3f}" for s in scores)
    print(f"| `{name}` | {after:.2f} | {cells} |")

"""Binning buys a linear model what it lacks, and takes from a network what it had.

One input, uniform on [0, 10]; a target that follows a sine of it, so the
relationship is smooth and strongly non-monotone. The feature is given to each
model raw, and then cut into 4, 10 and 32 quantile bins and one-hot encoded.

A straight line cannot follow a sine at all — R² of about zero — and binning
rescues it completely, because a one-hot block of bins turns a linear model
into a piecewise-constant one. A network needed no rescuing: it fits the curve
from the raw column, and every bin boundary you impose is resolution taken away
from it.

Most preprocessing advice about discretization was written for the first row.

Printed as a markdown table, in identifiers only, so one artifact serves both
the English and the Portuguese page.
"""

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score
from sklearn.neural_network import MLPRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import KBinsDiscretizer, StandardScaler

N, BIN_COUNTS = 2000, (4, 10, 32)

rng = np.random.default_rng(7)
X = rng.uniform(0, 10, N).reshape(-1, 1)
y = 3 * np.sin(1.6 * X.ravel()) + 0.3 * rng.normal(size=N)

models = {
    "linear_regression": lambda: LinearRegression(),
    "small_network": lambda: MLPRegressor(hidden_layer_sizes=(32,), max_iter=3000, random_state=0),
}
cv = KFold(5, shuffle=True, random_state=0)
score = lambda pipe: cross_val_score(pipe, X, y, cv=cv, scoring="r2").mean()

print("| `model` | `raw_column` | `" + "` | `".join(f"{b}_bins" for b in BIN_COUNTS) + "` |")
print("|---|---:|" + "---:|" * len(BIN_COUNTS))
for name, make in models.items():
    row = [score(make_pipeline(StandardScaler(), make()))]
    for bins in BIN_COUNTS:
        cut = KBinsDiscretizer(n_bins=bins, encode="onehot-dense", strategy="quantile",
                               quantile_method="averaged_inverted_cdf")
        row.append(score(make_pipeline(cut, make())))
    best = max(row)                                   # mark whichever actually won
    cells = " | ".join(f"**{v:.3f}**" if v == best else f"{v:.3f}" for v in row)
    print(f"| `{name}` | {cells} |")