# --- Setup: load our helpers and the data --------------------------------
import sys, pathlib
for _p in [pathlib.Path.cwd(), pathlib.Path.cwd().parent]:
    if (_p / 'lib' / 'statslab.py').exists():
        sys.path.insert(0, str(_p / 'lib')); break
import statslab as sl

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
sl.use_course_style()

# The SAME two datasets as the core lesson:
faithful = sl.load_csv('faithful.csv',
                       url='https://vincentarelbundock.github.io/Rdatasets/csv/datasets/faithful.csv')
loans = sl.load_csv('loans.csv',
                    url='https://www.openintro.org/data/csv/loans_full_schema.csv')

print(f'Old Faithful: a {type(faithful).__name__} with shape {faithful.shape}  (rows, columns)')
print(f'Lending Club: shape {loans.shape}')

Old Faithful: a DataFrame with shape (272, 3)  (rows, columns)
Lending Club: shape (10000, 55)

print('shape  :', faithful.shape)
print('columns:', list(faithful.columns))
print('\ndtypes (the type stored in each column):')
print(faithful.dtypes)
faithful.head()

shape  : (272, 3)
columns: ['rownames', 'eruptions', 'waiting']

dtypes (the type stored in each column):
rownames       int64
eruptions    float64
waiting        int64
dtype: object

e = faithful['eruptions']                  # eruption length, minutes — a Series
print('e is a', type(e).__name__, 'of', len(e), 'values; first 5:')
print(e.head().to_list())

# A boolean mask in action: the "long" eruptions.
long_mask = e > 3
print(f'\neruptions longer than 3 min: {long_mask.sum()} of {len(e)}  ({long_mask.mean()*100:.0f}%)')
print('(.sum() counts the Trues; .mean() turns them into a proportion.)')

e is a Series of 272 values; first 5:
[3.6, 1.8, 3.333, 2.283, 4.533]

eruptions longer than 3 min: 175 of 272  (64%)
(.sum() counts the Trues; .mean() turns them into a proportion.)

edges = np.arange(1.5, 5.5 + 0.5, 0.5)          # bin edges every half-minute: 1.5, 2.0, ..., 5.5
print('bin edges:', edges)

counts, _ = np.histogram(e, bins=edges)         # how many eruptions land in each bin
print('counts (np.histogram):', counts, '  total =', counts.sum())

# The pandas way to the SAME tallies, as a labelled table:
table = pd.cut(e, bins=edges, right=False).value_counts().sort_index()
print('\ncounts (pd.cut):', table.to_list(), '  total =', int(table.sum()))
print('same counts?', np.array_equal(counts, table.to_numpy()))

freq = table.rename('count').to_frame()
freq['bar'] = freq['count'].apply(lambda c: '#' * c)   # a tiny text histogram
print('\n' + freq.to_string())

bin edges: [1.5 2.  2.5 3.  3.5 4.  4.5 5.  5.5]
counts (np.histogram): [51 41  5  7 30 73 61  4]   total = 272

counts (pd.cut): [51, 41, 5, 7, 30, 73, 61, 4]   total = 272
same counts? True

            count                                                                        bar
eruptions                                                                                   
[1.5, 2.0)     51                        ###################################################
[2.0, 2.5)     41                                  #########################################
[2.5, 3.0)      5                                                                      #####
[3.0, 3.5)      7                                                                    #######
[3.5, 4.0)     30                                             ##############################
[4.0, 4.5)     73  #########################################################################
[4.5, 5.0)     61              #############################################################
[5.0, 5.5)      4                                                                       ####

centers = (edges[:-1] + edges[1:]) / 2          # midpoint of each bin
width   = np.diff(edges)                        # bin widths (all 0.5 here)

fig, ax = plt.subplots()
ax.bar(centers, counts, width=width * 0.95, color='#4c72b0', edgecolor='white')
ax.set_xlabel('eruption length (minutes)'); ax.set_ylabel('number of eruptions')
ax.set_title('A histogram we drew ourselves: ax.bar over a frequency table')
plt.show()

print('Each bar is one row of the frequency table. ax.hist() just does the')
print('count-then-bar for you in one call — there is no magic underneath.')

Each bar is one row of the frequency table. ax.hist() just does the
count-then-bar for you in one call — there is no magic underneath.

income = loans['annual_income']                 # dollars, self-reported
x = income.to_numpy()

skew_by_hand = np.mean(((x - x.mean()) / x.std(ddof=0)) ** 3)   # the formula, verbatim
skew_scipy   = stats.skew(x)

print(f'skew by hand        : {skew_by_hand:.6f}')
print(f'scipy.stats.skew(x) : {skew_scipy:.6f}')
print('match?', np.isclose(skew_by_hand, skew_scipy))
print(f'\n-> skew = {skew_by_hand:.2f}: a big positive number means a long RIGHT tail.')
print(f'   (mean ${income.mean():,.0f} sits well above median ${income.median():,.0f} — the tail drags the mean up.)')

# A friendly footnote: pandas' Series.skew() applies a small sample-size bias
# correction, so it comes out a hair different. Same idea, slightly different number.
print(f'\npandas income.skew(): {income.skew():.6f}  (bias-corrected, so it differs slightly — that is fine)')

skew by hand        : 9.033136
scipy.stats.skew(x) : 9.033136
match? True

-> skew = 9.03: a big positive number means a long RIGHT tail.
   (mean $79,222 sits well above median $65,000 — the tail drags the mean up.)

pandas income.skew(): 9.034492  (bias-corrected, so it differs slightly — that is fine)

