% PROGRAM NAME: twoquestions.m
% This code simulates answers to two questions
% The answer to each is 1(SD) to 5(SA)
%=============================================
%   DATA GENERATION CODE
%---------------------------------------------
p = 0.2*ones(1,5); % probabilities re: question #1
% Matirx of cond. probs for #2 (col) given #1(row)
P = [.60 .25 .11 .03 .01;
     .25 .45 .25 .03 .02;
     .10 .25 .30 .25 .10;
     .02 .03 .25 .45 .25;
     .01 .03 .11 .25 .60];
 Psum = zeros(5,5);
 Psum(:,1) = P(:,1);
 for k = 2:5
     Psum(:,k) = Psum(:,k-1) + P(:,k); 
 end
 msize = 100; % Sample size
 nsim = 1; %Number of simulations
 x = 5*ones(msize,nsim);
 y = 5*ones(msize,nsim);
 for n = 1:nsim
     for m = 1:msize
         ux = rand(1,1);
         if ux < .2
             x(m,n)=1;
         elseif ux < .4
             x(m,n) = 2;
         elseif ux < .6
             x(m,n)= 3;
         elseif ux < .8
             x(m,n) = 4;
         end
     end
     for m = 1:msize
         xval = x(m,n);
         uy = rand(1,1);
         if uy < Psum(xval,1)
             y(m,n) = 1;
         elseif uy < Psum(xval,2)
             y(m,n) = 2;
         elseif uy < Psum(xval,3)
             y(m,n) = 3;
         elseif uy < Psum(xval,4)
             y(m,n) = 4;
         end
     end
 end
 %===============================================
 %   PROBABILITY ESTIMATION CODE
 %-----------------------------------------------
 % Joint Probabilities
 Pxyhat = zeros(5,5,nsim);
 for n = 1:nsim
     for m = 1:msize
         xval = x(m,n);
         yval = y(m,n);
         Pxyhat(xval,yval,n) = Pxyhat(xval,yval,n) +1;
     end
    Pxyhat(:,:,n) = Pxyhat(:,:,n)/msize;
    % Marginal Probabilities
    pxhat = sum(Pxyhat,2);
    pyhat = sum(Pxyhat,1);
    % Conditional Probabilities of Y given X
    Phat = zeros(5,5);
    for k = 1:5
     Phat(k,:) = Pxyhat(k,:)./pxhat(k);
    end
 end
 %=================================================
 P
 Phat
 
 
         
