There is a huge library for Probability and Statistics in Mathematica.  We will explore some simple subjects.  

To truly understand this subject  requires fairly advanced Mathematics.  It is not difficult, it just requires patience and time.  

As with all things, learn the basics in detail, know definitions in detail, and then build up to more complex subjects.

Factorials come into use when evaluating the number of ways to arrange objects.  
Factorial[5]= 5 ! = 5*4*3*2*1

5!

Binomial[n, m]  gives the binomial coefficient    n !/(m ! (n - m)!) . 

Suppose you have 4 boys in the class, and you must pick 2 to run across the school to get you some colored pencils and chaulk.  

Let’s call the boys “Ketevi”, “Bob”,  “Yash”, and “Johannes” 
How many ways can you get 2 boys out of 4.

listofboys = {"Ketevi","Bob","Yash","Johannes"};  
chooseboys= Subsets[listofboys,{2}]

Let’s analyze this somewhat.  There are 4 ways to pick the first boy and 3 ways to pick the second. 
That gives me 12. But   since 2 boys are being picked they can be arranged in 2 ways and therefore there will be repetition of 2. 

4*3/2  =>      4*3*2/(2*2)    =   4!/(2!)(4-2)!    = 6 

This is an easy way to understand the formula, but it is applicable to any number n and m.

And so,  This should be 6 ways or   Binomial [4,2]

Binomial[4,2]

To randomly pick from these 6 choices we can use RandomChoice which picks an element from a given list.

RandomChoice[chooseboys]

Now about  probability density functions.  
First some definitions.  
Any choice can be mapped on the real continuous or discrete, integer number.  

Choice of red,  green, white  balls can be choice between 1, and 2, and 3.  

After they are mapped, the choices can be ordered.  3 comes after 2 which is after 1 for example.  

Now probabilities can be assigned. 

P(1) = 0.1  
P(2) = 0.3 
P(3) = 0.6  

But now the sum of all must add to 1.  P(1)+P(2)+P(3) = 1  

P(n) is actually called the Probability Mass Function (PMF). anf if the variable is continuous, P(x) is called the Probability Density Function (PDF). 

The stepwise function   P(1),  P(1)+P(2), P(1)+P(2)+P(3) is called the probability distribution.   For continuous variable this is an integral ! 



The numbers (1,2,3) are called the domain of this probability distribution.  

Now, unfortunately, this definition is often switched, even in rather famous works. Just make sure you understand what the definition is every time !

Pn[n_]:= Which[n==1, 0.1, n==2, 0.3, n==3, 0.6];

pdist =  Accumulate[Table[Pn[n],{n, 1,3}]];

ListPlot[Table[{n, Pn[n]}, {n, 1, 3}] ]

BarChart[pdist]

Distributions can be specified by a name and its parameters.  Often 
in physics the PMF and PDF are called the Probability Distribution. There is not much confusion here usually, except when precise language is needed.  

BinomialDistribution[n, p] where n is the number of elements and p is the probability of choosing an element.  
The random variable k is the number of chosen elements.  

k  which is from 0 to n  forms the domain of the distribution.  

PDF[BinomialDistribution[n, p], k]  is the function of k. 
We can make a plot of this function

PDF[BinomialDistribution[n,p],k]

Table[{k,PDF[BinomialDistribution[5,1/2],k]},{k,0,5}]//TableForm

?ListPlot

ListPlot[Table[{k,PDF[BinomialDistribution[5,1/2],k]},{k,0,5}],Frame->True,FrameLabel->{"Successes","Probability"},Filling->Axis, FrameStyle->Bold,LabelStyle->Directive[Bold,Medium]]

Some Continuous probability density functions.  
Normal of Gaussian distribution is 

NormalDistribution[mean, sigma]

PDF[NormalDistribution[0,1],x]

Plot[PDF[NormalDistribution[0,1],x],{x,-5,5},Filling->Axis]

PDF[ExponentialDistribution[1],x]

Plot[PDF[ExponentialDistribution[1],x],{x,0,5},Filling->Axis]

What is the chance that none of you 20 people have the same birthday ? 
First person has some birthday out of d=365 days  
Chance that second person does not share the birth of first is:  (365-1)/365   (d-1)/d  
Chance that third person does not share with first two is:  (d-2)/d  
etc.

q1= Product[(d-i)/d,{i,1,19}]

q1/.d->365//N//ScientificForm

This should be the same as  (which is actually reasonable, but gets small quickly.  )  Let’s plot it.

Q1[n_,d_]:=d !/((d-n)!*d^n)

Q1[20,365]//N

