%% Fórmula de Bayes: caso de estudio "tigre oculto"
% © 2024, Antonio Sala Piqueras. Universitat Politècnica de València. Todos los derechos reservados.
% Presentaciones en vídeo:
% http://personales.upv.es/asala/YT/V/tiger6br.html
% http://personales.upv.es/asala/YT/V/tiger7rbs.html
%
% Este código ejecutó en Matlab R2023b
% Objetivos: Ilustrar la fórmula de Bayes en el caso "Tigre oculto detrás de dos puertas", donde se deben acumular mediciones ruidosas sobre la ubicación del tigre para mejorar las posibilidades de supervivencia (evaluar correctamente dónde está el tigre antes de "actuar" abriendo una puerta).
%
% El "decidir" la acción óptima NO será objetivo de este material, centrado sólo en actualizar la "creencia" (belief en inglés) conforme se reciben observaciones.
%
%
%
%
%
%
% Tabla de Contenidos
% Planteamiento del problema
%     Construcción/Interpretación de tabla de probabilidad
% Caso 0: Se cierra reja del pasillo central y escucho UN rugido
%     Actualización de la creencia (posterior): fórmula de Bayes
% Revisión fórmula de Bayes
% Fórmula de Bayes para DOS rugidos (NO recursiva)
%     Fórmula de Bayes para número arbitrario de rugidos (NO recursiva)
% Fórmula de Bayes recursiva
% Simulación de formula de Bayes recursiva (tigre NO se mueve)
%% Planteamiento del problema
% Tenemos una situación donde un empleado de zoológico debe "adivinar" dónde está el tigre, para escapar por el lado donde no esté sin ser devorado por él:
%
% El tigre puede estar a la izquierda o a la derecha, y sus rugidos pueden escucharse "más fuerte" por la puerta izquierda o por la derecha. Simplificamos a:
% Medidas (observaciones) binarias:  m\in \{hear_left (1, HL), hear_right (2, HR) \} .
% Estados tigre binarios: t\in \{tiger_left (1, TL), tiger_right (2, TR) \}.

hear_L=1;  hear_R=2; %El número únicamente significa posición en listas
tiger_L=1; tiger_R=2; %El número únicamente significa posición en listas

%% Construcción/Interpretación de tabla de probabilidad
% Definimos matriz PM, donde PM(m,e) es la probabilidad de medir "m" cuando el estado del tigre es "t". O sea, es una tabla de probabilidad condicional.
% PM= [ p(HL|TL)  p(HL|TR) ;
%           p(HR|TL)  p(HR|TR)]
% Filas representan "observaciones", columnas representan "estados".
% Propiedades: COLUMNAS suman 1. No tiene por qué ser simétrica.

%PM=[.99 .9;.01 .1]; %non symmetric... different measurement noise at each state... most times I hear "left" even if it is "right".
%PM=[.8 .1;.2 .9];
%PM=[.8 .2;.2 .8];
PM=[.6 .4;.4 .6];
%PM=[.5 .0001;.5 .9999];
PM

% Será nuestro "modelo del mundo tigre" y haremos "inferencia estadística" a partir de él (estimar dónde está el tigre, decidir qué puerta abrir, decidir cuántos rugidos escuchar hasta estar "suficientemente seguros" para abrir la puerta  -- eso depende del miedo que tenga, que se modela con funciones de utilidad/coste adicionales--)...
%
%% Caso 0: Se cierra reja del pasillo central y escucho UN rugido
% La creencia ("belief") se denotará con b. Será la probabilidad "a priori" de tiger_right. Obviamente, la probabilidad a priori de tiger_left será 1-b.
% Si lo hacemos "en letra" lo de la marginal, tenemos:

syms b real % b := p(TR) a priori
prob_hearL(b)= PM(hear_L,tiger_R)*b + PM(hear_L,tiger_L)*(1-b);
vpa(prob_hearL)
prob_hearR(b)= PM(hear_R,tiger_R)*b + PM(hear_R,tiger_L)*(1-b);
vpa(prob_hearR)

