# Stat01-Geometric.wxmx 

## Preface

In Stat01-Geometric.wxmx we discuss the discrete Geometric (p) probability distribution and its 
application, using Maxima tools and methods. 

Edwin L. (Ted) Woollett
https://home.csulb.edu/~woollett/
April 1, 2024

## References

In our series Statistics with Maxima we have used some examples and explanations (with much 
editing and additions) from:

Ch. 3 Reagle & Salvatore [RS], Statistics and Econometrics, 2nd ed, Schaum's Outlines,
      2011,  McGraw Hill,

Ch. 8 Fred Senese [FS], Symbolic Mathematics for Chemists: A Guide for Chemists, 2019,
      Wiley,

Louis Lyons, Statistics for Nuclear and Particle Physics, 1986, Cambridge Univ. Press,

Luca Lista,  'Statistical Methods for Data Analysis in Particle Physics', 
                          Lecture Notes in Physics 909, 2016,  Springer-Verlag,

Frederick James, 'Statistical Methods in Experimental Physics', 2nd ed., 
                            2006, World Scientific.      

https://www.statology.org/geometric-distribution-real-life-examples/

https://online.stat.psu.edu/stat414/lesson/11/11.1

https://calcworkshop.com/discrete-probability-distribution/geometric-distribution/

https://www.statisticshowto.com/geometric-distribution/

https://minisham.redbrick.dcu.ie/CA/Notes/CA266/10_Geometric_Distribution.pdf

https://www.studysmarter.co.uk/explanations/math/statistics/geometric-distribution/

Homemade functions fll, head, tail, Lsum are useful for looking at long lists.

load (descriptive);
load (distrib); 
fpprintprec : 6$
ratprint : false$

fll ( aL) := [ first (aL), last (aL), length (aL) ]$
declare (fll, evfun)$
head(L) := if listp (L) then rest (L, - (length (L) - 3) ) else 
      error("Input to 'head' must be a list of expressions ")$
declare(head,evfun)$
tail (L) := if listp (L) then rest (L, length (L) - 3 )  else 
       error("Input to 'tail' must be a list of expressions ")$
declare(tail,evfun)$
Lsum (aList) := apply ("+", aList)$
declare (Lsum, evfun)$

## Discrete Distributions Defined

From www.investopdedia.com:

"A discrete probability distribution counts occurrences that have countable or finite outcomes.

Discrete distributions contrast with continuous distributions, where outcomes can fall anywhere
on a continuum.

Common examples of discrete distribution include the binomial, Poisson, and Bernoulli 
distributions.

These distributions often involve statistical analyses of "counts" or "how many times" an event 
occurs.

In finance, discrete distributions are used in options pricing and forecasting market shocks or 
recessions."
---------------------------------------------------------------------------------------------
From http://www.stat.yale.edu/Courses/1997-98/101/ranvar.htm:

"If a random variable can take only a finite number of distinct values, then it must be discrete.
Examples of discrete random variables include the number of children in a family, the Friday night 
attendance at a cinema, the number of patients in a doctor's surgery, the number of defective 
light bulbs in a box of ten."

## What Is a Discrete Random Variable?

Quoting Luca Lista (Sec. 1.1):
"Many processes in nature have uncertain outcomes. This means that their result cannot be 
predicted before the process occurs. A random process is a process that can be reproduced, 
to some extent, within some given boundary and initial conditions, but whose outcome is 
uncertain. This situation may be due to insufficient information about the process intrinsic 
dynamics which prevents to predict its outcome, or lack of sufficient accuracy in reproducing 
the initial conditions in order to ensure its exact reproducibility. Some processes like 
quantum mechanics phenomena have intrinsic randomness. This will lead to possibly different 
outcomes if the experiment is repeated several times, even if each time the initial conditions 
are exactly reproduced, within the possibility of control of the experimenter. Probability is a 
measurement of how favored one of the possible outcomes of such a random process is 
compared with any of the other possible outcomes."

A coin toss is "random" because we are ignorant of the 'initial conditions'. Repeated trials tell 
us something about how those initial conditions vary between trials

