({{, FCM1013: Discrete Mathematics}, {, Propotional Logic &Probability with Mathematica}})

March 16, 2016

References
[1] K.H. Rosen, Discrete Mathematics and Its Applications, 7th ed., : McGraw-Hill, 2012.

2 - Probability

Introduction to Dicrete Probability

Bayes’ Theorem

Expected Value and Variance

2.1 Introduction to Discrete Probability

2.1.1  Discrete Probability

Summary

When we think of probability we think of chance, such as the chance of flipping two heads in a row from two coin tosses, or the chance that a face card will come up in a deck of cards.  
The probability is a number that defines the “chance”.  We call all the possible outcomes such as HH, HT, TH, TT for flipping two coins the sample space.  
Each of these four are called outcomes, and any individual one or collection of individual ones is called and event.

Binomial

To find the probability of an event in a finite sample space, you calculate the number of times the event occurs and divide by the total number of possible outcomes (the size of the sample space).

Binomial

Example

We calculate the probability of winning a lottery, where we need to choose 6 numbers correctly out of 40 possible numbers.
The total number of ways to choose 6 numbers is:

Binomial[40,6]

Since there is one winning combination, the probability is 1 divided by that value.

1/Binomial[40,6]

create a definition

Continuing with this type of example, we define a function that computes the probability of winning a lottery where 6 numbers must be matched out of n possible numbers, provided that n is at least 6.

lottery[n_Integer]/;n≥6:=1.0/Binomial[n,6]

Then the probabilities above can be computed with the function.

{lottery[40],lottery[50]}

Clear[lottery]

Visualize the data

Now we can use the Table function to look at a list of probabilities for a range of values of n and graph these values with ListPlot to visualize how the number of possible values in the lottery affects the probability of choosing the correct values.

Table

ListPlot

lotteryValues=Table[{n,lottery[n]},{n,40,50}]

TableForm[lotteryValues]

ListPlot[lotteryValues]

D.I.Y

Exercise 1

Calculate the probability of winning a lottery, where we need to choose 5 numbers correctly out of 20 possible numbers.
The total number of ways to choose 5 numbers is:

winLottery[n_Integer]/;n≥5:=1.0/Binomial[n,5]

(*Tabulate data*)
TableForm[Table[{n,winLottery[n]},{n,10,20}],TableHeadings->{None, {"Possible numbers","Probability"}}]

(*Visualize data*)
ListPlot[Table[{n,winLottery[n]},{n,10,20}],AxesLabel->{"n","Probability"}]

2.1.2  Probability Distributions

Summary

In order to compute probabilities of events in terms of random variables, we must specify the distribution of the random variable.

Computing Probabilities from Discrete Distributions

Mathematica provides many different probability distributions, including several discrete distributions. To compute the probability of an event involving a random variable distributed according to a built-in distribution, you can use the Probability or NProbability functions. These functions are essentially identical, except NProbability applies numerical methods to approximate a numerical answer, while Probability uses symbolic methods. Using either function requires that we have a distribution, so we begin by introducing Mathematica’s implementation of the binomial distribution.

Probability Distributions - BinomialDistribution

The BinomialDistribution takes two parameters: the number of independent Bernoulli trials and the probability of a success.

For example, to calculate the probability that eight of the bits in a string of ten bits are 0s if the probability of a 0 bit is 0.9, the probability of a 1 bit is 0.1, and the bits are generated independently. We interpret success to be a 0 bit, so the probability of success is 0.9. The binomial distribution with these parameters would be represented as shown below. Note that BinomialDistribution is not a function and evaluating the line below will do nothing.

X\[Distributed]BinomialDistribution[10,0.9]

This is done using the Distributed (\[Distributed]) function or operator. As a function, Distributed requires two arguments: the symbol representing the random variable and the representation of the distribution. To use the operator form, enter EscdistEsc to produce the \[Distributed] symbol.

Example

The expression below computes the probability of 8 successes for a binomial distribution with 10 trials and probability of success 0.9.

Probability[X==8,X\[Distributed]BinomialDistribution[10,0.9]]

Probability Distributions - BernoulliDistribution

Related to the binomial distribution is the probability distribution of a single Bernoulli trial. The BernoulliDistribution takes only one parameter, the probability of success. The following creates the random variable associated to a single trial with probability of success 0.9. The possible values of a Bernoulli random variable are 0 or 1, representing failure and success, respectively.

