"""Exemplos de máxima verossimilhança. Requer numpy e matplotlib."""
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

OUT = Path(__file__).resolve().parent
plt.rcParams.update({'font.size': 12, 'axes.spines.top': False, 'axes.spines.right': False})
# Exemplo 1: sequência de 10 lançamentos independentes, 7 caras e 3 coroas.
p = np.linspace(0, 1, 1001)
L = p**7 * (1-p)**3
p_estimado = 7/10
fig, ax = plt.subplots(figsize=(9, 4.7), layout='constrained')
ax.plot(p, L/L.max(), lw=3, color='#2463a6')
ax.axvline(p_estimado, ls='--', color='#bb4d20')
ax.scatter([p_estimado], [1], color='#bb4d20', s=65, zorder=3)
ax.annotate('Máximo em p = 0,70', (.7, 1), xytext=(.43, 1.08), arrowprops={'arrowstyle':'->'})
ax.set(xlabel='p: probabilidade de cara', ylabel='Verossimilhança relativa ao máximo',
       title='Dados fixos: 7 caras e 3 coroas — variamos p', ylim=(0, 1.2))
ax.grid(alpha=.2)
fig.savefig(OUT/'verossimilhanca_moeda.png', dpi=180)
plt.close(fig)
# Exemplo 2: normal com desvio padrão conhecido igual a 2.
dados = np.array([8., 10., 12.])
sigma = 2.
mu_estimado = dados.mean()
x = np.linspace(2, 18, 1000)
def densidade(x, mu):
    return np.exp(-.5*((x-mu)/sigma)**2)/(sigma*np.sqrt(2*np.pi))
fig, axes = plt.subplots(2, 1, figsize=(9, 8), layout='constrained')
for mu, color in [(6, '#888888'), (10, '#2463a6'), (14, '#bb4d20')]:
    axes[0].plot(x, densidade(x, mu), label=f'μ = {mu}', color=color, lw=2)
    axes[0].scatter(dados, densidade(dados, mu), color=color, s=30)
axes[0].plot(dados, np.zeros(3), '|', color='black', ms=17, label='Dados: 8, 10, 12')
axes[0].set(title='Modelos candidatos: σ = 2, diferentes médias', xlabel='Valor x', ylabel='Densidade f(x | μ)')
axes[0].legend(loc='upper right', fontsize=10)
mus = np.linspace(2, 18, 1001)
logL = -.5*np.sum(((dados[None,:]-mus[:,None])/sigma)**2, axis=1)
relativa = np.exp(logL-logL.max())
axes[1].plot(mus, relativa, color='#2463a6', lw=3)
axes[1].axvline(mu_estimado, color='#bb4d20', ls='--', label='Máximo: μ = 10')
axes[1].set(title='Produto das densidades nos três valores observados', xlabel='Média candidata μ', ylabel='Verossimilhança relativa ao máximo')
axes[1].legend()
for ax in axes: ax.grid(alpha=.2)
fig.savefig(OUT/'verossimilhanca_normal.png', dpi=180)
plt.close(fig)
print(f'Estimativa para moeda: {p_estimado:.2f}')
print(f'Estimativa para média normal: {mu_estimado:.2f}')