import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from io import StringIO

# The target of California housing drawn three ways: 10 bins, 50 bins, and a
# Gaussian kernel density estimate (KDE) whose bandwidth comes from the usual
# rule of thumb for Gaussian kernels, h = 1.06 * std * n**(-1/5).
URL = "https://raw.githubusercontent.com/hsandmann/biblio/refs/heads/main/ml/aula02/housing.csv"
v = pd.read_csv(URL)['median_house_value'].to_numpy() / 1000       # in US$ thousand

BLUE, ORANGE, RED = '#58a6ff', '#f0883e', '#ff7b72'
plt.rcParams.update({'text.color': 'gray', 'axes.labelcolor': 'gray', 'axes.edgecolor': 'gray',
                     'xtick.color': 'gray', 'ytick.color': 'gray'})
fig, axes = plt.subplots(1, 3, figsize=(10, 3.3))

for ax, bins in zip(axes[:2], [10, 50]):
    counts, edges, bars = ax.hist(v, bins=bins, color=BLUE, edgecolor='none', rwidth=0.9)
    bars[-1].set_facecolor(ORANGE)
    ax.set_title(f'{bins} bins (width {edges[1] - edges[0]:.1f}k)', fontsize=9.5)
    ax.annotate(f'last bin: {counts[-1]:.0f}\nneighbor: {counts[-2]:.0f}',
                xy=(edges[-2] + (edges[1] - edges[0]) / 2, counts[-1]), xytext=(0.52, 0.80),
                textcoords='axes fraction', fontsize=8, color=ORANGE,
                arrowprops=dict(arrowstyle='->', color=ORANGE, lw=0.8))
    ax.set_ylabel('districts', fontsize=8.5)

# KDE by hand: the average of one Gaussian bump per district
h = 1.06 * v.std(ddof=1) * len(v) ** (-1 / 5)                     # rule of thumb
grid = np.linspace(0, 650, 651)
dens = np.zeros_like(grid)
for chunk in np.array_split(v, 20):
    dens += np.exp(-0.5 * ((grid[:, None] - chunk[None, :]) / h) ** 2).sum(axis=1)
dens /= len(v) * h * np.sqrt(2 * np.pi)
beyond = dens[grid > 500.001].sum() * (grid[1] - grid[0])

ax = axes[2]
ax.hist(v, bins=50, density=True, color=BLUE, alpha=0.35, edgecolor='none', rwidth=0.9)
ax.plot(grid, dens, color=ORANGE, lw=1.8)
ax.fill_between(grid, dens, where=grid > 500.001, color=RED, alpha=0.5, lw=0)
ax.axvline(500.001, color='gray', ls='--', lw=0.8)
ax.set_title(f'KDE, bandwidth {h:.1f}k (rule of thumb)', fontsize=9.5)
ax.set_ylim(0, dens.max() * 1.5)
ax.text(0.97, 0.95, f'{beyond:.1%} of the curve lies above\nthe 500k ceiling, where no\ndistrict exists',
        transform=ax.transAxes, ha='right', va='top', fontsize=8, color=RED)
ax.set_ylabel('density', fontsize=8.5)
ax.set_yticks([])

for ax in axes:
    ax.set_xlabel('median_house_value (US$ thousand)', fontsize=8.5)
    ax.tick_params(labelsize=8)
fig.suptitle('The same 20,640 values: the 500,001 ceiling looks like a fat tail, a spike, or a bump that spills past it', fontsize=10)

plt.tight_layout()

buffer = StringIO()
plt.savefig(buffer, format="svg", transparent=True)
print(buffer.getvalue())

plt.close()