% En notacion matricial resulta "matriz de probabilidades condicionales" x "vector de probabilidades a priori":

vpa(  PM*[(1-b); b]  , 4)
fplot([prob_hearL;prob_hearR],[0 1]), grid on, ylim([0 1])
legend("p(Hear Left)","p(Hear Right)") %OJO interpretacion, b no es "evento"...
xlabel("creencia a priori b :=  p(Tiger Right)"), title("Prob. Marginales: hearleft, hearright (en func. de p(TR))")

%%
%% Actualización de la creencia (posterior): fórmula de Bayes
% Si midiese hear_left, lo podemos escribir a partir de la tabla "condicional" y la probabilidad "a priori" del lado donde se situa el tigre:
% p(TR|hearleft)=\frac{p(hearleft|TR)\cdot p(TR)}{p(hearleft|TR)\cdot p(TR)+p(hearleft|TL)\cdot p(TL)}=\frac{p(hearleft|TR)\cdot b}{p(hearleft|TR)\cdot b+p(hearleft|TL)\cdot (1-b)}
% Antes lo hemos hecho para un valor de b numérico dado... ahora lo escribo en simbólico:

posteriori_ifHL(b) = PM(hear_L,tiger_R)*b/prob_hearL(b); %dividimos por la marginal
vpa(posteriori_ifHL)

% Si midiese hear_right,
% p(TR|hearright)=\frac{p(hearright|TR)\cdot p(TR)}{p(hearright|TR)\cdot p(TR)+p(hearright|TL)\cdot p(TL)}=\frac{p(hearright|TR)\cdot b}{p(hearright|TR)\cdot b+p(hearright|TL)\cdot (1-b)}

posteriori_ifHR(b) = PM(hear_R,tiger_R)*b/prob_hearR(b);%dividimos por la marginal
vpa(posteriori_ifHR)

% Representación gráfica:
% El proceso "conceptual" es que a partir de una "creencia" (probabilidad a priori) y una observación, calculo una "creencia actualizada" (probabilidad a posteriori): b_+=Bayes(b,m). b\in[0,1], m\in\{HL,HR\}. Gráficamente:

fplot([posteriori_ifHL;posteriori_ifHR],[0 1],LineWidth=2.5), grid on
hold on, fplot(b,[0 1],'-.',LineWidth=1), hold off
xlabel('b (Prob. right)'), ylabel('b+ (posterior)')
legend("New belief if HearLeft","New belief if HearRight",'identity',Location="best")

