function p=smoke(x,y)

% x= # cigarettes per day
% y= angina severity (y=0 no angina, y=1 mild angina, etc)

if y==0
   if x==0
      p=0.7;
   else
      lambda=10;
      z=exp(lambda);
      p=.3*lambda^x/(prod(1:x)*(z-1));
   end
elseif y==1
   if x==0
      p=0.6;
   else
      lambda=15;
      z=exp(lambda);
      p=.4*lambda^x/(prod(1:x)*(z-1));
   end
elseif y==2
   if x==0
      p=0.5;
   else
      lambda=20;
      z=exp(lambda);
      p=.5*lambda^x/(prod(1:x)*(z-1));
   end
elseif y==3
   if x==0
      p=0.4;
   else
      lambda=25;
      z=exp(lambda);
      p=.6*lambda^x/(prod(1:x)*(z-1));
   end
elseif y==4
   if x==0
      p=0.3;
   else
      lambda=30;
      z=exp(lambda);
      p=.7*lambda^x/(prod(1:x)*(z-1));
   end
end