For the purposes of calculating things for experimental physics, we need physical probability. 
In particular we need 'frequentist probability':  
   "Probability is the frequency with which a particular outcome occurs in repeated trials."
    P = (number of occasions on which that outcome occurs)/(total number of measurements).

Quoting L. Lyons, Sec. 2.1:
"In many situations we deal with experiments in which the essential circumstances are kept
constant, and yet repititions of the experiment produce different results. Thus the result of an
individual measurement or trial may be unpredictable, and yet the possible results of a series 
of such measurements have a well defined distribution."

What about events that can't be repeated? They don't have probabilities.

Quoting [RS] Sec. 3.3:

"A random variable is a variable whose values are associated with some probability of being 
observed. A discrete (as opposed to continuous)  random variable is one that can assume only 
finite and distinct values.  The set of all possible values of a random variable and its associated 
probabilities is called a probability distribution. The sum of all probabilites equals 1."

Quoting 
https://saylordotorg.github.io/text_introductory-statistics/s08-discrete-random-variables.html,

"The probability distribution of a discrete random variable X is a listing of each possible value x 
taken by X along with the probability P(x) that X takes that value in one trial of the experiment.

The mean μ of a discrete random variable X is a number that indicates the average value of X
over numerous trials of the experiment. It is computed using the formula μ=Σx P(x).

The variance σ^2 and standard deviation σ of a discrete random variable X are numbers that 
indicate the variability of X over numerous trials of the experiment. They may be computed using 
the formula σ^2 = (Σx^2 P(x) ) − μ^2, taking the square root to obtain σ."

## Mean and Variance of Discrete Data Set

Consider a data set that contains M unique discrete values x_k, and assume the value x_k 
occurs with frequency f_k. Let N equal the sum of the frequencies.
    N = sum (f_k, k, 1, M).
The mean 
    <x> = sum (f_k*x_k, k, 1, M) / N.
The variance
       Var(x) = sum ( f_k* (x_k - <x>)^2, k, 1, M )/ N.
The standard deviation is the square root of the variance.

## The Discrete Geometric (p) distribution

Quoting 
https://www.statology.org/geometric-distribution-real-life-examples/

"The Geometric distribution is a probability distribution that is used to model the probability of 
experiencing a certain amount of failures before experiencing the first success in a series of 
Bernoulli trials."

"A Bernoulli trial is an experiment with only two possible outcomes – “success” or “failure” – 
and the probability of success is the same each time the experiment is conducted. An example 
of a Bernoulli trial is a coin flip. The coin can only land on two sides (we could call heads a
“success” and tails a “failure”) and the probability of success on each flip is 0.5, assuming the 
coin is fair."

"If a random variable X follows a geometric distribution, then the probability of experiencing k 
failures before experiencing the first success can be found by the following formula:

P(k) = p* (1-p)^k

where:

k: number of failures before first success
p: probability of success on each trial"
Note that a probability of an event is a number p such that 0 < p < 1.
The percent probability is 100*p.

The theoretical mean μ of a Geometric(p) random variable is (1-p)/p.
The theoretical variance of a Geometric(p) random variable is (1-p)/p^2.
One standard deviation σ of a Geometric(p) random variable is the square root of the variance.

Quoting 
https://stats.libretexts.org/Bookshelves/Probability_Theory/
              Probability_Mathematical_Statistics_and_Stochastic_
              Processes_(Siegrist)/11%3A_Bernoulli_Trials/11.03%3A_The_Geometric_Distribution

"..... if the first success has not occurred by trial number m, then the remaining number of 
trials needed to achieve the first success has the same distribution as the trial number of the 
first success in a fresh sequence of Bernoulli trials. In short, Bernoulli trials have no memory. 
This fact has implications for a gambler betting on Bernoulli trials (such as in the casino games 
roulette or craps). No betting strategy based on observations of past outcomes of the trials can 
possibly help the gambler."

### pdf_geometric (k, p)