%
% No hay un punto "fijo" en el que la creencia nueva sea igual a la antigua excepto la creencia total a la derecha (b=1) o izquierda (b=0). Si estoy "100% seguro" de algo, experimentar no va a aportar nueva información.
%
%% Revisión fórmula de Bayes
% Llamemos, por compactar notación
% PM= \left(\begin{array}{cc}p(HL|TL) &p(HL|TR)
%  \\ p(HR|TL) &p(HR|TR)\end{array}\right)= \left(\begin{array}{cc}q &r
%  \\ 1-q &1-r\end{array}\right)
% a las condicionales HL HR (filas) dado TL y TR (columnas), respectivamente. Si b=p(TR), la fórmula de Bayes dirá que, según observe HL o HR, aplicaré una de estas dos fórmulas:
% b_{+}=p(TR\,|\,HL)=\frac{r\cdot b}{q\cdot (1-b)+r\cdot b}
% b_{+}=p(TR\,|\,HR)=\frac{(1-r)\cdot b}{(1-q)\cdot (1-b)+(1-r)\cdot b}
%
%
%% Fórmula de Bayes para DOS rugidos (NO recursiva)
% Si me espero hasta escuchar DOS rugidos (en abstracto observación o_1 \in \{HL,HR\} y o_2 \in \{HL,HR\}), suponiendo "independencia estadística" (bueno, en rigor se denomina "independencia condicional" porque estamos hablando de productos de probabilidades condicionales), entonces las nuevas condicionales respecto a TL son:
%
% p((HL,HL)|TL) es    p(HL|TL)·p(HL|TL) = q^2,
% p((HL,HR)|TL)=p((HR,HL)|TL)= p(HL|TL)p(HR|TL) = q\cdot (1-q) ,
% p((HR,HR)|TL)=p(HR|TL)p(HR|TL)=(1-q)^2
%
% y lo mismo con las condicionales respecto a TR:
% p((HL,HL)|TR) es r^2,
% p((HL,HR)|TR)=P((HR,HL)|TR)=r\cdot (1-r) ,
% p((HR,HR)|TR)=(1-r)^2
%
%
%
% Con lo cual mi nueva tabla condicional para DOS rugidos para aplicar Bayes será:
% PM2 = \left(\begin{array}{cc}p(HLL|TL)&p(HLL|TR)
% \\ p(HLR|TL) &p(HLR|TR) \\ p(HRL|TL) &p(HRL|TR) \\ p(HRR|TL) &p(HRR|TR)\end{array}\right) = \left(\begin{array}{cc}q^2 &r^2
% \\ q(1-q) &r(1-r) \\ q(1-q) &r(1-r) \\ (1-q)^2 &(1-r)^2\end{array}\right)
% *HLL es una abreviatura de  ( o_1= HL, o_2=HL) y similar para las otras combinaciones de observaciones.
%
%
%
%
% Nota: El significado de o_2\wedge o_1 se refiere a observaciones en distinto instante de tiempo... No podemos observar o_1=HL, o_2=HR si se refiriesen al mismo instante, obviamente. Los problemas que estudian variables (aleatorias) en diferentes instantes de tiempo se dice que estudian la "dinámica" de un sistema.
% Por ello, formalmente, debería ser o_1 \in \{HL_{t=1},HR_{t=1}\} y o_2 \in \{HL_{t=2},HR_{t=2}\}, aunque omitiremos el subíndice temporal asociado a cada evento por simplicidad de notación.
%
% Fórmula de Bayes para la observación (HL_{t=1},HL_{t=2})... que abrevio a HLL:
% b_{++}=p(TR\,|\,HLL)=\frac{p(HLL \,|\, TR)p(TR)}{p(HLL \,|\, TL)p(TL)+p(HLL \,|\,TR)p(TR)}\\ =\frac{r^2\cdot b}{q^2\cdot (1-b)+r^2\cdot b}
% deberíamos completar todas las opciones:
% p(TR\,|\,HLL), p(TR\,|\,HLR), p(TR\,|\,HRL), p(TR\,|\,HRR)
% *Se deja como ejercicio.
% las p(TL|\dots) serán 1-la correspondiente de p(TR|\dots)...Esas no hacen falta guardar.
%
% Para 3 rugidos (sin que el tigre pueda cambiar de lado, reja central cerrada), suponiendo "independencia condicional", la tabla para la fórmula de Bayes cuando observemos algo sería:
% PM3=\left(\begin{array}{cc}p(HLLL|TL)&p(HLLL|TR)
% \\ p(HLLR|TL) &p(HLLR|TR) \\ p(HLRL|TL) &p(HLRL|TR) \\ p(HLRR|TL) &p(HLRR|TR) \\
% p(HRLL|TL)&p(HRLL|TR)
% \\ p(HRLR|TL) &p(HRLR|TR) \\ p(HRRL|TL) &p(HRRL|TR) \\ p(HRRR|TL) &p(HRRR|TR) \end{array}\right) = \left(\begin{array}{cc}q^3 &r^3
% \\ q^2(1-q) &r^2(1-r) \\ q^2(1-q) &r^2(1-r) \\ q(1-q)^2 &r(1-r)^2 \\
% q^2(1-q) &r^2(1-r)
% \\ q(1-q)^2 &r(1-r)^2 \\ q(1-q)^2 &r(1-r)^2 \\ (1-q)^3 &(1-r)^3\end{array}\right)
%
% Reducción de tabla "sólo contar"
% Si los eventos son "contar 2L y 1R" (el orden no importa), entonces podemos resumir la tabla en:
% PM3b=\left(\begin{array}{cc}p(HL^3|TL)&p(HL^3|TR)
% \\ p(HL^2R|TL) &p(HL^2|TR)\\  p(HLR^2|TL) &p(HLR^2|TR)
% \\ p(HR^3|TL) &p(HR^3|TR) \end{array}\right) = \left(\begin{array}{cc}q^3 &r^3
% \\ 3q^2(1-q) &3r^2(1-r)\\ 3q(1-q)^2 &3r(1-r)^2 \\
%  (1-q)^3 &(1-r)^3\end{array}\right)
%
%
%
%
%  Las fórmulas de Bayes acabarán dando el mismo resultado: dividir cada elemento de la fila asociada a la observación experimentada por la suma de la fila resulta lo mismo esté o no multiplicada por 3. Ello hace que, en este tipo de problemas, el orden de las observaciones sea obviado en la mayoría de casos, porque así se reduce el número de posibles eventos y las tablas son más pequeñas (¡imagina lo grande que sería todo con 20 rugidos!).
%
%% Fórmula de Bayes para número arbitrario de rugidos (NO recursiva)
% Con la combinatoria adecuada, todo se puede resumir en que si de un total de m+n rugidos escuchamos m por la izquierda y n por la derecha, haríamos Bayes con:
% \left(\begin{array}{cc}\ p(HL^mR^n\,|\,TL)&p(HL^mR^n\,|\,TR)\
%  \end{array}\right)  \qquad \\ \qquad = \left(\begin{array}{cc}\ \frac{(m+n)!}{m!\,n!}q^m(1-q)^n \  &\frac{(m+n)!}{m!\,n!}r^m(1-r)^n \
% \end{array}\right)
% Los libros de estadística (busca distribución "binomial" en tus libros de texto) suelen hacerlo así.
%
%
%
%
%
% Dirán que las probabilidades condicionales tienen distribución binomial, en concreto:
% \bullet el número de rugidos "izquierda" de un total de "Tot" observaciones (condicionadas a TL) es una variable aleatoria con distribución binomial de parámetros (Tot,q).
% \bullet el número de rugidos "izquierda" de un total de "Tot" observaciones (condicionadas a TR) es una variable aleatoria con distribución binomial de parámetros (Tot,r).
%
%
% Ejemplo numérico:
% Leamos la tabla de probabilidad condicional con la que comenzamos todo:

