Probability: The Logic of Science

Jaynes presents a few core ideas and requirements for his rational system. Probability emerges as the representation of circumstances in which any given realization of a process is either TRUE or FALSE but both are possible and can be expressed by probabilities.

• that sum to one for all events
• are greater than or equal to zero for any given event

Convex Combinations

The formal definition is a linear combination with non-negative coefficients that sum to one. Sound familiar? This reasoning applies pretty broadly.

Weighted averages….

Decision Trees

Characterizing expected value for risk neutral agents. Decision trees [radiant has one] formalize this idea. We can also just compute them.

Representing Probability Distributions

Probability Distributions of Two Forms

Our core concept is a probability distribution just as above. These come in two forms for two types [discrete (qualitative)] and continuous (quantitative)] and can be either:

• Assumed, or
• Derived

The Poster and Examples

• Distributions are nouns.

• Sentences are incomplete without verbs – parameters.

• We need both; it is for this reason that the former slide is true.

• We do not always have a grounding for either the name or the parameter.

• For now, we will work with univariate distributions though multivariate distributions do exist.

Continuous vs. Discrete Distributions

The differences are sums versus integrals. Why?

• Histograms or
  
• Density Plots

The probability of exactly any given value is zero on a true continuum.

Functions

Probability distributions are mathematical formulae expressing likelihood for some set of qualities or quantities.

• They have names: nouns.
  
• They also have verbs: parameters.

Like a proper English sentence, both are required.

Models

What is a model?

Our Applications

• The uniform is defined by a minimum and maximum.
  
• The normal with mean \mu and standard deviation \sigma or variance - \sigma^2.
• The Poisson will be defined an arrival rate \lambda – lambda.
  
• Bernoulli trials: Two outcomes occur with probability \pi and 1-\pi.
  • The binomial distribution will be defined by a number of trials n and a probability \pi.
    
  • The geometric distribution defines the first success of n trials with probability \pi.
    
  • The negative binomial distribution defines the probability of k successes in n trials with probability \pi. It is related to the Poisson.

The Uniform Distribution

• Is flat, each value is equally likely.
  
• Defined on 0 to 1 gives a random cumulative probability X \leq x.
  
• \uparrow It’s a random probability in 0 and 1.
  
• Is also known as the rectangular distribution.

Uniform(0,1)

library(patchwork)
Unif <- data.frame(x=seq(0, 1, by = 0.005)) %>% mutate(p.x = punif(x), d.x = dunif(x))
p1 <- ggplot(Unif) + aes(x=x, y=p.x) + geom_step() + labs(title="Distribution Function [cdf/cmf]") + theme_minimal()
p2 <- ggplot(Unif) + aes(x=x, y=d.x) + geom_step() + labs(title="Density Function [pdf/pmf]") + theme_minimal()
p2 + p1

The Normal [Gaussian]

f(x|\mu,\sigma^2 ) = \frac{1}{\sqrt{2\pi\sigma^{2}}} \exp \left[ -\frac{1}{2} \left(\frac{x - \mu}{\sigma}\right)^{2}\right]

Is the workhorse of statistics. Key features:

• Is self-replicating: sums of normals are normal.
  
• If X is normal, then Z = \frac{(X - \mu)}{\sigma} is normal.

The Normal [Plotted]

library(patchwork)
Unif <- data.frame(x=seq(0, 1, by = 0.005)) %>% mutate(p.x = punif(x), d.x = dunif(x))
p1 <- ggplot(Unif) + aes(x=x, y=p.x) + geom_step() + labs(title="Distribution Function [cdf/cmf]") + theme_minimal()
p2 <- ggplot(Unif) + aes(x=x, y=d.x) + geom_step() + labs(title="Density Function [pdf/pmf]") + theme_minimal()
p2 + p1

Sample z-scores

The scale command in R does this for a sample.

z = \frac{x - \overline{x}}{s_{x}} where \overline{x} is the sample mean of x and s_{x} is the sample standard deviation of x.

In samples, the 0 and 1 are exact; these are features of the mean and degrees of freedom. If I know the mean and any n-1 observations, the n^{th} observation is exactly the value such that the deviations add up to zero/cancel out.

Z and symmetry

z is an easy way to assess symmetry.

• The mean of z is always zero but the distribution of z to the left and right of zero is informative. If they are roughly even, then symmetry is likely.
  
• If the signs are uneven, then symmetry is unlikely.
  
• In R, z is automated with the scale() command. The last line uses a table and the sign command to show the positive and negative z.

# Generate random normal income
DataF <- data.frame(Hypo.Income = rnorm(1000, 55000, 10000)) %>%
# z-transform income [mean 55000ish, std. dev. 10000ish]
mutate(z.Income = scale(Hypo.Income))
# Show the data.frame
head(DataF)
table(sign(DataF$z.Income))

  Hypo.Income   z.Income
