% GEND file gend.m   Marginal and joint dbn for integer compound demand
% Version of 5/21/97
% Calculates marginal distribution for compound demand D
% and joint distribution for {N,D} in the integer case
% Do not forget zero coefficients for missing powers
% in the generating functions for N, Y
disp('Do not forget zero coefficients for missing powers')
gn = input('Enter gen fn COEFFICIENTS for gN  ');
gy = input('Enter gen fn COEFFICIENTS for gY  ');
n  = length(gn) - 1;           % Highest power in gN
m  = length(gy) - 1;           % Highest power in gY
P  = zeros(n + 1,n*m + 1);     % Base for generating P
y  = 1;                        % Initialization
P(1,1) = gn(1);                % First row of P (P(N=0) in the first position)
for i = 1:n                    % Row by row determination of P
   y  = conv(y,gy);            % Successive powers of gy
   P(i+1,1:i*m+1) = y*gn(i+1); % Successive rows of P
end
PD = sum(P);                   % Probability for each possible value of D
a  = find(gn);                 % Location of nonzero N probabilities
b  = find(PD);                 % Location of nonzero D probabilities
P  = P(a,b);                   % Removal of zero rows and columns
P  = rot90(P);                 % Orientation as on the plane
N  = 0:n;
N  = N(a);                     % N values with positive probabilites
PN = gn(a);                    % Positive N probabilities
Y  = 0:m;                      % All possible values of Y
Y  = Y(find(gy));              % Y values with positive probabilities
PY = gy(find(gy));             % Positive Y proabilities
D  = 0:n*m;                    % All possible values of D
PD = PD(b);                    % Positive D probabilities
D  = D(b);                     % D values with positive probabilities
gD = [D; PD]';                 % Display combination
disp('Results are in N, PN, Y, PY, D, PD, P')
disp('May use jcalc or jcalcf on N, D, P')
disp('To view distribution for D, call for gD')