%% Prueba hipótesis de media (unilateral)
% © 2019, Antonio Sala Piqueras. Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo en: http://personales.upv.es/asala/YT/V/thmdml.html
% Este código ejecutó correctamente en Matlab R2018b.
%% Motivación/objetivos
% Comprender el significado de la teoría en un ejemplo de detección de que una variable está, en media, por encima de un umbral, con un 99% de probabilidad. Esto tendría aplicación a "detección de fallos".
%% Contenido
% Table of Contents
% Motivación/objetivos
%     Contenido
% Modelado
% Simulación de proceso de toma de muestras y detección
% Gráficas de resultados
%
%% Modelado
% Por las características de un proceso, se conoce su desviación típica (2 unidades).
% Se desea hacer un test para detectar si la media supera un valor de 10 unidades y disparar una alarma cuando ocurra.

media_umbral=10;sigma=2;

media_real=10.6; %esto es para simular, obviamente es desconocido para mí

%% Simulación de proceso de toma de muestras y detección
% Empezamos a tomar muestras y calcular:

sum_muestras=0;
MaxMuestras=500;
grafZ=zeros(1,MaxMuestras);grafY=zeros(1,MaxMuestras);grafMedia=zeros(1,MaxMuestras);
for nmuestras=1:MaxMuestras
    y=randn()*sigma+media_real;%generamos una muestra nueva
    
    sum_muestras=sum_muestras+y;
    mediamuestral=sum_muestras/nmuestras;
    desvtipicamediamuestral=sigma/sqrt(nmuestras); 
    Z=(mediamuestral-media_umbral)/desvtipicamediamuestral; %esto debería ser N(0,1)
    
    grafY(nmuestras)=y; grafMedia(nmuestras)=mediamuestral;  grafZ(nmuestras)=Z; %gráficas
end
plot([grafY;grafMedia]'), grid on, legend('Muestras','media')

%% Gráficas de resultados
% Objetivo del test: equivocarse menos de un 1% de las veces que rechacemos que la media está por debajo de media_umbral, esto es, que afirmemos que está por encima.

alpha_unilateral=0.99; 
norminv(alpha_unilateral)
%con este código hubiese salido lo mismo limiteZ=sqrt(2)*erfcinv(2*(1-alpha_unilateral))
plot(grafZ), grid on
yline(limiteZ,'-.m','Umbral Z');
title('Estadístico Z conforme se acumulan muestras')
tiempodeteccion=find(grafZ>limiteZ, 1 )
xline(tiempodeteccion,'-.r','Disparo alarma');

% El gráfico más intuitivo es la cota de la media muestral a partir de la cual descarto la hipótesis nula: media >Z_{limite}*\frac{sigma}{\sqrt{n}}+\mu

plot([limiteZ*sigma./sqrt(1:MaxMuestras)+media_umbral;grafMedia]'), legend('umbral','media datos')
xline(tiempodeteccion,'-.r','Disparo alarma'); yline(media_umbral,'-m','Media hipótesis'); grid on

% Dibujemos la probabilidad (logarítmica) de que un Z mayor o igual que el observado se produzca con una media real menor de 10 (hipótesis 0 cierta):

grafPvalue=1-normcdf(grafZ)
%grafPvalue=0.5-0.5*erf(grafZ/sqrt(2)) sintaxis alternativa para lo mismo.
semilogy(grafPvalue), grid on
yline(1-alpha_unilateral,"-.m",'Umbral 1-alpha');
xline(tiempodeteccion,'-.r','Disparo alarma');
title('Probabilidad falsa alarma')

%