1    55994.88  0.1198684
2    46199.58 -0.8675543
3    44239.77 -1.0651147
4    49542.55 -0.5305633
5    62587.74  0.7844678
6    52400.56 -0.2424600

 -1   1 
495 505 

Why Normals?

• The Central Limit Theorem
• They Dominate Ops [6\sigma]
• Normal Approximations Abound

Normals

The Generic Bernoulli Trial

Suppose the variable of interest is discrete and takes only two values: yes and no. For example, is a customer satisfied with the outcomes of a given service visit?

For each individual, because the probability of yes (1) \pi and no (0) 1-\pi must sum to one, we can write:

f(x|\pi) = \pi^{x}(1-\pi)^{1-x}

The Binomial

BinomialR

A Nicer Plot

library(plotly)
p <- ggplotly(Plot1)
p

The Rule of Five

• This is called the Rule of Five by Hubbard in his How to Measure Anything.

Geometric Distributions

How many failures before the first success? Now defined exclusively by p. In each case, (1-p) happens k times. Then, on the k+1^{th} try, p. Note 0 failures can happen…

Pr(y=k) = (1-p)^{k}p

Example: Entrepreneurs

Suppose any startup has a p=0.1 chance of success. How many failures?

Negative Binomial Distributions

How many failures before the r^{th} success? In each case, (1-p) happens k times. Then, on the k+1^{th} try, we get our r^{th} p. Note 0 failures can happen…

Pr(y=k) = {k+r-1 \choose r-1}(1-p)^{k}p^{r}

Needed Sales

I need to make 5 sales to close for the day. How many potential customers will I have to have to get those five sales when each customer purchases with probability 0.2.

library(patchwork)
DF <-  data.frame(Customers = c(0:70)) %>% 
  mutate(m.Customers = dnbinom(Customers, size=5, prob=0.2), 
         p.Customers = pnbinom(Customers, size=5, prob=0.2)) 
pl1 <- DF %>% ggplot() + aes(x=Customers) + geom_line(aes(y=p.Customers)) 
pl2 <- DF %>% ggplot() + aes(x=Customers) + geom_point(aes(y=m.Customers))

Events: The Poisson

Poisson

Take a binomial with p very small and let n \rightarrow \infty. We get the Poisson distribution (y) given an arrival rate \lambda specified in events per period.

f(y|\lambda) = \frac{\lambda^{y}e^{-\lambda}}{y!}

Examples: The Poisson

• Walk in customers
• Emergency Room Arrivals
• Births, deaths, marriages
• Prussian Cavalry Deaths by Horse Kick
• Fish?

Air Traffic Controllers

FAA Decision: Expend or do not expend scarce resources investigating claimed staffing shortages at the Cleveland Air Route Traffic Control Center.

Essential facts: The Cleveland ARTCC is the US’s busiest in routing cross-country air traffic. In mid-August of 1998, it was reported that the first week of August experienced 3 errors in a one week period; an error occurs when flights come within five miles of one another by horizontal distance or 2000 feet by vertical distance. The Controller’s union claims a staffing shortage though other factors could be responsible. 21 errors per year (21/52 errors per week) has been the norm in Cleveland for over a decade.

Deaths by Horse Kick in the Prussian cavalry?

library(vcd)
data(VonBort)
head(VonBort)

  deaths year corps fisher
1      0 1875     G     no
2      0 1875     I     no
3      0 1875    II    yes
4      0 1875   III    yes
5      0 1875    IV    yes
6      0 1875     V    yes

mean(VonBort$deaths)

[1] 0.7

Simulation: A Powerful Tool

In the geometric example, I was concerned with sales. I might also want to generate revenues because I know the rough mean and standard deviation of sales. Combining such things together forms the basis of a Monte Carlo simulation. More on this in a bit, it is known as calibration.

An Example

Customers arrive at a rate of 7 per hour. You convert customers to buyers at a rate of 85%. Buyers spend, on average 600 dollars with a standard deviation of 150 dollars.

Sim <- 1:1000
Customers <- rpois(1000, 7)
Buyers <- rbinom(1000, size=Customers, prob = 0.85)
Data <- data.frame(Sim, Buyers, Customers)
Data <- Data %>% group_by(Sim) %>% mutate(Revenue = sum(rnorm(Buyers, 600, 150))) %>% ungroup()

A Summary

Distributions are how variables and probability relate. They are a graph that we can enter in two ways. From the probability side to solve for values or from values to solve for probability. It is always a function of the graph.

Distributions generally have to be sentences.

• The name is a noun but it also has
• parameters – verbs – that makes the noun tangible.

Monte Carlo Simulation Generally

• Election Forecasts and Calibration
• More generally, on calibration
• Sensitivity analysis
• One way to think about much of this is that we have many spreadsheet models, the twist is to think about the proper role of what we are uncertain about and how uncertain we are.