[X]=e446

function [X]=e446
%
%  Example 4.4.6 - the coupon collection example
%
m=52;
x=ones(1,52); x=cumsum(x); pi=(m-(x-1))/m;
qi=1-pi; Q=diag(qi); p1=pi; pi(52)=[];  A=diag(pi,1);
Q=Q+A; alpha=zeros(1,52); alpha(1,1)=1;
r=zeros(52,1); r(m,1)=p1(m); I=eye(m);
Q1=inv(I-Q); E1=sum(alpha*Q1);    %  1st moment
EF2=2*sum(alpha*Q*Q1*Q1);         %  2nd factorial moment
E2=EF2+E1;                        %  by (1.2.4)
Var=E2-E1*E1; Sdv=sqrt(Var);
fprintf('  E[T] = %12.4f  \n',E1);
fprintf('Var[T] = %12.4f  \n',Sdv);
pk=[]; S=alpha*I;
for k=1:500 
   x=S*r; S=S*Q; pk=[pk x]; 
end
X=zeros(500,2); x=ones(1,500); x=cumsum(x); X(:,1)=x';
X(:,2)=pk';