%% Estimación de Variable Aleatoria
% Copyright (c) 2023 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
% Métodos Matemáticos para Telecomunicaciones.
%
% En este ejemplo se exploran diferentes relaciones estadísticas que pueden aparecer entre dos variables aleatorias y cómo dichas relaciones determinan la forma de estimar el valor de una variable aleatoria a partir del conocimiento del valor de la otra.
%% Generación de las realizaciones de la v.a. bidimensional
% Puede modificar el tipo de relación estadística entre las dos variables aleatorias que se generarán sin más que modificar el segundo argumento de la llamada a la función genera_va_bidimensional. Si desea explorar otros ejemplos y/o relaciones puede modificar el código de dicha función al final del Live Script.

N = 1000000;     % Número de realizaciones
% [x,y] = genera_va_bidimensional(N,'independent');   % Vectores con N realizaciones de la v.a. bidimensional
% [x,y] = genera_va_bidimensional(N,'correlated');   % Vectores con N realizaciones de la v.a. bidimensional
[x,y] = genera_va_bidimensional(N,'cuadratic');   % Vectores con N realizaciones de la v.a. bidimensional
% [x,y] = genera_va_bidimensional(N,'cubic');   % Vectores con N realizaciones de la v.a. bidimensional

%% Estimación de la Función Densidad de Probabilidad Conjunta
% Se estima la fdp a partir del histograma bidimensional

histogram2(x,y,'Normalization','pdf');
xlabel('x');ylabel('y');zlabel('estima f_X_Y(x,y)');title('Función Densidad de Probabilidad Conjunta')

% Observe en la figura inferior (fdp conjunta representada como mapa de calor) cómo para un valor concreto dado X=x en el eje horizontal, se puede estimar visualmente cúal es el valor de Y esperado, es decir, evaluar aproximadamente E[Y|X=x].

subplot(111);h=histogram2(x,y,'Normalization','probability','DisplayStyle','tile');
xlabel('x');ylabel('y');title('Función Densidad de Probabilidad Conjunta');colorbar;axis equal

%% Estimación de la Y a partir de X
% Se estima la variable aleatoria Y a partir del valor de X mediante el estimador sin restricciones de mínimo MSE, es decir,\hat{Y}=g(X)=E[Y|X=x].
% Para realizar esta estima de forma empírica (a partir de los datos) procedemos de la siguiente forma: se ordenan de menor a mayor según el eje x las N parejas de valores (x,y), se agrupan en bloques (clusters) de tamaño ClusterSize y se calcula el centroide (la media) de los valores x ey de cada cluster. En el límite, cuando tanto N como el número de clusters tiendan a infinito, la función representada tenderá a E[Y|X=x].

ClusterSize = 1000;
[x_ordered,indexes] = sort(x);y_ordered=y(indexes);
Xvalues = mean(reshape(x_ordered,ClusterSize,N/ClusterSize));
meanYgivenX = mean(reshape(y_ordered,ClusterSize,N/ClusterSize));

subplot(111);plot(x,y,'.')
xlabel('x');ylabel('y');title('Scatter plot');axis equal
hold on;plot(Xvalues,meanYgivenX,'r');hold off
legend('Data points','Estimated E[Y|X=x]')

%
%% Función matlab genera_va_bidimensional
% Modifique el código de esta función para generar la variable aleatoria bidimensional deseada.

function [x,y] = genera_va_bidimensional(N,tipo_relacion)
% Función de generación de variables aleatorias bidimensionales
% con cierta relación entre ellas 
% Puede modificar el código para generar la variable aleatoria bidimensional deseada
    
    switch(lower(tipo_relacion))
        case {'independent','independientes'}
            z = 4*randn(1,N)+2;w = randn(1,N)+5;x = z;y = w;  % v.a.'s independientes
        case {'correlated','correladas'}
            z = 4*randn(1,N)+2;w = randn(1,N);x = 4*z+w+2;y = 2*z-2*w;  % v.a.'s correladas
        case {'cuadratic','cuadratica','cuadrática'}
            z = 4*randn(1,N)+2;w = randn(1,N);x = 4*z+w+2;y = 0.2*z.^2-2*w;  % v.a.'s con relación cuadrática
        case {'cubic','cubica','cúbica'}
            z = 4*randn(1,N)+2;w = randn(1,N);x = -4*z-2*w+2;y = 0.02*(z-2).^3-0.2*(z).^2-2*w;  % v.a.'s con relación cúbica
        case {'problema2'}
            y = rand(1,N); g = rand(1,N);x = g.*y;            
        otherwise
            z = 4*randn(1,N)+2;w = randn(1,N);x = z;y = w;  % v.a.'s independientes
    end
end