S=1000; %Number of simulation rounds
p=[.1 .4 .3 .2]; %PMF for x
P=[.1 .5 .8 1]; %CDF for x
X=[0 1 2 3]; %Possible values of x

%a)
Sample=zeros(S,1); %Initialize the sample vector

for k=1:S
    r=rand(); %Random number from Uni(0,1)
    counter=1; %Start looking from the first value of X
    while (r>P(counter)) %While r is greater than CDF at current value of X... 
        counter=counter+1; %We go to the next value of X
    end %When r is lower than the CDF at the current value of X...
    Sample(k)=X(counter); %We have found the value of X corresponding to r
end

%b)
SampleMean=mean(Sample)
SampleStd=std(Sample)

TrueMean=p*X'
TrueStd=sqrt(p*(X.^2)'-TrueMean^2)

%As an extra thing, let's calculate the ratios of different order sizes in
%the sample:
ratios=zeros(1,4);
for i=1:4
    ratios(i)=sum(Sample==i-1)/length(Sample);
end
ratios