function [sampleout,distout] = prob_convolve( sample1,dist1,sample2,dist2 )
%PROB_CONVOLVE gets the dist fcn and sample space of the sum of two exps.

%   sample1=the first experiment's sample space

%   dist1=the first experiments's distribution. It must be specified so
%        that dist1(i) is the probability of sample1(i).

%   sample2=the second experiment's sample space

%   dist2=the second experiments's distribution. It must be specified so
%        that dist2(i) is the probability of sample2(i).

%The actual code here mimics the convolution of two random variables as
%described at the start of Chapter 7. This has not yet been formally
%discussed in class. If you want to understand the math behind this, you
%will need to read the first 2 pages of section 7.1 in our text. Its not
%hard, but as the name suggests ("convolve") it's a bit tedious.

%Gets the size of the two sample spaces.
n1=length(sample1);
n2=length(sample2);

%Makes a place for us to examine all possible sums of sample1 and sample2.
sums=zeros(n1,n2);

%Calculates all possible sums of the outcomes in the two experiments. There
%very well may be duplicate entries. Indeed, if sample1=sample2, this is a
%symetric matrix.
for i=1:n1
    for j=1:n2
        sums(i,j)=sample1(i)+sample2(j);
    end
end

%finds all unique entries in the sums matrix. This is of course the sample
%space of the sum of the experiments.
sampleout=unique(sums);

%initializes a place for the distribution of this new sample space.
distout=zeros(size(sampleout));

%this loop figures out the distribution. The math mimicks convolution in
%section 7.1 of the text. 

%Basically we run through all the outcomes of the sum
for i=1:length(sampleout)
    
    %for each outcome of the sum, find all possible places it occurs in
    %sums. These coorespond to times when an outcome in sample1 and and
    %outcome in sample2 sum to an outcome in sampleout
    [rowinds,colinds]=find(sums==sampleout(i));
    
    %we need to sum the products of the associated entries in the
    %distribution found in the above code. We initialize a variable to hold
    %the information needed for this sum.
    temp=zeros(size(rowinds));
    
    %we now get the information needed for this sum. This is EXACTLY
    %mimicking the convolution process described in 7.1
    for j=1:length(rowinds)
        temp(j)=dist1(rowinds(j))*dist2(colinds(j));
    end
    
    %sums the distribution products found above. When this loop finishes
    %the program is done.
    distout(i)=sum(temp);
end
        

end