PM
q=PM(hear_L,tiger_L)
r=PM(hear_L,tiger_R)

% Si escucho rugidos "m" veces a la izquierda y "n" a la derecha, entonces:

m=4; %number of HL events, num. rugidos izquierda
n=4; %number of HR events, num. rugidos derecha
Tot=m+n; %total roaring events, total rugidos
p_LmRn_cond_TL = nchoosek(Tot,m) * q^m * (1-q)^n
p_LmRn_cond_TR = nchoosek(Tot,m) * r^m * (1-r)^n

% Con lo que aplicando Bayes, la probabilidad a posteriori de que el tigre esté a la izquierda o a la derecha será:

b=0.5; %probabilidad a priori de TR
p_TL_cond_LmRn = p_LmRn_cond_TL*(1-b) / ( p_LmRn_cond_TL*(1-b) + p_LmRn_cond_TR*b )
p_TR_cond_LmRn = p_LmRn_cond_TR*b / ( p_LmRn_cond_TL*(1-b) + p_LmRn_cond_TR*b)

%
%
%
%
%% Fórmula de Bayes recursiva
% ¿Podemos aplicar recursivamente la fórmula con la tabla de UN rugido original para obtener el mismo resultado ante varios rugidos? La respuesta es AFIRMATIVA. Veamos primero caso "dos" rugidos y luego el caso general.
% Esta es la fórmula para UN rugido, válida:
% b_{+}=p(TR\,|\,HL)=\frac{r\cdot b}{q\cdot (1-b)+r\cdot b}
% Conjetura... ¿Qué producirá esto?
% b_{++}=p(TR\,|\,HLL)=\frac{r\cdot b_+}{q\cdot (1-b_+)+r\cdot b_+}
% Vamos a sustituir:
% b_{++}=p(TL|LL)=\frac{r*\frac{r\cdot b}{q\cdot (1-b)+r\cdot b}}{q*\frac{q\cdot (1-b)}{q\cdot (1-b)+r\cdot b}+r*\frac{r\cdot b}{q\cdot (1-b)+r\cdot b}}
% Entonces, multiplicando numerador y denominador por q\cdot (1-b)+r\cdot b, queda:
% b_{++}=\frac{r^2b}{q^2(1-b)+r^2b}
% que coincide con la fórmula de Bayes y la tabla de 4 posibles observaciones PM2.
% Veamos el caso general.
% Teorema: en general si unas observaciones pasadas O_N=\{o_1,o_2,\dots,o_N\}\in\{HL,HR\}^N han dado lugar a una creencia (probabilidad a posteriori):
% b_{N}=p(TR|O_N)=\frac{p(O_N|TR)\cdot b}{p(O_N|TL)\cdot (1-b)+p(O_N|TR)\cdot b}
% entonces, ante una nueva observación o_{N+1}\in\{HL,HR\}:
% b_{N+1}=p(TR\,|\, O_{N+1})=\frac{p(o_{N+1}|TR)\cdot b_N}{p(o_{N+1}|TL)\cdot (1-b_N)+p(o_{N+1}|TR)\cdot b_N}
% siendo O_{N+1}=O_N\wedge o_{N+1}=\{o_1,o_2,\dots,o_N,o_{N+1}\}\in\{HL,HR\}^{N+1}, si la nueva observación es condicionalmente independiente de las anteriores.
%
% Demostracion: en efecto, sustituyendo la fórmula de b_N resulta:
% b_{N+1}=p(TR\,|\, O_N\wedge o_{N+1}) \\ \ =\frac{p(o_{N+1}|TR)\cdot \left(\frac{p(O_N|TR)\cdot b}{p(O_N|TL)\cdot (1-b)+p(O_N|TR)\cdot b}\right)}{p(o_{N+1}|TL)\cdot \left(\frac{p(O_N|TL)\cdot (1-b)}{p(O_N|TL)\cdot (1-b)+p(O_N|TR)\cdot b}\right)+p(o_{N+1}|TR)\cdot \left(\frac{p(O_N|TR)\cdot b}{p(O_N|TL)\cdot (1-b)+p(O_N|TR)\cdot b}\right)}
% y quitando denominadores:
% b_{N+1}=\frac{p(o_{N+1}|TR)p(O_N|TR)\cdot b}{p(o_{N+1}|TL)p(O_N|TL)\cdot (1-b)+p(o_{N+1}|TR)p(O_N|TR)\cdot b}
% con lo que si O_N y o_{N+1} son CONDICIONALMENTE INDEPENDIENTES, esto es, se verifica quep(o_{N+1}\,|\,TR)p(O_N\,|\,TR)=p(O_N\wedge o_{N+1}\,|\,TR)=p(O_{N+1}\,|\,TR),  resulta
% b_{N+1}=\frac{p(O_{N+1}\,|\,TR)\cdot b}{p(O_{N+1}\,|\,TL)\cdot (1-b)+p(O_{N+1}\,|\,TR)\cdot b}
% esto es, la fórmula de Bayes con la probabilidad del "evento conjunto" formado por las viejas observaciones junto a la nueva.
%
% En resumen, para eventos condicionalmente independientes, es lo mismo aplicar Bayes con la "tabla ampliada" que aplicar Bayes RECURSIVAMENTE.
%
% Utilidad práctica de fórmula Bayes recursiva:
% En algunos problemas no sabemos a priori cuántos datos (rugidos, observaciones, experimentos) necesitamos... Iremos acumulando "hasta que tengamos suficiente evidencia"... en la fórmula recursiva sólo manejamos "un único número que se actualiza con cada observación", sin números combinatorios, distribución binomial, etc.
%
% Esta es la idea principal detras de "filtros bayesianos recursivos", que se discuten en otros materiales. También, describir un modelo estadístico mediante tablas de probabilidad condicional entre variables condicionalmente independientes es la base de las "redes bayesianas".
%
%
%% Simulación de formula de Bayes recursiva (tigre NO se mueve)