q1_np, q3_np = np.percentile(x, [25, 75])           # NumPy way
q1_pd, q3_pd = income.quantile([0.25, 0.75])        # pandas way
print(f'Q1: np.percentile = ${q1_np:,.0f}   Series.quantile = ${q1_pd:,.0f}   match? {np.isclose(q1_np, q1_pd)}')
print(f'Q3: np.percentile = ${q3_np:,.0f}   Series.quantile = ${q3_pd:,.0f}   match? {np.isclose(q3_np, q3_pd)}')

iqr = q3_np - q1_np
lower_fence = q1_np - 1.5 * iqr
upper_fence = q3_np + 1.5 * iqr
print(f'\nIQR = ${iqr:,.0f}')
print(f'lower fence Q1 - 1.5*IQR = ${lower_fence:,.0f}')
print(f'upper fence Q3 + 1.5*IQR = ${upper_fence:,.0f}')

outlier_mask = (income < lower_fence) | (income > upper_fence)   # the filter-and-count idiom
print(f'\noutliers flagged: {outlier_mask.sum():,} of {len(income):,}  ({outlier_mask.mean()*100:.1f}%)')
print(f'  ... all on the high side: {(income > upper_fence).sum():,} above, {(income < lower_fence).sum():,} below')

Q1: np.percentile = $45,000   Series.quantile = $45,000   match? True
Q3: np.percentile = $95,000   Series.quantile = $95,000   match? True

IQR = $50,000
lower fence Q1 - 1.5*IQR = $-30,000
upper fence Q3 + 1.5*IQR = $170,000

outliers flagged: 508 of 10,000  (5.1%)
  ... all on the high side: 508 above, 0 below

# The loop way (how you might first think of it):
count_loop = 0
for v in income:
    if v < lower_fence or v > upper_fence:
        count_loop += 1

# The vectorized way (what NumPy/pandas is for): a True/False column, summed.
count_vec = ((income < lower_fence) | (income > upper_fence)).sum()

print(f'loop count       : {count_loop:,}')
print(f'vectorized count : {count_vec:,}   <- same number, one line, and far faster on big data')
print('match?', count_loop == count_vec)

loop count       : 508
vectorized count : 508   <- same number, one line, and far faster on big data
match? True

# pd.cut: each eruption labelled with its half-minute interval, then counted.
cut_counts = pd.cut(e, bins=edges, right=False).value_counts().sort_index()
print('pd.cut bin counts :', cut_counts.to_list())

# np.digitize: each eruption gets the index (1..8) of the bin it lands in.
bin_index = np.digitize(e, edges)
idx_counts = [int((bin_index == k).sum()) for k in range(1, len(edges))]
print('np.digitize counts:', idx_counts)
print('same tallies as the frequency table?', idx_counts == cut_counts.to_list())

pd.cut bin counts : [51, 41, 5, 7, 30, 73, 61, 4]
np.digitize counts: [51, 41, 5, 7, 30, 73, 61, 4]
same tallies as the frequency table? True

def flag_outliers(series, k=1.5):
    '''Return (lower_fence, upper_fence, mask) for the k*IQR rule on a numeric Series.

    The mask is a boolean Series: True where the value is past either fence.
    Pass a larger k (e.g. 3.0) for a stricter "far out" flag.
    '''
    s = series.dropna()
    q1, q3 = s.quantile([0.25, 0.75])
    iqr = q3 - q1
    lower = q1 - k * iqr
    upper = q3 + k * iqr
    mask = (series < lower) | (series > upper)
    return lower, upper, mask

# Reuse it on income — should reproduce the Section-3 numbers exactly:
lo, hi, mask = flag_outliers(income)
print(f'income : fence (${lo:,.0f}, ${hi:,.0f})  ->  {mask.sum():,} outliers')

# ... and on a totally different column, no copy-paste of the formula:
lo2, hi2, mask2 = flag_outliers(loans['interest_rate'])
print(f'interest_rate : fence ({lo2:.2f}%, {hi2:.2f}%)  ->  {mask2.sum():,} outliers')

# Dial k up to 3.0 — a stricter "really far out" flag finds fewer:
lo3, hi3, mask3 = flag_outliers(income, k=3.0)
print(f'income (k=3.0) : fence (${lo3:,.0f}, ${hi3:,.0f})  ->  {mask3.sum():,} outliers')

income : fence ($-30,000, $170,000)  ->  508 outliers
interest_rate : fence (1.00%, 23.48%)  ->  364 outliers
income (k=3.0) : fence ($-105,000, $245,000)  ->  187 outliers

Tool	What it is	What we use it for
pandas	tables (`DataFrame`) and columns (`Series`)	loading the CSVs, selecting, quantiles
NumPy	fast math on whole arrays at once	counting bins, the skew formula, the fences
matplotlib	plotting	a histogram we draw ourselves with `ax.bar`
SciPy	a scientific toolbox	one trusted `skew()` to check ours against

The code behind it: histograms, skew & the IQR rule in code¶

1) The goal¶

2) The data in code¶

Selecting a column and a row mask¶

3) Build it step by step¶

A frequency table from scratch¶

A histogram is just bars over that table¶

Skew, from its formula¶

The 1.5 × IQR outlier fence, from the formula¶

4) The idiom — a mask vs. a loop, and binning¶

Counting outliers: vectorized mask vs. a Python loop¶

Binning: which bin does each value fall in?¶

5) Refactor into a reusable function¶

6) Now you code¶

What you learned¶

	rownames	eruptions	waiting
0	1	3.600	79
1	2	1.800	54
2	3	3.333	74
3	4	2.283	62
4	5	4.533	85