ListPlot[Table[Q1[n,365],{n, 2, 100}], Frame -> True,  FrameStyle -> Bold, FrameLabel->{"Number of Friends", "Chance of Not Sharing B'day"}, LabelStyle -> {Large, Bold, Italic}]

Notice that the complement of Q1 is 1-Q1, it is the chance that one or more share a birthday

What is the chance that out of 20 people 2 or more share a birthday. It is amazingly large.  Nearly 41 %.

p1 = 1-  Q1[20,365] //N

Statistics is the subject of extracting information about the parameters of the underlying probability distribution from  data. 

We can put data into Mathematica and do simple and very advanced analysis on it.  

Let’s first create some data.

mydata = {2,1,3,5,3.2,1.1,5.6,2,1.,1.1,3.6,9.4,5.,5.1,5.4,4.9};

Mean[mydata]

StandardDeviation[mydata]

ListPlot[mydata] (* makes a plot of the data versus its index number *)

We can use Histogram function to plot the frequency of data.  How many in each bin.  
Histogram[data, bins]

?Histogram

Histogram[mydata,20]

How about pairs of data points {x,y}

pairdata =  {{-2.,-1.9},{-1.,-0.9},{0.,0.2},{1,1.01},{2., 2.5},{3,2.6},{5,5.9}}

dataplot = ListPlot[pairdata,Frame->True,FrameLabel->{"X","Y"},FrameStyle->Bold,LabelStyle->Large]

How about making a linear fit to this data ?

?Fit

myfit =  Fit[pairdata,{1,x}, x]

Show[dataplot,  Plot[myfit,{x,-2,5}]]

Let’s talk a little bit about simulations.  Often Monte-Carlo simulations are used for generating random data that can be tested.  

How do we generate a sequence of random numbers that mimics a probability density function ? 

Take the example of an Exponential function.  

P(x) =  ({{{, ({{(e^(-x)), x≥0}, {0, elsewise}})}})

Integrate of this is the CDF  (Cumulate Density Function.  Also called the “Distribution Function” )  

CDF[x] =  ({{{, ({{1-(e^(-x)), x≥0}, {0, elsewise}})}}) 

Now if I randomly throw a number between 0 and 1 which is the range of the PDF and calculate the corresponding x, I should have a properly sampled random number ~ P(x)  that corresponds to P(x) as the PDF.

CDF[ExponentialDistribution[1],x]

Plot[{PDF[ExponentialDistribution[1],x],CDF[ExponentialDistribution[1],x] },{x,0,5},Filling->Axis]

px =  RandomReal[]  (* gives a random number between 0 and 1 *)

Xcall = -Log[1-px]  (* this gives the answer *)

dataexp = Table[-Log[1-RandomReal[]], {n, 1, 100}]   (* Let;s make 100 of these *)

Histogram[dataexp]

All computer programs have inbuilt random number generations.  Mathematica gives built in way to generate any of the PDFs that are in-built.

?RandomVariate

makenormal = RandomVariate[NormalDistribution[0,1],10000] ; 
Histogram[makenormal]

This final lesson may be the most important lesson of  all.  

This is called the Central Limit Theorem.  It is at the heart of all measurements and all of your daily work.  

When you observe some natural occurrence,  your mind naturally attempt to take an average. 
This average is a very good thing to rely on for the future.  It is always distributed around some mean value with a Normal probability density function  !  

The observation needs to have some mean, and needs to have some sort of tendency to be around that mean, but it could be distributed in any way.    

The distribution of the means must be Normal !  

Let’s do this exercise.  let’s make 1000 means of a simple number thrown 5 times each.  

RandomReal[{-1,1},5]  produces a number anywhere from {-1,1},  5 times.   
I take the mean and produce 1 million events.  They should be distributed in a normal manner. 

Play with this idea and try to see where the rule breaks down.  And what is actually applicable.

mydata = Table[ Mean[RandomReal[{-1,1},10]] ,{i, 1, 1000000}] ; 

c3= Histogram[mydata,{0.05}] ;  (* histogram with bin size of 0.05 *)  

c2= Plot[1000000*0.05*PDF[NormalDistribution[0,2/Sqrt[12]/Sqrt[10]],x],{x,-1.,1.}];

Show[c2,c3]   (* remember that histograms have a bin width and so this needs to be adjusted when making the plot *)

Mean[mydata]

2/Sqrt[12]/Sqrt[20]//N

StandardDeviation[mydata]

?RandomReal

The best way to learn  is by trying simple things first. Get a deeper understanding by using variations on some simple problem for which you already know the answer.  

Use the Mathematica website to learn about functions. Most simple problems have been solved by someone and the solution is available on the web.  

Complicated problems can be pieced together by breaking them up in simple ones.  Setup each simple one and make sure it works, and then move on to put it all together.