How to calculate probabilities with SciPy distributions

Many statistical questions reduce to an area under a model: the chance of observing a value below a limit, above a limit, or inside a range. A parameterized SciPy distribution calculates those areas without evaluating or integrating the density by hand.

A frozen scipy.stats.norm object stores a mean of 70 and a standard deviation of 8. Its cdf() method returns lower-tail probability, sf() returns upper-tail probability, and the difference between two cumulative values gives the probability inside an interval.

The ppf() method moves in the opposite direction by converting a cumulative probability into a cutoff. For a continuous distribution, pdf() is density rather than probability at one exact value, so threshold and interval questions belong to the cumulative methods.

1. Create probability_distribution.py with the imports and frozen normal distribution.

   probability_distribution.py

   import numpy as np
   from scipy.stats import norm
    
   scores = norm(loc=70, scale=8)

   loc sets the normal distribution's mean, and scale sets its standard deviation. Freezing them once keeps every later calculation on the same model.

2. Add lower-tail and upper-tail probability calculations below the distribution.

   below_75 = scores.cdf(75)
   above_85 = scores.sf(85)

   sf() evaluates the upper tail directly and can be more accurate than subtracting cdf() from 1 when the cumulative probability is close to 1.

3. Add the probability between scores 65 and 80 below the tail calculations.

   between_65_and_80 = scores.cdf(80) - scores.cdf(65)

   For a continuous distribution, subtracting the cumulative probability at the lower boundary from the cumulative probability at the upper boundary gives the area inside the interval.

4. Add the 90th-percentile cutoff below the interval calculation.

   percentile_90 = scores.ppf(0.90)

   ppf(0.90) returns the score with 90 percent of the distribution at or below it.

5. Add a cumulative round-trip assertion below the percentile calculation.

   probability_levels = np.array([0.1, 0.5, 0.9])
   round_trip = scores.cdf(scores.ppf(probability_levels))
   np.testing.assert_allclose(round_trip, probability_levels, atol=1e-12)

   The assertion raises an error if converting the probability levels to cutoffs and back does not recover the original values within floating-point tolerance.

6. Append the calculated probabilities and cutoff below the assertion.

   print(f"P(X <= 75): {below_75:.4f}")
   print(f"P(X > 85): {above_85:.4f}")
   print(f"P(65 < X <= 80): {between_65_and_80:.4f}")
   print(f"90th percentile: {percentile_90:.2f}")
   print("CDF round trip:", np.round(round_trip, 1))

7. Run the completed probability calculation to verify the cumulative round trip.

   $ python3 probability_distribution.py
   P(X <= 75): 0.7340
   P(X > 85): 0.0304
   P(65 < X <= 80): 0.6284
   90th percentile: 80.25
   CDF round trip: [0.1 0.5 0.9]

   The command exits before printing the results if the cumulative and inverse-cumulative round trip fails.