(* Solution to Problem 4 in Lecture 8.*)

(* Here is the solution to Part a):*)

b[x_, p_,k_]:=Binomial[k,x]p^x(1-p)^(k-x)

(* Here is an illustration. The N@ causes the output to be written in deciaml form. We are letting p=1/4 and k=4. *)

Table[{x,N@b[x, 1/4, 4]}, {x, 0, 4}]

(* Here is the solution to Part b).  Here p is still 1/4, but we are letting k=100.  We are making 100 flips.  *)

Sum[N@b[x, 1/4, 100], {x,31, 40}]

(************************************************)

(* Here is what we did previously in Problem 3:*)

f[x_,n_,s_,k_]:=Binomial[s,x]Binomial[n-s, k-x]/Binomial[n,k]

(* Here is an illustration.  We are drawing 4 items without replacement from 1000 items, of which 1/4 are special. *)

Table[{x,N@f[x, 1000, 250,4]}, {x, 0, 4}]

(* Now in Part b) of Problem 3 (Lecture 7), we calculated the following.*)

Sum[N@f[x,1000, 250, 100], {x,31, 40}]

(* Since we are drawing without replacement, the draws are not independent.  But the larger the population in the urn, the less they influence one another.  Drawing from an urn with 100 million beads, of which 25 million are special is bey much like flipping a coin that lands on its special side with probability 1/4.  The following calculations demonstrate this. *)

Sum[N@f[x,10000, 2500, 100], {x,31, 40}]

Sum[N@f[x,100000, 25000, 100], {x,31, 40}]

Sum[N@f[x,1000000, 250000, 100], {x,31, 40}]

Sum[N@f[x,10000000, 2500000, 100], {x,31, 40}]

Sum[N@f[x,100000000, 25000000, 100], {x,31, 40}]

(************************************************)
(************************************************)
(************************************************)
(************************************************)

(* Lecture 8, Example 2.*)

buyabox[]:=RandomInteger[{1,6}]

(*Illustration*)

Table[buyabox[], {20}]

shoptillcomplete[]:=(
boxesbought=0;(* Sets the variable "boxesbought" equal to 0. *)
couponsowned={};(* Sets the variable "couponsowned" equal to {}. *)
Do[(* Look up "Do" in Mathematica Help. *) 
thiscoupon=buyabox[] ;(* Sets the variable "thiscoupon" equal to a randon integer from {1,2,3,4,5,6}. *)
boxesbought=boxesbought+1;(* Records the fact that we have bought another box. *)
couponsowned=Union@Append[couponsowned, thiscoupon];(* Adds thiscoupon to our list of coupons owned, if our list does not contain it already.  Look up Union and Append. *)
If[Length[couponsowned]==6, Break[]],(* Checks to see if we have six numbers in our list, couponsowned.  If we do, we break out of the Do loop. *) 
{1000}(* This limits the number of boxes we buy to 1000.  It is possible---but very improbable---that we might not have all six coupons even after buying 1000 boxes.*) 
];(* Here ends the Do loop.*)
boxesbought (* Calling the variable boxesbought now causes Mathematica to return its value when the Do loop ended.*)
)

shoptillcomplete[](*If we input this function, the output reports the number of boxes we had to purcahse in one simulation to get all six coupons.*)

mylist=Table[shoptillcomplete[], {1000}] (* We send out 1000 shoppers, each of whom buys cereal boxes until he/she has a complete set of coupons.  We get a list of the number of boxes each shopper had to buy.*)

N[Mean[mylist]](*This computes the average number of boxes purchased by the 1000.*)

N[Mean[Table[shoptillcomplete[], {10000}]]](* Let's turn 10000 shoppers loose, and then average their buys. *)

(* A mystery! Look at this!  Is it a coincidence that this number is close the the average, above? *)

N[1+6/5+6/4+6/3+6/2+6]