x = linspace(-pi, pi, 4000);
Ns = [10, 20, 50, 100];

figure

for k = 1:length(Ns)
    
    N = Ns(k);
    
    A = randn(1,N);
    B = randn(1,N);
    
    f = zeros(size(x));
    for n = 1:N
        f = f + A(n)*cos(n*x) + B(n)*sin(n*x);
    end
    
    subplot(2,2,k)
    plot(x, f, 'DisplayName', ['N = ' num2str(N)])
    xlim([-pi pi])
    ylim([-30 30])
    legend show
    grid on
    
end

sgtitle('Serie de Fourier con coeficientes N(0,1)')

import numpy as np
import matplotlib.pyplot as plt

# -------------------------
# Parámetros del modelo
# -------------------------
L = np.pi
N = 30
seed = 0

M = 3000
ells = np.array([1e-3, 2e-3, 5e-3, 1e-2, 2e-2, 5e-2, 0.1, 0.2, 0.4, 0.7, 1.0, L])

# Malla adaptativa (densidad constante)
points_per_unit = 120
Jmin, Jmax = 11, 500

rng = np.random.default_rng(seed)

n = np.arange(1, N + 1)

def sample_coeffs():
    A0 = rng.normal(0.0, 1.0)
    A = rng.normal(0.0, 1.0, size=N)
    B = rng.normal(0.0, 1.0, size=N)
    return A0, A, B

def eval_f_on_grid(xgrid, A0, A, B):
    """
    Evalúa la serie de Fourier en los puntos xgrid
    
    Args:
        xgrid: array de puntos donde evaluar (J,)
        A0: término constante (escalar)
        A: coeficientes coseno (N,)
        B: coeficientes seno (N,)
    
    Returns:
        f: array (J,) con los valores de la función
    """
    # Crear matriz theta: (J, N) donde theta[j, n-1] = (π/L) * xgrid[j] * n
    theta = (np.pi / L) * np.outer(xgrid, n)
    
    # Calcular la serie de Fourier sin usar @
    # Para cada punto x, sumamos sobre n: A_n*cos(theta_n) + B_n*sin(theta_n)
    cos_terms = np.cos(theta) * A  # Multiplicación elemento a elemento, luego se suma sobre n
    sin_terms = np.sin(theta) * B
    
    # Sumar sobre el eje de los términos (axis=1)
    f = A0 + np.sum(cos_terms, axis=1) + np.sum(sin_terms, axis=1)
    
    return f

# -------------------------
# Monte Carlo para cada longitud
# -------------------------
lengths = ells
p_hats = []

for ell, length in zip(ells, lengths):
    a, b = -ell, ell

    Jell = int(np.clip(np.ceil(points_per_unit * length), Jmin, Jmax))
    xgrid = np.linspace(a, b, Jell)

    success = 0
    for _ in range(M):
        A0, A, B = sample_coeffs()
        f = eval_f_on_grid(xgrid, A0, A, B)
        if np.all(f >= 0):
            success += 1

    p_hat = success / M
    p_hats.append(p_hat)
    print(f"longitud={length:.4g} (ℓ={ell:.4g}) | J={Jell:3d} | p_hat={p_hat:.4f} | {success}/{M}")

p_hats = np.array(p_hats, dtype=float)

# -------------------------
# Plot: Probabilidad vs Longitud
# -------------------------
plt.figure(figsize=(8, 5))
plt.plot(lengths, p_hats, marker="o", linewidth=1.8)

plt.title("Probabilidad estimada vs L")
plt.xlabel("L")
plt.ylabel(r"$\mathbb{P}(f_\sigma(x) ≥ 0, \forall x, \in [-L,L])$")
plt.grid(True, alpha=0.3)


plt.tight_layout()
plt.show()

%Fourier N(0,1/n)

x = linspace(-pi, pi, 4000);
Ns = [100, 200, 500, 1000];

figure

for k = 1:length(Ns)
    
    N = Ns(k);
    n = 1:N;
    
    % Desviación típica = 1/n
    A = (1./n) .* randn(1,N);
    B = (1./n) .* randn(1,N);
    
    f = zeros(size(x));
    for j = 1:N
        f = f + A(j)*cos(j*x) + B(j)*sin(j*x);
    end
    
    subplot(2,2,k)
    plot(x, f, 'DisplayName', ['N = ' num2str(N)])
    xlim([-pi pi])
    ylim([-5 5])
    legend show
    grid on
    
end

sgtitle('Serie de Fourier con A_n,B_n ~ N(0, 1/n)')


%Fourier N(0,1/n^2)
x = linspace(-pi, pi, 4000);
Ns = [100, 200, 500, 1000];

figure

for k = 1:length(Ns)
    
    N = Ns(k);
    n = 1:N;
    
    % Desviación típica = 1/n
    A = (1./n.^2) .* randn(1,N);
    B = (1./n.^2) .* randn(1,N);
    
    f = zeros(size(x));
    for j = 1:N
        f = f + A(j)*cos(j*x) + B(j)*sin(j*x);
    end
    
    subplot(2,2,k)
    plot(x, f, 'DisplayName', ['N = ' num2str(N)])
    xlim([-pi pi])
    ylim([-5 5])
    legend show
    grid on
    
end

sgtitle('Serie de Fourier con A_n,B_n ~ N(0, 1/n^2)')