Posterior_ifL=matlabFunction(posteriori_ifHL); %symbolic to numeric
Posterior_ifR=matlabFunction(posteriori_ifHR); %symbolic to numeric

real_tiger=tiger_L; %actual location of tiger

belief_ini=0.5; %initial belief, non-informative prior
N=2000; %number of samples per simulation, first sample is "prior", so N=observations+1
Nsimulations=100;
belief_plot=zeros(N,Nsimulations);
observ_plot=zeros(N,Nsimulations);
t_fin=zeros(Nsimulations,1);

for Nsim=1:Nsimulations
endflag=false;
belief=belief_ini;
belief_plot(1,Nsim)=belief_ini;
observ_plot(1,Nsim)=NaN; %prior does not have observation

%simulate, if so wished, a given sequence of observations:
%obs_test=[hear_L hear_L hear_L hear_L hear_R hear_R hear_R hear_R];

    for i=2:N
        % Measurement simulation (uncertain)
        % Simulate Bernouilli (unfair coin) distribution from uniform
        if rand()<PM(hear_L,real_tiger) %rand returns uniform in [0,1]
            obs=hear_L;
        else
            obs=hear_R;
        end

        %obs=obs_test(i-1); %If I wish to simulate a given sequence of
        % observations instead of generating them randomly

        % Recursive Bayes Filter start
        if obs==hear_L %hear left
            belief=Posterior_ifL(belief); %belief update based on measuring LEFT
        else %hear right
            belief=Posterior_ifR(belief); %belief update based on measuring RIGHT
        end
        % Recursive Bayes Filter END

        %check termination time if N is very large
        if abs(belief-(real_tiger-1))<1/5000 && endflag==false
            t_fin(Nsim)=i;
            endflag=true;
        end
        observ_plot(i,Nsim)=obs; %store observation data
        belief_plot(i,Nsim)=belief; %store posterior probability (belief) data
    end
