Dice and Distributions:

Introduction and Initialization:

We'll start by initializing the kernel.  We'll do that the standard way, by  asking Mathematica to solve a trivial arithmetic problem.  Click your cursor on the text below and type shift-return to execute the comand.

2+2

Die rolling and Random Numbers:

Die rolling is tedious, as you learned in lab.  We'd like to roll hundreds and thousands of dice,  so we will call on our friend the computer.  First, you need  to know how to generate a random number.  The function for that is (oddly  enough) Random.  The command below will generate a random integer between  1 and 6, inclusive:

RandomInteger[{1,6}]

We'd like a bunch of random numbers, making a long list or table.  In  Mathematica we can do this with the table command:  Look carefully at the  following command, and then execute it.

RandomInteger[{1,6},12]

You'll notice that this is equivalent to rolling the dice twelve times.  We now have  a list of 12 elements.  You will learn a lot about lists in Mathematica.  They are basic  to the way the program works, and often a pain to manipulate.    One nice trick  is to give your list a name, as we do below:

fred=RandomInteger[{1,6},12]

This both executes the command and simultaneously assigns the name ``fred''  to this list of elements.  (Always use lower case letters in the names of variables.   Mathematica assumes that upper case letters are used for functions.)

(1)  Write a command that creates a list of 10 random integers ranging from 
	        1 to 20, and assigns this list the name "foo".

Fluctuations:

We'd like to do some analysis on these random numbers, in a fashion similar to  what we did with the dice.  One handy command will be BinCounts, which is  defined in the DataManipulation package we loaded above.  First, let's see how  Mathematica defines the command.   We can get some simple information on any Mathematica command by typing ?Command.

(2)  Execute the command ?BinCounts and read its description.

This will be handy for sorting the dice.  The command below will sort the list  "fred" that we defined above, and count how many 1's, 2's, 3's, etc. we have.   Make sure you understand both the syntax of the command, and the output.

BinCounts[fred,{0,6,1}]

(3) Perform the same operation on the data set you named "foo".  How many bins
	     should you have?

Of course we want to roll more that 12 dice!  We'd like to roll several hundred.  The command below makes a list (table) of 120 random numbers between 1 and 6,  and then bins them in integer bins from 1 to 6..  Execute it.

BinCounts[RandomInteger[{1,6},120],{0,6,1}]

This distribution is P(n), which you calculated in lab.  You can get a histogram of  if with the following command:

BarChart[BinCounts[RandomInteger[{1,6},120],{0,6,1}]]

One would expect that each number occurs 1/6 of the time, or 20 times out of 120.  You'll notice that each number does not have the same number of counts - some  occur  more often, some less.   From your work in lab rolling dice, you should have  expected this.   There will be fluctuations about the average value.   Let's see if  we can be quantitative.
	
	First let's generate a set of numbers, bin them, and then assign a name to the list  of hits.

betsy=BinCounts[RandomInteger[{1,6},120],{0,6,1}]

Now, you could caluclate the mean  and standard deviation as you did when you  rolled the dice.  But there's a simple way in Mathematica.  (These commands were  loaded in the package DescriptiveStatistics.)

Mean[betsy]

StandardDeviation[betsy]

Oops!    Let's force Mathematica to turn that into a decimal number.

N[StandardDeviation[betsy]]

Ok, that's better.   We now now the average number of times each number appears,  and the fluctuations.   Now wasn't that much easier than doing it by hand?  The   relative fluctuation is the standard deviation/mean.

N[((StandardDeviation[betsy])/(Mean[betsy]))]

However, we know that for  six sided dice, the mean is just n/6.  

Below is a command that will calculate the standard deviation of (10^n) die rolls, for  n between 1 and 4.  Note that the inner Table command is the one in which the random numbers are generated.  The ``outer''  Table command loops over the different values of n.  We can then study how the standard deviation increases as we increase the number of dice rolled, or random numbers generated.

mike=Table[StandardDeviation[BinCounts[RandomInteger[{1,6},(10^n)],{0,6,1}]],{n,1,4}]

The ListPlot command allows us to plot this set of numbers.

ListPlot[mike,Joined->True]