The Maxima function pdf_geometric (k, p) returns the value at k of the probability function of a 
Geometric(p) random variable, with 0 < p <= 1.

The probability function is defined as p*(1 - p)^k. This is interpreted as the probability of k failures 
before the first success.

### Statology 1: Coin Tosses

This example is from:
       https://www.statology.org/geometric-distribution-real-life-examples/

"Suppose we want to know how many times we’ll have to flip a fair coin until it lands on heads. 
We can use the following formulas to determine the probability of experiencing 0, 1, 2, 3 failures, 
etc. before the coin lands on heads; the coin can experience 0 “failures” if it lands on heads on 
the first flip." We are defining "success" to be: the coin lands on heads. Which means "failure" 
is defined as tails.

For each coin flip, p = P(success), (1- p) = P(failure), the joint probability of k failures followed 
by one success factors into the product of the probability of each event (whether failure or 
success), since the flipping of the coin event is completely independent of the outcome of any 
of the other flipping events. 

Hence P(k failures and 1 success)  = P(F)*P(F)* ....*P(F)*P(S) = (1-p)^k * p.

Let's define a formal Maxima function of this form.

P(k, p) := p*(1 - p)^k$

For a balanced unbiased coin, p = 0.5.
We will compare the results produced by pdf_geometric (k,p) with P(k,p).

p:0.5$
mypoints :makelist ([k,  pdf_geometric (k, p)], k, 0, 4);
makelist ( [k, P(k, p)], k, 0, 4);

Note that p*(1- p)^0 = p = 0.5 for any value of p.

Since (1-p) = 0.5 = p for this example, starting with P(0) = 0.5, we get 
P(1) = 0.5*p = 0.5^2 = 0.25, and so on.

Let k = number of failures before the first success, x = number of trials needed to finally find success.
x = 1,2,3,....,  k = x - 1 = 0, 1, 2, 3...

The probability of experiencing 0 failures before coin lands on heads is 50%. (H) (x = 1, k = 0)
The probability of experiencing 1 failure before coin lands on heads is 25%. (T, H)  (x = 2, k = 1)
The probability of experiencing 2 failures before coin lands on heads is 12.5%. (T,T,H)  (x = 3, k = 2)
The probability of experiencing 3 failures before coin lands on heads is 6.25%.  (T,T,T,H)

wxdraw2d ( xrange = [-1, 5], yrange = [0, 0.6], points_joined = impulses,
    xlabel = "Number of Failures", ylabel = "Probability", 
    title = "Distribution of Failures Before Coin Lands on Heads",   
    background_color = light_gray, grid = true,key = "p = 0.5",
    line_width = 4, color = red, points (mypoints) )$

Let x be the number of independent trials required to finally get a success. 
Let k be the number of failures required to finally get success on the next trial.
Then x = 1, 2, 3, ... and k = 0, 1, 2, 3, ... and k = x - 1.

What is the probability P( k >=  0)? Answer: 1 (100%). P (k >= 0) should be interpreted as
the probability we will need a value of k greater than or equal to 0 before getting a success,
which is the "area under the curve" from k = 0 to k = inf. For a discrete distribution, the
continuous distribution concept "area under the curve" becomes the sum of the discrete 
probability values given by pdf_geometric (k, p) in some range [k1, k2].

The probability x >= 1 is unity: P(x >= 1) = 1. Note P(x >= 1) should be interpreted as the probability
we will need a value of trials (x) greater than or equal to 1, but x starts at x = 1 and continues to 
arbitrarily large integral values. 

Now x >= 1 ==>  x - 1 >= 0 ==> k >=  0. So P(k >=0) = P(x >= 1)

What is the probability P( k > 0)? What is the probability we will need one or more failures before the
next trial will result in success? 
Answer:  P(k>0) =  P(k >= 0) - P(k = 0)  = 1 - 0.5 = 0.5 (50%).

What is the probability P (k < 2)?  What is the probability we will need less than two failures before the
next trial will result in success? This is "the area under the curve" from k = 0 to k = 1, ie., the sum
of the discrete probabilities P(0) and P(1).

