JointDistCommands.m

%% generate the data
n = 1000;
% which die was used first and second roll - fair, hi, lo
dice1 = randi(3,[n,1]);
dice2 = randi(3,[n,1]);
% did the 'loading' happen for die 1 and die 2
load1 = randi(2,[n,1]);
load2 = randi(2,[n,1]);
r1 = roll(dice1,load1,randi(5,[n,1]),randi(6,[n,1]));
r2 = roll(dice2,load2,randi(5,[n,1]),randi(6,[n,1]));
D = [dice1,dice2,r1,r2];
%% There are lots of ways to visualize it
% show this as an image
imagesc(D);
% sort and then show it as an image 
S = sortrows(D,1:2);
imagesc(S);
% show a histogram
D34 = D(:,3:4);
hist3(D34,[6,6]);
% get the counts and 'centers' that define the histogram
[H,C] = hist3(D34,[6,6]);
P = H/1000;
SP = (H + (1/36))/1001;
[I,J]=meshgrid(1:6,1:6);
surf(I,J,SP);
% we could also 'jitter' the points and do a scatter plot
E = D34 + randn(1000,2)*0.1;
plot(E(:,1),E(:,2),'r*');
% or combine that with a surface....
surf(I,J,SP);
hold on;
plot3(E(:,1),E(:,2),zeros(1000,1)+0.1,'r*');
%% Reasoning with the joint
% easy thing is to just look up interesting cases in the data and see how
% frequent they are
sum(D(:,3)==D(:,4))/1000;   % doubles
sum(D(:,3)+D(:,4)==7)/1000; % sevens
% find all the indices of the data - in this case these are the different
% outcomes of the joint experiment
[I,J]=find(SP);  %all possible indices of SP
IJ = [I,J];      %index 1,index 2
% then we can find where - at what indices - we get sevens or doubles conveniently...
Sevens = IJ(I==J,:);
Doubles = IJ(I==J,:);
sum(SP(sub2ind(size(SP),Sevens(:,1),Sevens(:,2))));
sum(SP(sub2ind(size(SP),Doubles(:,1),Doubles(:,2))));
% or even ask a hard question like: what is Pr(both die fair|roll>10)?
fairAnd11Or12 = sum((D(:,1)==1) & (D(:,2)==1) & (D(:,3)+D(:,4) >= 10));
is11Or12 = sum((D(:,3)+D(:,4) >= 10));
fairAnd11Or12/is11Or12;