<< Graphics`Graphics`
<< Statistics`NormalDistribution`

Here we create the data set for the weight of 500 lab rats:

data = RandomArray[NormalDistribution[300  Pi, 3], 500];

And now we will plot the data:

Histogram[data, Frame -> True, 
      FrameLabel -> {"Weight of Rat", "Frequency"}]

We can easily change the bin size:

Histogram[data, HistogramCategories -> 50, Frame -> True, 
      FrameLabel -> {"Weight of Rat", "Frequency"}]

Histogram[data, HistogramCategories -> 5000, Frame -> True, 
      FrameLabel -> {"Weight of Rat", "Frequency"}]

Histogram[data, HistogramCategories -> 3, Frame -> True, 
      FrameLabel -> {"Weight of Rat", "Frequency"}]

Histogram[data, HistogramCategories -> 2, Frame -> True, 
      FrameLabel -> {"Weight of Rat", "Frequency"}]

Histogram[data, HistogramCategories -> 20, 
      HistogramScale -> 1, Frame -> True, 
      FrameLabel -> {"Weight of Rat", "Probability"}]

Discrete Probability

Here is the data for my coin toss game.

numclass = 4;
num = 10;
data2  =  {3, 5, 2, 0};

\(Histogram[data2, FrequencyData -> True, Frame -> True, 
      FrameLabel -> {"Number of Heads", "Frequency"}, 
      FrameTicks -> {{{0.5, "0"}, {1.5, "1"}, {2.5, "2"}, {3.5, "3"}}, Automatic, {}, {}}]\)

\(Histogram[data2, FrequencyData -> True, HistogramScale -> 1, 
      Frame -> True, 
      FrameLabel -> {"Number of Heads", "Probability"}, 
      FrameTicks -> {{{0.5, "0"}, {1.5, "1"}, {2.5, "2"}, {3.5, "3"}}, Automatic, {}, {}}]\)

Here is the check that what we have is a probability distribution:

Sum[Part[data2, i]/num, {i, 1, numclass}]

The expected value:

Sum[(i - 1)*Part[data2, i]/num, {i, 1, 
          numclass}];
EX =  N[%]

The variance:

Sum[((i - 1) - EX)^2*Part[data2, i]/num, {i, 1, 
          numclass}];
VX = N[%]

We can have Mathematica run many trials for us using the random number generator to simulate tossing the coin. We will assume that 0 represents a head.

data3  =  {0, 0, 0, 0};
num  = 100000;

Do[
{numheads = 
          0, 

Do[
{toss  = Random[Integer], 
If[toss == 0, 
              numheads  = numheads  +  1]}, {j, 1, 
            3}], 

data3[[numheads + 
                1]]  =  data3[[numheads + 1]] + 1}, {i, 1, 
        num}]

data3

\(Histogram[data3, FrequencyData -> True, HistogramScale -> 1, 
      Frame -> True, 
      FrameLabel -> {"Number of Heads", "Probability"}, 
      FrameTicks -> {{{0.5, "0"}, {1.5, "1"}, {2.5, "2"}, {3.5, "3"}}, 
          Automatic, {}, {}}]
\)
\(EX\  = \ 
      N[Sum[\((i - 1)\)*Part[data3, i]/num, {i, 1, 
            Length[data3]}]]
\)
VX = 
      N[Sum[((i - 1) - EX)^2*Part[data3, i]/num, {i, 1, 
            Length[data3]}]]

Here is the analysis changed to to four tosses and count the number of heads:

data3  =  {0, 0, 0, 0, 0};
num  = 100000;

Do[
{numheads = 
          0, 

Do[
{toss  = Random[Integer], 
If[toss == 0, 
              numheads  = numheads  +  1]}, {j, 1, 
            4}], 

data3[[numheads + 
                1]]  =  data3[[numheads + 1]] + 1}, {i, 1, 
        num}]

data3

\(Histogram[data3, FrequencyData -> True, HistogramScale -> 1, 
      Frame -> True, 
      FrameLabel -> {"Number of Heads", "Probability"}, 
      FrameTicks -> {{{0.5, "0"}, {1.5, "1"}, {2.5, "2"}, {3.5, "3"}, {4.5, "4"}}, 
          Automatic, {}, {}}]
\)
\(EX\  = \ 
      N[Sum[\((i - 1)\)*Part[data3, i]/num, {i, 1, 
            Length[data3]}]]
\)
VX = 
      N[Sum[((i - 1) - EX)^2*Part[data3, i]/num, {i, 1, 
            Length[data3]}]]

Continuous Probability

f[x_]  =  4 x^2 Exp[-2 x]
Plot[f[x], {x, 0, 5}, Frame -> True, 
      FrameLabel -> {"x", "Probability"}]

Integrate[f[x], {x, 0, Infinity}]

EX  =  Integrate[x f[x], {x, 0, Infinity}]

VX  =  Integrate[(x - EX)^2 f[x], {x, 0, Infinity}]

Distributions

a  =  2
b  =  3
f[x_]  =  If[x > b, 0, UnitStep[x - a]]
Plot[f[x], {x, 0, 5}, Frame -> True, 
      FrameLabel -> {"x", "Uniform Distribution"}]

<<Statistics`NormalDistribution`

ndist  = NormalDistribution[12, Pi]
f[x_]  =  PDF[ndist, x]
Plot[f[x], {x, 0, 24}, Frame -> True, 
      FrameLabel -> {"x", "Normal Distribution"}]

<< Statistics`DiscreteDistributions`

bdist = BinomialDistribution[34, 0.3]
f[k_]  =  PDF[bdist, k]
list  =  Table[{k, f[k]}, {k, 0, 34}]
ListPlot[list, Frame -> True, 
      FrameLabel -> {"k", "Binomial Distribution"}]