Answer:  P(k < 2) = P(0) + P(1) = 0.5 + 0.25 = 0.75 (75%).

### cdf_geometric (k, p)

The 'cumulative distribution function' (cdf): cdf_geometric (k, p) returns the sum of the values P(k)
for a Geometric (p) random variable, with 0 < p < 1,  beginning at k=0 and ending at k.
P (k <= n) = cdf_geometric (n, p)

Continuing with the fair coin toss example, we then get P(k  < 2) = P(0) + P(1)  = 0.5 + 0.25 = 0.75
from cdf_geometric (1, 0.5).

cdf_geometric (1, 0.5);

What is the probability P(k <=2)?  Answer: P(0) + P(1) + P(2) = 0.75 + 0.125 = 0.875 (87.5%).

cdf_geometric (2, 0.5);

What is the probability P(k > 2)?  Answer: 1 - P(k <= 2) = 1 - 0.875 = 0.125 (12.5%).

1 - cdf_geometric (2, 0.5);

What is the probability P(k>=2)?  Answer:  1 - P(k < 2) = 1 - 0.75 = 0.25 (25%).

1 - cdf_geometric (1, 0.5);

In summary, we have:

P(k = 2): 0.12500   pdf_geometric (2, 0.5);

P(k < 2): 0.75000   cdf_geometric (1, 0.5);

P(k ≤ 2): 0.87500    cdf_geometric (2, 0.5);

P(k > 2): 0.12500     1 - cdf_geometric (2, 0.5);

P(k ≥ 2): 0.25000     1 - cdf_geometric (1, 0.5);

### Maxima Function    geomCalc (k, p)

We can incorporate the methods in the above calculation to create a Maxima function
geomCalc (k, p). We could, of course, have multiplied by 100, to list percent probabilities instead of
probabilities.

geomCalc (kk, pp) := 
block ( 
    if not numberp (pp) then (
           print ("pp must be a number such that 0 < pp < 1"),
           return (false) ),
    if pp <= 0 or pp >= 1 then (
          print ("Need 0 < pp < 1"),
          return (false) ),
    pp : float (pp),
    print (" "),
    print (sconcat (" For p = ", pp, " the probabilities are:")),
    print ("  "),
    print (sconcat (" P ( k = ", kk, " )   : ", pdf_geometric (kk, pp) )), 
    print (sconcat (" P ( k < ", kk, " )   : ", cdf_geometric (kk - 1 , pp) )), 
    print (sconcat (" P ( k <= ", kk, " ) : ", cdf_geometric (kk, pp) )), 
    print (sconcat (" P ( k > ", kk, " )   : ", 1 -   cdf_geometric (kk, pp) )), 
    print (sconcat (" P ( k >= ", kk, " ) : ", 1 -  cdf_geometric (kk - 1, pp) )), 
    done )$

geomCalc (2, 0.5)$

geomCalc (1, 0.5)$

geomCalc (0, 0.5)$

geomCalc (2, x)$

geomCalc (2, 3)$

These results agree with the results returned by
      https://www.statology.org/geometric-distribution-calculator/

Instead of using cdf_geometric is various ways to construct calcGeom, we could have used formulas
directly related to p, such as P (k <= m) = 1 - (1 - p)^(m+1), P (k < m) = 1 - (1 - p)^m,
P (k >= m) = (1 - p)^m, P (k > m) = (1 - p)^(m+1). Note that cdf_geometric (m, p)  = P (k <= m).

### Statology 2: Number of Bankruptcies

This example is from:
       https://www.statology.org/geometric-distribution-real-life-examples/

"Suppose it’s known that 4% of individuals who visit a certain bank are visiting to file bankruptcy. 

Suppose a  banker wants to know the probability that he will meet with 10 people 
before encountering someone who is filing for bankruptcy.  We are defining "success" to be:
 the banker meets someone wishing to file bankruptcy. We are defining "failure" to be: the 
banker meets someone not wishing to file for bankruptcy. So in this possible scenario, the 
10th person the banker meets wants to file for bankrupty, meaning the preceding 9 persons are
"failures". What is the probability of that happening?