end

% El belief (creencia que el tigre está a la derecha) tras la evidencia acumulada es:

belief

% Tiempo en estar seguro 1 contra 5000:

[min(t_fin) mean(t_fin) max(t_fin)]
stairs(0:(N-1),belief_plot,LineWidth=1.5), grid on
if Nsimulations==1
    hold on
    plot(0:(N-1),observ_plot-1,'.',MarkerSize=14)
    xline(t_fin,'-.')
    hold off
end
hold off, ylabel("creencia a posteriori p(TR|Observs.)"),xlabel("N. observaciones")

%
% Si el tigre "de verdad" está a la izquierda, entonces en media, la siguiente creencia será:

fplot([posteriori_ifHL;posteriori_ifHR],[0 1],'--',LineWidth=1,Color=[.65 .75 .75]),hold on
fplot([posteriori_ifHL*PM(hear_L,tiger_L)+posteriori_ifHR*PM(hear_R,tiger_L)],[0 1],LineWidth=4), grid on
hold on, fplot(b,[0 1],'-.',LineWidth=1), hold off
xlabel('b (Prob. right)'), ylabel('b+'), title("Average next belief if tiger is LEFT")

% Si el tigre "de verdad" está a la derecha, entonces en media, la siguiente creencia será:

fplot([posteriori_ifHL;posteriori_ifHR],[0 1],'--',LineWidth=1,Color=[.65 .75 .75]),hold on
fplot([posteriori_ifHL*PM(hear_L,tiger_R)+posteriori_ifHR*PM(hear_R,tiger_R)],[0 1],LineWidth=3), grid on
 fplot(b,[0 1],'-.',LineWidth=1), hold off
xlabel('b (Prob. right)'), ylabel('b+'), title("Average next belief if tiger is RIGHT")

%
%%
%%