(4) Copy the above command and r epeat the above analysis for 12, 48,  96, 192, 384, 768, 1536, 2072, 6144, and  12288 rolls of 12 sided  ``dice'' calculating the relative fluctuations as a function of the  number of rolls..    (These numbers are  the set 12 (2^n) for n ranging  from 0 to 10.)    Use the command ListPlot to plot the relative fluctuations  as a function of ``n''.  (You may want to execute ?ListPlot to learn about  this function.)    Can you see any systematic dependence?   Can you guess its  functional form?  You may wish to try LogListPlot or LogLogListPlot.   There is a lot of  scatter to this graph, due to the random nature of the data.   This is an exploration question;   you should try several different plots.

Partitions and Distributions

Recall the function N(n) you calculated in the lab.  This gave the number of singles,  pairs, triples, etc. you rolled.   We'll show you how to do this here, using Mathematica.
	
	Recall our list "fred".    We can run in through our binning command, as we did above,  and store the list in a string called "bob".

bob=BinCounts[fred,{0,6,1}]

Now look at this output and think about what it means.  A "2" means that a given  number appeared twice in the list of random numbers.  The number of "2's" in this  list is therefore the number of pairs in the set.  The number of "3's" is the number of triples, and so on.  We just need to count the number of times, "1" occurs in the  list to get N(1), the number of times 2 appears to get N(2), etc.  But we know how  to do this!  This is just the BinCounts command.

BinCounts[bob,{0,12,1}]

Make sure you understand the output of the above command before proceeding.  You should scroll back and look at the list "fred" and confirm that the first element really is the number of singles, the second  is the number of pairs, etc.  If you are confused,  please talk with your TA.
	
	Below is a command which you can alter to change the number of sides on each  die, the number of dice you roll at a time, and the number of total rolls you make.  You would not be expected to come up with this on your own, but you should be able  to generally understand that it does exactly what you did above.  
	
	The case below consists of 20 dice that each have 40 sides.  You roll the set of  20 dice 100 times.  The list N(n) is given the name "dist".  Execute the command.

nside=40;
ndice=20;
nroll=100;
dist=(∑_(n=1)^nroll)BinCounts[BinCounts[RandomInteger[{1,nside},ndice],{0,nside,1}],{0,ndice,1}]

Ok, that's nice, but let's plot this.   We'll use the ListPlot command from above,  and will tell Mathematica to join the points.

ListPlot[dist,Joined->True]

Ugh!  Mathematica was pretty stupid.  It didn't even pick the limits right!  We  can get a nicer plot by scaling the y-axis logarithmically.  This is easily done:

ListLogPlot[dist,Joined->True]

Mathematica complains about the zeroes in our list, since it can't take the  log of zero.  But you'll notice that we get a nice straight line.

(5)  Copy and change the above command to roll 10 20-sided dice a 
	     thousand times.  Again, make a log plot as you did above.   From this 
	     curve make a guess at the form of N(n).

You might guess that the shape of N(n) depends upon the values of ndice and  nside.  (The variable nroll determines the overall scale of the curve, but not its  true shape.)   Recall our analogy in the handout.   The variable ndice corresponds  to the number of packets of energy you have.  You assign a packet of energy to  a particle depending upon the number you roll on the die.  Thus nside corresponds  to the number of particles.

(6)  Copy and change the above command, so that ndice is twice nside.  
	    (You may fiddle with the exact values of nside, ndice, and nroll.)
	     Calculate and plot the distribution on a log scale.  Repeat this with ndice 
	     equal to four times nside.   How does these distributions compare with the 
	     one you obtained in (5)?  Can you explain the changes qualitatively?

We note in passing that the curves you generate above are for a system with a fixed amount of energy.  We usually work at a fixed temperature.   The exact  relationship between two cases goes beyond this course.  But the general  feature you observe in (5) are true:  there are far  fewer particles with a lot of  energy.

Extra Credit

You'll notice that we didn't calculate the distribution S(n).   This was the sum  of the dice you rolled.   For extra credit you may use Mathematica to calculate  S(n) for six 8-sided dice.   You'll want to use the Sum command.  First write a command to generate the six random numbers.  Then alter it to calculate the sum of them.  Next, make a table  of a large number of such sums, using the Table command as we did above. .  Finally, run that list through BinCounts.    

	How does the function S(n) look?  How does its width change as you change  the number of dice?   How about its width relative to the peak value?   Can you  guess a functional form for it?