cdf_geometric (9, 0.04), numer;

Thus there is a 33.5% probability the banker will meet with 10 people before encountering someone
wishing to file for bankruptcy.

geomCalc (9, 0.04)$

### PSU Example 11.1

We adapt a problem from:
https://online.stat.psu.edu/stat414/lesson/11/11.1

A reporter is looking on the streets of Kansas City, selecting people until he finds the first person who
attended the most recent Kansas City Chiefs home game. Suppose the probability of selecting such
a person is p = 0.20 = P(success)  for each person met. Then the probability of failure for each person
met is (1 - p)  = P(failure) = 0.80. 

Let x = number selected (with the last one success), and k = number of failures before success.
x = 1, 2, 3, ..., and k = 0, 1, 2, .....
So k = x - 1. The condition x > 6 ==>  x - 1 > 6 - 1 ==>  k > 5, etc.

A.) What is the probability he must select 4 people before finding such a person. Selecting 4 
people one after the other before finding success means the first three people selected 
represent failures. x = 4 ==> k = 3.  (F,F,F,S)

pdf_geometric (3, 0.20), numer;

Hence about a 10% chance the reporter must select 4 people before finding one who attended the
last home game.

B.) What is the probability that the reporter must select more than 6 people before he finds one who 
attended the last home football game?

Ans:  P (x > 6) = 1 - P (x <= 6) = 1 - P (k <= 5) = 1 - cdf_geometric (5, 0.2)

1 - cdf_geometric (5, 0.2), numer;

There is about a 26% chance that the reporter would have to select more than 6 people before finding
one who attended the last home game.

geomCalc (5, 0.2)$

### mean_geometric (p)

The Maxima function mean_geometric (p) returns the mean of a Geometric(p) random variable, 
with 0 < p <= 1. The theoretical mean of Geometric(p) [with k = number of failures before 
success,   k = 0, 1, 2,...] is given by E(k) =  (1 - p)/p, which is returned by Maxima's mean_geometric
function.
The theoretical mean number of **trials** x requred (the last being success) is E(x) =  1/p.

Calculating E(k) by hand as (1 - p)/p, with p = 0.2, E(k) = 4, E(x) = 1/p = 5.

0.8/0.2;

Using Maxima's mean_geometric function:

mean_geometric (0.2);

The mean value of a geometric distribution with p = 0.2 is k  = 4 failures, or x = 5 independent trials.
We should expect the reporter to have to select 5 people before he finds one who attended the last
home game. Of course, on any given try, it may take 1 person or it may take 10, but 5 is the 
average number.

### var_geometric (p)

The Maxima function var_geometric (p) returns the variance of a Geometric(p) random variable, 
with 0 < p <= 1. The theoretical variance of a Geometric (p) random variable [with k = number 
of failures before success,  k = 0,1,2,...]  is Var (k) = (1 - p)/p^2. This is also Var (x).

Calculating the variance in the number of failures for p = 0.2 by hand:

0.8/ (0.2)^2;

var_geometric (0.2);

The square root of the variance V(k) is one standard deviation about the mean.

sqrt (%);

### std_geometric (p)

 The Maxima function std_geometric (p) returns the standard deviation of a Geometric(p) 
random variable, with 0 < p <= 1. This is the square root of the variance, or sqrt(1-p)/p. 

sqrt(0.8)/0.2;

std_geometric (0.2);

### Statology 3:  Supporters of a Law

This example is from:
       https://www.statology.org/geometric-distribution-real-life-examples/

"Suppose a researcher is waiting outside of a library to ask people if they support a certain law. 
The probability that a given person supports the law is p = 0.2."

"What is "the probability of interviewing 0, 1, 2 people, etc. before the researcher speaks with 
someone who supports the law?"  We are defining success as finding a person who supports 
the particular law.

x = number of people interviewed. P(x = 1) = P(k = 0)
Let k be the number of failures before finding a person who supports the law when the chance
that any given person interviewed supports the law is 20%.

