from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from math import erf, sqrt, pi
OUT=Path(__file__).resolve().parent
plt.rcParams.update({'font.size':12,'axes.spines.top':False,'axes.spines.right':False})
def pdf(x,mu=10,sigma=2): return np.exp(-.5*((np.asarray(x)-mu)/sigma)**2)/(sigma*sqrt(2*pi))
def cdf(x,mu=10,sigma=2): return .5*(1+erf((x-mu)/(sigma*sqrt(2))))
x=np.linspace(3,17,1500)
fig,axs=plt.subplots(1,2,figsize=(12,4.5),layout='constrained')
for ax in axs:
    ax.plot(x,pdf(x),color='#2463a6',lw=2)
    ax.set(xlabel='Valor x',ylabel='Densidade f(x)',ylim=(0,.23))
    ax.grid(alpha=.2)
axs[0].vlines(10,0,pdf(10),color='#bb4d20',lw=2)
axs[0].scatter([10],[pdf(10)],color='#bb4d20')
axs[0].set_title('Altura: f(10) ≈ 0,1995')
axs[0].text(3.7,.21,'P(X = 10) = 0')
a,b=9.5,10.5
xx=np.linspace(a,b,400)
axs[1].fill_between(xx,0,pdf(xx),color='#2463a6',alpha=.3,label='Área: probabilidade')
axs[1].plot([a,a,b,b,a],[0,pdf(10),pdf(10),0,0],color='#bb4d20',ls='--',label='Retângulo: aproximação')
axs[1].set_title('Intervalo de largura Δx = 1')
axs[1].legend(fontsize=10,loc='upper right')
fig.savefig(OUT/'01_altura_area.png',dpi=170)
plt.close(fig)
dados=np.array([8.,10.,12.]);mus=np.linspace(5,15,801)
logL=-.5*np.sum(((dados[None,:]-mus[:,None])/2)**2,axis=1)-3*np.log(2*sqrt(2*pi))
L=np.exp(logL)
fig,axs=plt.subplots(1,2,figsize=(12,4.5),layout='constrained')
axs[0].plot(mus,L,lw=2,label='Densidade conjunta L(μ)')
for dx in [1,.5,.1]:axs[0].plot(mus,L*dx**3,label=f'L(μ) × ({dx:g})³',ls='--')
axs[0].set(xlabel='Média candidata μ',ylabel='Valor da expressão',title='A largura altera a escala')
axs[0].legend(fontsize=10)
axs[1].plot(mus,L/L.max(),lw=3,label='Todas as curvas normalizadas')
for dx in [1,.5,.1]:
    exact=np.array([np.prod([cdf(v+dx/2,mu)-cdf(v-dx/2,mu) for v in dados]) for mu in mus])
    axs[1].plot(mus,exact/exact.max(),ls='--',label=f'Probabilidade exata: Δx = {dx:g}')
axs[1].axvline(10,color='#bb4d20',ls=':')
axs[1].set(xlabel='Média candidata μ',ylabel='Valor relativo ao próprio máximo',title='Máximo em μ = 10 neste exemplo')
axs[1].legend(fontsize=9)
for ax in axs:ax.grid(alpha=.2)
fig.savefig(OUT/'02_largura_maximo.png',dpi=170)
plt.close(fig)
for dx in [1,.1,.01]:
    exact=cdf(10+dx/2)-cdf(10-dx/2)
    print(dx,exact,pdf(10)*dx)
for mu in [6,10,14]:
    print('mu',mu,'L',np.prod(pdf(dados,mu)),'aprox dx=.1',np.prod(pdf(dados,mu))*.1**3)