%Generates de conditional histogram for P(L2|L1)
%where L1 and L2 are (simulated) linear filter responses
clear;
binsx=25; %Resolution in L1
binsy=41; %Resolution in L2
samples=100000; %Sample pairs (number of filter outputs) for simulation
L1=rand(1,samples)-0.5; %Assuming we have enough examples to generate P(L2|L1)
L2=randn(1,samples); %Generate L2 independently
L2=sqrt(abs(L1)).*L2; %Introducing dependence on variance of L1
jhist=zeros(binsy,binsx);
minL1=abs(min(L1));
bwL1=(max(L1)-min(L1))/binsx;
minL2=abs(min(L2));
bwL2=(max(L2)-min(L2))/binsy;
for i=1:length(L1)  %Calculating joint histogram
    j=fix((L1(i)+minL1)/bwL1)+1;
    k=fix((L2(i)+minL2)/bwL2)+1;
    if j==binsx+1 j=binsx; end
    if k==binsy+1 k=binsy; end
    jhist(k,j)=jhist(k,j)+1;
end
jhist=jhist+1; %Regularizing
for i=1:binsx
%    jhist(:,i)=jhist(:,i)/sum(jhist(:,i)); %Calculating conditional histograms
    jhist(:,i)=jhist(:,i)/max(jhist(:,i)); %Adjusting for easy viewing
end
framex=0.1; %Defining frame
framey=0.3;
%Generating image:
imagesc(jhist(round(binsy*framey):binsy-round(binsy*framey),round(binsx*framex):binsx-round(binsx*framex)));colormap(gray);axis equal;axis off;
