# -*- coding: utf-8 -*-
"""
Created on Fri Oct  3 12:41:34 2025

@author: Moritz Romeike
"""

# ------------------------------------------------------------------------
# Programmcode: Theil-Index & Atkinson-Index + Histogramm (Python)
# ------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt

# Marktanteile in Prozent (wie im R-Skript)
lieferanten_anteile_pct = np.array([
    12.5, 10.3, 9.8, 8.5, 7.6, 6.9, 6.5, 6.2, 5.8, 5.2,
    4.9, 4.6, 4.2, 3.9, 3.7, 3.5, 3.3, 3.1, 2.9, 2.7,
    2.5, 2.3, 2.1, 1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3,
    1.2, 1.1, 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3,
    0.3, 0.3, 0.3, 0.2, 0.2, 0.2, 0.1, 0.1, 0.1, 0.1
], dtype=float)

# Auf Anteile (0..1) normieren – wie in R: x / sum(x)
p = (lieferanten_anteile_pct / lieferanten_anteile_pct.sum()).astype(float)

# --- Theil-Index ---------------------------------------------------------
# T = sum_i p_i * ln(p_i / mu), mit mu = 1/n (bei Anteilen mit Summe 1)
# Vereinfachung: T = sum_i p_i * ln(n * p_i)
def theil_index(p: np.ndarray) -> float:
    p = np.asarray(p, dtype=float)
    p = p[p > 0]                       # 0-Anteile ignorieren, sonst log(0)
    n = p.size
    return float(np.sum(p * np.log(n * p)))

# --- Atkinson-Index ------------------------------------------------------
# A(ε) = 1 - ( (1/n * sum_i p_i^(1-ε))^(1/(1-ε)) ) / μ,    μ = 1/n
# Für ε = 0.5 (wie im R-Skript)
def atkinson_index(p: np.ndarray, epsilon: float = 0.5) -> float:
    p = np.asarray(p, dtype=float)
    n = p.size
    mu = p.mean()                      # = 1/n bei normierten Anteilen
    if epsilon == 1.0:
        # Grenzfall ε→1 (log-Form) – hier nicht benötigt
        raise NotImplementedError("Atkinson mit ε=1 ist separat zu behandeln.")
    term = (np.mean(np.power(p, 1.0 - epsilon))) ** (1.0 / (1.0 - epsilon))
    A = 1.0 - (term / mu)
    return float(A)

T = theil_index(p)
A = atkinson_index(p, epsilon=0.5)

print(f"Theil-Index:     {T:.3f}")
print(f"Atkinson-Index:  {A:.3f}  (ε = 0.5)")

# --- Histogramm der Prozentanteile --------------------------------------
plt.figure(figsize=(8,4))
plt.hist(lieferanten_anteile_pct, bins=np.arange(0, np.ceil(lieferanten_anteile_pct.max())+1, 1),
         edgecolor="black")
plt.title("Verteilung der Lieferantenanteile")
plt.suptitle(f"Theil-Index: {T:.3f} | Atkinson-Index (ε=0.5): {A:.3f}", y=0.98, fontsize=10)
plt.xlabel("Marktanteil (%)")
plt.ylabel("Häufigkeit")
plt.tight_layout()
plt.show()
# ------------------------------------------------------------------------