The probability of 0 failures before finding a person who supports the law is then

pdf_geometric (0,0.2);

P(k = 0) = P(x = 1) = 0.2 (20% chance).

P(x = 2) = P(k = 1).

pdf_geometric (1, 0.2);

and so on.

mypoints : makelist ([k, pdf_geometric (k, 0.2)], k, 0, 4);

The probability that the first person (0 failures) the researcher speaks to supports the law is 
20% (1 in 5).
The probability that the researcher must interview 2 people (one failure) before finding 
someone who  supports the law is 16%.
The probability that the researcher must interview 3 people (two failures) before finding 
someone who  supports the law is 12.8%.

wxdraw2d ( xrange = [-1, 5], yrange = [0, 0.25], points_joined = impulses,
    xlabel = "Number of Failures", ylabel = "Probability", 
    title = "Distribution of Failures Before Finding Supporter of Law",   
    key  = "p = 0.2",
    background_color = light_gray, grid = true,
    line_width = 4, color = red, points (mypoints) )$

We can use our Maxima function geomCalc (k, p) defined in Example 1:

### Statology 4: Chances of Network Failures

This example is from:
       https://www.statology.org/geometric-distribution-real-life-examples/

"Suppose it’s known that the chance that a certain company experiences a network failure 
in a given week is 10%. Suppose the CEO of the company would like to know the chance 
that the company can go 5 weeks or longer without experiencing a network failure."

p = P(success) = 0.1 = probability that a week will go by with a network failure.
(1-p) = P(failure) = 0.9 = probability that a week will go by without a network failure.
We need to calculate P (k >= 5).

geomCalc (5, 0.1)$

Thus the chance the company can go five weeks or longer WITHOUT experiencing a 
network failure is 59%.

### Calcwork Ex. 1

We summarize a problem from:

https://calcworkshop.com/discrete-probability-distribution/geometric-distribution/

"Suppose Max owns a lightbulb manufacturing company and determines that 3 out of every 75 bulbs 
are defective. What is the probability that Max will find the first faulty lightbulb on the 6th one that 
he tested?"

If P(success) = P(finding defective lightbulb) = p = 3/75 = 0.04, the P(failure) = (1 - p) = 0.96.
P(x = 6) = P (k = 5) = pdf_geometric (5, 0.04).

pdf_geometric (5, 0.04), numer;

About a 3% chance of finding the first faulty lightbulb on the 6th one tested.

"Now, what if Max wants to know the likelihood that it takes at least six trials until he finds the 
first defective lightbulb?"

P (x >= 6) = P (k >= 5) = 1 - cdf_geometric (4, 0.04).

1 - cdf_geometric (4, 0.04);

geomCalc (5, 0.04)$

"... there is an 0.815 chance of Max needing at least six trials until he finds the first defective lightbulb."

"And using this same example, let’s determine the number lightbulbs we would expect Max to inspect 
until he finds his first defective, as well as the standard deviation."

Expectation value for the number of lightbulbs one must test before finding a faulty one is E(x) = 1/p.
Expectation value for the number of failures before the next lightbulb tested is faulty is
E(k) = (1 - p)/p

1/0.04;
0.96/0.04;

mean_geometric (0.04);

std_geometric (0.04);

"This shows us that we would expect Max to inspect 25 lightbulbs before finding his first defective, 
with a standard error of 24.49."

### Stat-howto Ex. 1

A reporter is looking for people who voted in the last election and voted as an independent.
The probability of finding such a voter is P(success) = p  = 0.2. What is the probability you would
meet an independent voter on your third try (x = 3, k = 2) 

pdf_geometric (2, 0.2), numer;

### Inverse type of problem.

We summarize Ex. 3 from
https://minisham.redbrick.dcu.ie/CA/Notes/CA266/10_Geometric_Distribution.pdf

"It is known that 20% of products on a production line are defective. Products are inspected until first
defective is encountered. Let x = number of inspections to obtain first defective.
What is the minimum number of inspections that would be necessary so that the probability of 
observing a defective is more that 75%?