Probability[X==0,X\[Distributed]BernoulliDistribution[.9]]

Probability Distributions - GeometricDistribution

In Mathematica, the value of a random variable distributed geometrically is the number of failures that occur before the first success. The GeometricDistribution in Mathematica requires one argument, the probability of success. 
The following computes the probability that if a coin with probability of heads 0.3 is flipped until it comes up heads, there will be at most 5 tails that appear.

Probability[X≤5,X\[Distributed]GeometricDistribution[.3]]

Use Boolean operators when building the event. For example, the probability that it requires at most 2 or at least 7 flips of a coin with probability of heads 0.3 before a heads comes up is:

Probability[X≤2||X≥7,X\[Distributed]GeometricDistribution[.3]]

Probability Distributions - Multiple Distributions

Suppose you have two coins, one with probability of heads 0.3 and the other with probability of heads 0.25. The probability that flipping each coin once will result in at least one head can be computed as follows.

Probability[X==1||Y==1,{X\[Distributed]BernoulliDistribution[.3],Y\[Distributed]BernoulliDistribution[.25]}]

Note that the second argument is a list whose elements define the distributions for X and Y.

2.2 Bayes’ Theorem

2.2.1  Introduction

(A) From Wolfram MathWorld ( http://mathworld.wolfram.com - the web’s most extensive mathematics resource )

(B) From Wolfram|Alpha

WolframAlpha["bayes theorem",{{"FormulasPod:FamousMathProblem",1},"FormulaData"},PodStates->{"HistoryPod:FamousMathProblem__More"}]

ReleaseHold[Hold[P[(Subscript[A, i])|A]==((P[(Subscript[A, i])] P[A|(Subscript[A, i])])/((∑_(j=1)^N)P[(Subscript[A, j])] P[A|(Subscript[A, j])]))]]

P[(Subscript[A, i])|A]==((P[A|(Subscript[A, i])] P[(Subscript[A, i])])/((∑_(j=1)^N)P[A|(Subscript[A, j])] P[(Subscript[A, j])]))//TraditionalForm

2.2.2  Applying BAYES’ Theorem

Bayes' theorem can be used to solve problems that arise in many disciplines. We will discuss an application of Bayes' theorem to medicine. In particular, we will illustrate how Bayes' theorem can be used to assess the probability that someone testing positive for a disease actually has this disease. The results obtained from Bayes' theorem are often somewhat surprising, as Example shows.

Example

Example (From K. H. Rosen: Discrete Mathematics and Its Applications 7th Edition: pg 471)

Suppose that one person in 100,000 has a particular rare disease for which there is a fairly accurate diagnostic test. This test is correct 99.0% of the time when given to a person selected at random who has the disease; it is correct 99.5% of the time when given to a person selected at random who does not have the disease. Given this information can we find 

(a) the probability that a person who tests POSITIVE for the disease has the disease?
(b) the probability that a person who tests NEGATIVE for the disease does not have the disease?

Should a person who tests positive be very concerned that he or she has the disease?

Solution:

F be the event that a person selected at random has the disease;
p(F) = 1/100,000 = 0.00001

Thus, the event that a person selected at random has no disease;
p((OverBar[F])) = 1 - p(F)    = 1 - 0.00001 = 0.99999

A person who has the disease tests positive 99% of the time;
p(E | F) = 0.99

Thus, a person who has the disease tests negative;
p((OverBar[E]) | F) = 1- p(E | F) = 0.01

a person who does not have the disease tests negative 99.5% of the time;
p((OverBar[E]) | (OverBar[F])) = 0.995

a person without the disease tests positive
p(E | (OverBar[F])) = 1- p((OverBar[E]) | (OverBar[F])) = 0.005

Use Bayes’ theorem to compute p(F | E).

(a) the probability that a person who tests positive for the disease has the disease = p(F |E)

({{({{p(F|E)=((p(E|F)p(F))/(p(E|F)p(F)+p(E|(OverBar[F]))p((OverBar[F]))))}, {= ((0.99 0.00001)/(0.99 0.00001+0.005 0.99999))}, {= 0.001976107069073924`}})}})

(b) the probability that a person who tests negative for the disease does not have the disease

({{({{p((OverBar[F])|(OverBar[E]))=((p((OverBar[E])|(OverBar[F]))p((OverBar[F])))/(p((OverBar[E])|(OverBar[F]))p((OverBar[F]))+p((OverBar[E])|F)p(F)))}, {=((0.995 0.99999)/(0.995 0.99999+0.01 0.00001))}, {≈0.9999999}})}})

For rare diseases, people tend to intuitively overestimate the probability of being sick after having received a positive test result. This probability is calculated using Bayes' s theorem from some of the conditional probabilities involved in the scenario. 

Understanding these probabilities, such as the probability of being sick when having received a positive test result, expressed as P (sick | positive test), can become easier using a contingency table like the one shown. 

Given properties of the test, the probability of the disease, and the size of the reference group, the table shows the number of people falling into the four different possible categories. The conditional probability of being sick after having received a positive test result is simply the ratio of the people that are sick and that have tested positive to the total number of people having tested positive.

Applying Bayes’ Theorem in Computable Document Format (CDF)

ProbabilityOfBeingSick[prior_,recogdisease_,recoghealth_,ntot_]:=Module[{res,nsick,nhealthy,rightpositive,rightnegative,
wrongpositive,wrongnegative, a11,a12,a13,a21,a22,a23,outstring},

a13=nsick= Round[ntot prior];
a23=nhealthy= Round[ntot- nsick];
a11=rightpositive= Round[nsick recogdisease ]; 
a12=wrongnegative= Round[nsick - rightpositive];
a22 =rightnegative = Round[nhealthy recoghealth]; 
a21=wrongpositive= Round[nhealthy - rightnegative];
If[a11+a21>0,outstring = "P(sick|positive test): "<>ToString[a11]<>"/"<>ToString[a11+a21]<>" = "<>ToString[N[a11/(a11+a21)]],
outstring = "P(sick|positive test): 0"
];

res=Grid[{{"Probability of being sick after having tested\npositive for a disease: P(sick|positive test)",SpanFromLeft},{"","positive test","negative test","Total"},{"sick",Item[a11,Background->LightRed],a12,a13},{"healthy",Item[a21,Background->Yellow],a22,a23},{"Total",a11+a21,a12+a22,a11+a12+a21+a22},{" ",SpanFromLeft},{Item[outstring],SpanFromLeft}},Frame->All,ItemSize->7];
Return[Text @ res];];

Manipulate[ProbabilityOfBeingSick[prior,recogdisease,recoghealth,ntot],
"How good is the test in detecting the disease?",{{recogdisease,0.95,"P(positive test| sick):"},0,1.,Appearance->"Labeled"},
"How good is the test in detecting the absence of the disease?",{{recoghealth,0.9,"P(negative test| healthy):"},0,1.,Appearance->"Labeled"},
Delimiter,
"How likely is the disease?",{{prior,0.01,"P(disease)"},0.,1.,Appearance->"Labeled"},
Delimiter,
"Number of people for reference",{{ntot,3000,"N"},1,10000,1,Appearance->"Labeled"},
ControlPlacement->Top,SaveDefinitions->True
]

http://demonstrations.wolfram.com/ProbabilityOfBeingSickAfterHavingTestedPositiveForADiseaseBa/
Contributed by: Frank Scherbaum
After work by: Gerd Gigerenzer

Even very good medical tests for diseases have—usually very small—probabilities for false positive and false negative results. For rare diseases, the number of healthy people getting a false positive test (yellow background) can become much larger than the number of sick people getting a positive test (pink background).

D.I.Y

2.2.3  Bayesian Spam Filters

For events E and F from a sample space S with p(E)≠0 and p(F)≠0, one has

p(F|E)=((p(E|F)p(F))/(p(E|F)p(F)+p(E|(OverBar[F]))p((OverBar[F]))))

The following text describes how to use this theorem to create a Bayesian spam filter. We will use Mathematica to implement such a filter.

A message is received containing the word w. The event S will be the event that the message is spam and the event E is the event that the message contains the word w. If we assume that a message is as likely to be spam as not, so that p(S)=p((OverBar[S]))=(1/2), then Bayes' Theorem tells us that the probability that the incoming message is spam given that it contains the word w is:

p(S|E)=((p(E|S))/(p(E|S)+p(E|(OverBar[S]))))

By estimating the conditional probabilities with empirical data, we can compute an estimate that the given message is spam.

Example (From K. H. Rosen: Discrete Mathematics and Its Applications 7th Edition: pg 477)
Suppose that a Bayesian spam filter is trained on a set of 500 spam messages and 200 messages that are not spam. The word “exciting” appears in 40 spam messages and in 25 messages that are not spam. Would an incoming message be rejected as spam if it contains the word “exciting” and the threshold for rejecting spam is 0.9?

p(S|E)=((40/500)/(40/500+25/200))
               = 0.3902439024390244` ≤0.9
Thus,will not be rejected

D.I.Y

2.3 Expected Value and Variance

2.3.1  Introduction

Summary

In this section we will explore Mathematica's abilities more closely and use probability distributions to explore the concepts of expected value and variance.
Mathematica provides the distribution GeometricDistribution, which takes one parameter, the probability of a “success”.

GeometricDistribution

Probability

We can use the Probability function to compute the probabilities of events. 
For example, the probability p(X=5) is computed by the following.

Probability[X==5,X\[Distributed]GeometricDistribution[1/4]]

Note that the Mathematica's definition of a geometric random variable differs slightly from the textbook's. 
The textbook defines the value of the geometric random variable, in terms of coin flips, to be the number of flips it takes to get a tails, where the parameter is the probability of tails. Mathematica's definition is that the value of the random variable is the number of heads that appear before tails comes up. So the probability p(X=k) is:

Probability[X==k,X\[Distributed]GeometricDistribution[1/4]]

Also, Mathematica have built-in function for Poisson Distribution.

PoissonDistribution

D.I.Y

Exercise 1

Probability[X==3,X\[Distributed]GeometricDistribution[1/5]]

Mean, Variance and Standard Deviation

Mathematica can also calculate the expected value, variance, and standard deviation of distributions. 
The Mean, Variance, and StandardDeviation functions all accept a distribution as the sole argument and return the statistic.

Mean

Variance

StandardDeviation

Example:

Mean[GeometricDistribution[1/4]]

Variance[GeometricDistribution[1/4]]

These functions can also compute symbolically:

StandardDeviation[GeometricDistribution[k]]

All three functions can also accept a list holding data. 
For example, to compute the mean of {1,3,4,7,11}:

Mean[{1,3,4,7,11}]

D.I.Y

Exercise 1

Mean[{1,2,3,7,11,13,15,17,19}]

Exercise 2

Variance[{1,2,3,7,11,13,15,17,19}]

Exercise 3

StandardDeviation[{1,2,3,7,11,13,15,17,19}]

Expectation

Alternately, you can use the Expectation function. Expectation is more flexible in that it allows you to form an expression using random variables as a first argument and specify the distributions of the random variables as the second argument. 
The syntax is similar to that of Probability, except the first argument is an expression, not an event.

Expectation

For example, the following will find the expected value of the geometric distribution with parameter 1/4.

Expectation[X,X\[Distributed]GeometricDistribution[1/4]]

Expectation allows us to find the expected value of more complicated expressions, such as E(X^2).

Expectation[X^2,X\[Distributed]GeometricDistribution[1/4]]

D.I.Y

Exercise 1

Expectation[k+1,k\[Distributed]GeometricDistribution[p]]

Recall

Recalling that the variance of a random variable is E((((X-E(X)))^2)), we can use Expectation to compute variance as well.

mean=Expectation[X,X\[Distributed]GeometricDistribution[1/4]]

Expectation[(X-mean)^2,X\[Distributed]GeometricDistribution[1/4]]

You can also mix distributions, as shown below.

distribution1=GeometricDistribution[1/4];

distribution2=BinomialDistribution[20,.3];

BinomialDistribution

mean2=Expectation[X+2*Y,{X\[Distributed]distribution1,Y\[Distributed]distribution2}]

Expectation[((X+2*Y)-mean2)^2,{X\[Distributed]distribution1,Y\[Distributed]distribution2}]

D.I.Y

Exercise 1

dist1=GeometricDistribution[1/5];
dist2=BinomialDistribution[30,0.2]

mean3=Expectation[Y-5 X,{X\[Distributed]dist1,Y\[Distributed]dist2}]

Expectation[((Y-5X)-mean3)^2,{X\[Distributed]dist1,Y\[Distributed]dist2}]