Choose n so that P(x ≤ n) ≥ .75

### quantile_geometric (q, p)

quantile_geometric (q, p) returns the q-quantile (an integer) of a Geometric (p) random variable, 
with 0 < q < 1, 0 < p < 1.

This is the inverse of cdf_geometric (k, p).

The probability from which the quantile is derived is defined as p*(1 - p)^k, interpreted as
the probability of k failures before the first success.

quantile_geometric (q, p) is not strictly the inverse of cdf_geometric (k, p), since
quantile_geometric (q, p) is only defined for values of q which belong to the set defined by
cdf_geometric (k, p) with k = 0, 1, ... (non-negative integers). In the following we construct that
set and display floating point values of qq as the second element of
  [k, cdf_geometric (k, 0.2), quantile_geometric (q, 0.2)].

Returning to the problem framed above: Choose n so that P(x ≤ n) ≥ .75

quantile_geometric (0.75, 0.2);

for k:4 thru 7 do (
    qq : float (cdf_geometric (k, 0.2)),
     print ([k, qq, quantile_geometric (qq, 0.2)] ) )$

We see that the minimum number of inspections such that the probability of finding a defective
product is equal to or greater than 0.75 is 6 (which value has 0.79 rather that 0.75).

Continuing with the example, What is the average number of inspections to obtain the first defective?

The average number of trials (inspections) E(x) = 1/p and the average number of failures before
finding a defective is E(k) = (1-p)/p

1/0.2;
0.8/0.2;

mean_geometric (0.2);

The average number of inspections before finding a defective is 5. (x = 5, k = 4 failures).

### Kidney Donor Scenario

A problem from
https://www.studysmarter.co.uk/explanations/math/statistics/geometric-distribution/

"A patient suffers kidney failure and requires a transplant from a suitable donor. The probability 
that a random donor will match this patient’s requirements is 0.2.

1.) Suppose that no donor matches the patient's requirements until a fifth donor comes in. 
            What is the probability of this scenario?

P(success) = p = 0.2 = probability a random donor will match patient's need. P(x = 5) = P(k = 4).

pdf_geometric (4, 0.2), numer;

The probability the fifth donor will provide the needed kidney is about 8.2%

2.) Find the probability of the patient requiring 10 or fewer donors until a match is found.

P(x <= 10) = P(k <= 9)

geomCalc (9, 0.2)$

About a 90% chance a match will be found after 10 donors.

3.)What is the expected number of donors required to get a match?

The expected number of donors E(x) = 1/p = 1/0.2 = 5

4.) Find the standard deviation of this scenario.

std_geometric (0.2);

### Rolling a Fair Dice

A problem from
https://www.studysmarter.co.uk/explanations/math/statistics/geometric-distribution/

"Suppose you roll a fair dice until you get a three as a result."

The probability of rolling a 3 is P(success) = p = 1/6 = 0.1667, P(failure) = 1-p = 5/6 = 0.8333.

1.) "What is the probability that you don't roll a three until your fourth roll?"

P(x=4) = P(k = 3)

pdf_geometric (3, 1/6), numer;

9.6% ~ 10% chance won't get a three until the fourth try.

2.) "Find the probability of getting the three you need in less than 10 rolls."

P (x < 10) = P (k < 9).

geomCalc (9, 1/6);

About 81% chance you will get the three you need in less than 10 rolls.

3.) "What is the expected number of rolls required to get your desired outcome?"

E(x) = 1/p = 1/(1/6) = 6 rolls.

4.) "Find the variance of this experiment."

var_geometric (1/6), numer;

sqrt(%);

### Winning a Prize From a Claw Machine

A problem from
https://www.studysmarter.co.uk/explanations/math/statistics/geometric-distribution/

"Suppose that the probability of winning an item from a claw machine is 0.05."

p = 0.05 = P(success) on each attempt to extract a prize using a claw machine.

1.) "What is the probability of winning an item on your first try?"

P(x = 1) = P(k = 0).

pdf_geometric (0, 0.05), numer;

2.) What is the probability of winning an item in less than 20 tries?"

P(x < 20) = P(k < 19).

geomCalc (19, 0.05);

About a 62% chance you will win a prize in less than 20 tries.

3.) Suppose you need to use a quarter for each try. What is the expected amount 
       of money spent for getting a prize?"

The average number of tries before winning a prize is E(x) = 1/p = 1/0.05 = 20.
(20 tries)*(0.25$/try) = 20 * (1/4) $ = 5$. Expect the pay about $5 per prize won.

### random_geometric (p), random_geometric (p, n)

The Maxima function random_geometric (p) returns the number of failures (k) in one 
Bernoulli trial in which the probability of "success" is p.

The Maxima function random_geometric (p,n) returns a list of n calls to random_geometric (p).

The probability from which the random sample is derived is defined as p (1 - p)^k. This is 
interpreted  as the probability of k failures before the first success.

Suppose p = 0.2; the probability of success in each trial is 20%. 

random_geometric (0.2);

for j thru 5 do print (random_geometric (0.2))$

To obtain the number of failures before success in each of 10 Bernoulli trials, given that the
chance of success in each such trial is 20%, we would use random_geometric (0.2, 10).

for j thru 5 do print (random_geometric (0.2, 10))$

### 100 Simulated Bernoulli Trials, p = 0.2

Let data be a list of the number of failures before success in 100 simulated Bernoulli trials,
assuming p = 0.2.

data : random_geometric (0.2, 100)$
fll (data);
head (data);
tail (data);

### discrete_freq (Data)

The Maxima function discrete_freq( aList) counts the number of unique discrete "readings" 
of some instrument recorded in the list aList and returns a new list [ list-of-unique-readings, 
list-of-frequency-of-each-unique-reading]. The elements of the list 'frequencies' corresponds 
to the elements of the list 'uniqueData', element by element. 

[uniqueData, frequencies] : discrete_freq (data)$
uniqueData;
frequencies;

length (uniqueData);

Out of a list of 100 semi-random integers, there is a lot of overlap, resulting in only a few unique
integers.  

nfrequencies : frequencies/Lsum (frequencies), numer;

Lsum (frequencies);

Lsum (nfrequencies);

lmax (nfrequencies);

mypoints : makelist ([uniqueData[j], nfrequencies[j]], j, 1, length (uniqueData))$
fll (mypoints);

wxdraw2d ( xrange = [-1, 26], yrange = [0, 0.2], points_joined = impulses,
    xlabel = "Number of Failures before Success", ylabel = "Probability", 
    title = "100 simulated Bernoulli trials with p = 0.2",       
    background_color = light_gray, grid = true,
    line_width = 4, color = red, points (mypoints),
    color = black, line_width = 2, 
    key = "pdf geometric (k, 0.2)",
    explicit (pdf_geometric (kk, 0.2), kk, 0, 24)  )$

### 1000 Simulated Bernoulli Trials, p = 0.2

We should get results closer to the theoretical prediction if we simulate 1000 Bernoulli trials
instead of 100. Let data be a list of the number of failures before success in 1000 simulated 
Bernoulli trials, assuming p = 0.2.

data : random_geometric (0.2, 1000)$
fll (data);

[uniqueData, frequencies] : discrete_freq (data)$
uniqueData;
frequencies;

nfrequencies : frequencies/Lsum (frequencies), numer;

Lsum (nfrequencies);

lmax (nfrequencies);

mypoints : makelist ([uniqueData[j], nfrequencies[j]], j, 1, length (uniqueData))$
fll (mypoints);

wxdraw2d ( xrange = [-1, 17], yrange = [0, 0.2], points_joined = impulses,
    xlabel = "Number of Failures before Success", ylabel = "Probability", 
    title = "random geometric (0.2, 1000)",       
    background_color = light_gray,    grid = true,
    line_width = 4, color = red, points (mypoints),
    color = black, line_width = 2, 
    key = "pdf geometric (k, 0.2)",
    explicit (pdf_geometric (kk, 0.2), kk, 0, 17)  )$