""" ============================================================================= CS230 · 第 3 讲工具模块：平面上的二分类数据（"花朵"、同心圆、两个月牙） ============================================================================= 对应讲义：lectures/lecture03_shallow_neural_networks.md（全讲都用这里的数据） ----------------------------------------------------------------------------- 【为什么用二维的玩具数据】 ----------------------------------------------------------------------------- 每个样本只有两个特征 (x₁, x₂)，所以整个数据集可以画在一张纸上（或者终端的一块字符画里）， 网络学到的**决策边界**（decision boundary）也能直接画出来。这让我们能"看见"： · 逻辑回归只能画一条直线； · 一层隐藏层的网络能画出弯曲的边界； · 隐藏单元越多，边界能弯得越复杂（也越容易去追逐噪声）。 高维数据（比如图片）做不到这一点。所以 Coursera 的编程作业 "Planar Data Classification with One Hidden Layer" 也选了二维数据。 ----------------------------------------------------------------------------- 【三个数据集】 ----------------------------------------------------------------------------- flower ：Coursera 作业里的"花朵"数据（load_planar_dataset），8 片花瓣，红蓝两类交替。 这里用 np.random.default_rng 重写了一个**等价**的生成器（公式完全照原作业， 只是随机数的来源不同，所以具体的点和原作业不一样）。 circles ：两个同心圆（外圈一类、内圈一类），仿照原作业 load_extra_datasets 里的 sklearn.datasets.make_circles(factor=0.5, noise=0.3)，用 numpy 重写。 moons ：两个交错的月牙，仿照 sklearn.datasets.make_moons(noise=0.2)，用 numpy 重写。 三者都是"直线分不开"的数据：circles 需要一个圈，moons 需要一条 S 形的曲线， flower 需要把平面切成好几块扇形。 另外还有一个"上帝视角"的工具：因为数据是我们自己生成的，可以算出 Bayes 最优分类器 （bayes_predict）——知道真实分布时最好能做到多准。它给"网络还差多少"提供了一把尺子（讲义 §8）。 【形状约定】和全课一样（第 1 讲 §0.4）：X.shape == (2, m)，每一列是一个样本；Y.shape == (1, m)，值是 0.0 或 1.0。 【注释标签】【思路】【参数】【返回】【边界】【DL 联系】【numpy 小知识】【Python 小知识】 直接运行本文件会打印三个数据集的基本统计： python3 planar_data.py """ import numpy as np # ============================================================================= # 1. 花朵数据 Flower（Coursera: load_planar_dataset） # ============================================================================= def load_planar_dataset(m=400, rng=None, a=4.0): """ 生成 Coursera 作业里的"花朵"数据。 【思路】原作业的生成方法（这里逐行照搬，只把 np.random.randn 换成 rng.normal）： 每一类有 N = m/2 个点。对第 j 类（j = 0, 1）： t = linspace(j·3.12, (j+1)·3.12, N) + 0.2·噪声 （角度：第 0 类占 [0, 3.12]，第 1 类占 [3.12, 6.24]） r = a·sin(4t) + 0.2·噪声 （半径：sin(4t) 让半径随角度起伏 → "花瓣"） 点 = (r·sin t, r·cos t) （极坐标 → 直角坐标） r 可以是负数：负半径相当于把点"翻"到原点的另一侧。正是这一点让两类的花瓣**交替**出现 —— 第 0 类的一部分花瓣被翻到了第 1 类那一半，反之亦然，于是平面被分成很多块红蓝相间的扇形。 两个 0.2 的噪声让花瓣有宽度、让两类在花心附近和花瓣边缘互相重叠， 所以**这个数据集没有 100% 正确的分类器**（任何模型都会错一些点）。 【参数】 m : 样本总数（偶数；两类各一半）。原作业是 400。 rng : np.random.Generator；None 时用种子 1（原作业用 np.random.seed(1)） a : 花瓣的最大半径（原作业是 4） 【返回】(X, Y)：X 形状 (2, m)，Y 形状 (1, m)，前 m/2 个样本是第 0 类，后 m/2 个是第 1 类 【边界】m 是奇数时，每类 m // 2 个点，总数会比 m 少 1。 【DL 联系】第 3 讲用它说明"逻辑回归画不出这种边界，一层隐藏层可以"（§1、§8）。 【numpy 小知识】 np.linspace(start, stop, N)：在 [start, stop] 之间均匀取 N 个点（包含两端）。 rng.normal(size=N)：N 个标准正态分布 N(0, 1) 的随机数（等价于老写法 np.random.randn(N)）。 np.zeros((2, m))：2 行 m 列的全 0 数组。注意形状要写成一个元组 (2, m)。 X[:, ix] = ...：把第 ix 这些**列**整块赋值。 """ if rng is None: rng = np.random.default_rng(1) N = m // 2 X = np.zeros((2, 2 * N)) Y = np.zeros((1, 2 * N)) for j in range(2): ix = slice(N * j, N * (j + 1)) # 第 j 类占的列：[N·j, N·(j+1)) t = np.linspace(j * 3.12, (j + 1) * 3.12, N) + rng.normal(size=N) * 0.2 # 角度 theta r = a * np.sin(4 * t) + rng.normal(size=N) * 0.2 # 半径（可以为负） X[:, ix] = np.vstack([r * np.sin(t), r * np.cos(t)]) # 两行：x₁ 和 x₂ Y[:, ix] = j return X, Y # ============================================================================= # 2. 两个对照数据集：同心圆、两个月牙（仿照 sklearn 的同名函数） # ============================================================================= def make_circles(m=400, rng=None, factor=0.5, noise=0.3): """ 两个同心圆：外圈（半径 1）是第 0 类，内圈（半径 factor）是第 1 类，每个点加上高斯噪声。 【思路】仿照 sklearn.datasets.make_circles：角度在 [0, 2π) 上均匀取点（不含 2π，避免首尾重合）， 再给两个坐标各加 noise 标准差的高斯噪声。Coursera 的 load_extra_datasets 用的参数是 factor = 0.5、noise = 0.3：噪声相对于两圈之间 0.5 的间距**很大**，两类有不少重叠， 所以这个数据集的最好准确率明显低于 100%。 【参数】m：样本总数（两类各 m // 2）；rng：随机数生成器（None 时用种子 0）；factor：内圈半径；noise：噪声标准差 【返回】(X, Y)，形状 (2, m)、(1, m) 【DL 联系】决策边界需要是一个"圈" —— 一条直线永远做不到。 【numpy 小知识】np.linspace(0, 2π, N, endpoint=False) 不包含终点 2π。 np.hstack([A, B]) 把两个数组左右拼接（列数相加）。 """ if rng is None: rng = np.random.default_rng(0) N = m // 2 ang = np.linspace(0.0, 2.0 * np.pi, N, endpoint=False) outer = np.vstack([np.cos(ang), np.sin(ang)]) # (2, N) 外圈，第 0 类 inner = factor * np.vstack([np.cos(ang), np.sin(ang)]) # (2, N) 内圈，第 1 类 X = np.hstack([outer, inner]) + rng.normal(scale=noise, size=(2, 2 * N)) Y = np.hstack([np.zeros((1, N)), np.ones((1, N))]) return X, Y def make_moons(m=400, rng=None, noise=0.2): """ 两个交错的月牙：上面的半圆是第 0 类，下面错开的半圆是第 1 类。 【思路】仿照 sklearn.datasets.make_moons： 第 0 类：(cos t, sin t)，t ∈ [0, π] （上半圆，圆心 (0, 0)） 第 1 类：(1 − cos t, 1 − sin t − 0.5)，t ∈ [0, π] （下半圆，圆心 (1, 0.5)，开口朝上） 两个月牙互相"咬"在一起，再加 noise 的高斯噪声。Coursera 用的是 noise = 0.2。 【参数】m：样本总数；rng：随机数生成器（None 时用种子 0）；noise：噪声标准差 【返回】(X, Y)，形状 (2, m)、(1, m) 【DL 联系】一条直线能分对大部分点（月牙的"主体"在两侧），但两个月牙的尖端互相伸进对方的地盘， 要把尖端也分对，边界必须弯成 S 形。所以它是"逻辑回归不差、但网络更好"的例子。 """ if rng is None: rng = np.random.default_rng(0) N = m // 2 t = np.linspace(0.0, np.pi, N) upper = np.vstack([np.cos(t), np.sin(t)]) lower = np.vstack([1.0 - np.cos(t), 1.0 - np.sin(t) - 0.5]) X = np.hstack([upper, lower]) + rng.normal(scale=noise, size=(2, 2 * N)) Y = np.hstack([np.zeros((1, N)), np.ones((1, N))]) return X, Y # ============================================================================= # 3. 统一的入口：训练集 + 独立生成的测试集 # ============================================================================= DATASETS = { "flower": load_planar_dataset, "circles": make_circles, "moons": make_moons, } # 每个数据集训练集 / 测试集用的种子。训练集和测试集用**不同的种子**独立生成： # 两者来自同一个分布，但没有一个点是共享的 —— 测试准确率才是"没见过的数据上的表现"。 # flower 的训练种子 1 是向原作业的 np.random.seed(1) 致敬（生成器不同，点并不相同）。 SEEDS = { "flower": (1, 101), "circles": (2, 102), "moons": (3, 103), } def get_dataset(name="flower", m_train=400, m_test=2000): """ 返回某个数据集的训练集和测试集。 【参数】name："flower" | "circles" | "moons"；m_train、m_test：训练 / 测试样本数 训练集默认 400 个（和原作业一样）；测试集默认 2000 个 —— 测试集越大，测出来的准确率越稳 （400 个测试样本时，准确率本身的随机误差约 ±1.5 个百分点；2000 个时约 ±0.7）。 【返回】(X_train, Y_train, X_test, Y_test)，形状 (2, m_train)、(1, m_train)、(2, m_test)、(1, m_test) 【边界】name 写错时抛出 KeyError（并列出可用的名字），而不是悄悄返回别的数据。 【Python 小知识】DATASETS[name] 取出的是一个**函数**（不带括号），后面加 (...) 才是调用它。 把函数放进字典，是一种"按名字选函数"的常用写法。 """ if name not in DATASETS: raise KeyError(f"unknown dataset {name!r}; choose from {sorted(DATASETS)}") make = DATASETS[name] seed_tr, seed_te = SEEDS[name] X_tr, Y_tr = make(m_train, np.random.default_rng(seed_tr)) X_te, Y_te = make(m_test, np.random.default_rng(seed_te)) return X_tr, Y_tr, X_te, Y_te # ============================================================================= # 4. "上限"：知道真实的生成过程时，最好能做到多准（Bayes 最优分类器） # ============================================================================= def _normal_cdf(x): """ 标准正态分布的累积分布函数 Φ(x) = P(N(0, 1) ≤ x)，逐元素计算。 【numpy 小知识】numpy 没有 erf 函数；np.vectorize(math.erf) 把只能处理一个数的 math.erf 包装成能逐元素处理数组的函数（本质上还是 Python 循环，慢，但这里数据量不大）。 """ import math return 0.5 * (1.0 + np.vectorize(math.erf)(np.asarray(x) / np.sqrt(2.0))) def class_likelihoods(name, X): """ 在**真实的生成过程**下，每个点属于第 0 类、第 1 类的（相对）概率密度 p(x | 类别)。 【思路】我们自己写了数据生成器，所以知道"上帝视角"的真相。两类各占一半时，Bayes 最优分类器就是 "哪一类的 p(x | 类别) 大就猜哪一类"——任何模型（包括无限大的网络）的期望准确率都不可能超过它。 · circles / moons：每个点 = 曲线上的一个位置 + 各向同性的高斯噪声。p(x | 类别) = 沿曲线 把高斯密度平均一遍（用很密的网格做数值积分）。 · flower：点 = (r sin t, r cos t)，t = 均匀分布 + N(0, 0.2²)，r = 4 sin 4t + N(0, 0.2²)。 给定 x，能生成它的 (t, r) 只有 t = θ + kπ、r = (−1)^k·|x| 这几组（θ 是 x 的极角）， 所以 p(x | 类别 j) ∝ Σ_k p_j(t = θ + kπ) · N((−1)^k |x|; 4 sin 4θ, 0.2²)。 （变量替换的 Jacobian 是 |r| = |x|，两类相同，比较大小时可以约掉。） 这里用"连续"的版本近似原作业里 linspace 取出的 N 个等间隔的 t（N 很大时两者几乎一样）。 【参数】name：数据集名字；X：(2, m) 【返回】(2, m) 的数组：第 j 行是 p(x | 类别 j)（只差一个两类共同的常数因子） 【边界】只对本模块默认参数（a = 4、factor = 0.5、noise = 0.3 / 0.2）生成的数据成立。 """ if name == "circles": ang = np.linspace(0.0, 2.0 * np.pi, 720, endpoint=False) curves = [np.vstack([np.cos(ang), np.sin(ang)]), 0.5 * np.vstack([np.cos(ang), np.sin(ang)])] sigma = 0.3 elif name == "moons": t = np.linspace(0.0, np.pi, 1000) curves = [np.vstack([np.cos(t), np.sin(t)]), np.vstack([1.0 - np.cos(t), 1.0 - np.sin(t) - 0.5])] sigma = 0.2 elif name == "flower": theta = np.arctan2(X[0], X[1]) # x₁ = r sin t、x₂ = r cos t ⇒ t = atan2(x₁, x₂) rho = np.hypot(X[0], X[1]) out = np.zeros((2, X.shape[1])) for j in range(2): a, b = 3.12 * j, 3.12 * (j + 1) for k in range(-1, 4): # t = θ + kπ 覆盖了 t 可能的取值范围 t = theta + k * np.pi p_t = (_normal_cdf((t - a) / 0.2) - _normal_cdf((t - b) / 0.2)) / (b - a) r = rho if k % 2 == 0 else -rho p_r = np.exp(-(r - 4.0 * np.sin(4.0 * theta)) ** 2 / (2 * 0.2 ** 2)) out[j] += p_t * p_r return out else: raise KeyError(f"unknown dataset {name!r}") out = np.zeros((2, X.shape[1])) for j, C in enumerate(curves): d2 = (X[0][:, None] - C[0][None, :]) ** 2 + (X[1][:, None] - C[1][None, :]) ** 2 # (m, K) 距离平方 out[j] = np.mean(np.exp(-d2 / (2 * sigma ** 2)), axis=1) return out def bayes_predict(name, X): """ Bayes 最优分类器的预测：p(x | 第 1 类) > p(x | 第 0 类) 就猜第 1 类（两类先验各一半）。 【参数】name：数据集名字；X：(2, m) 【返回】(1, m) 的 0.0 / 1.0 数组（和 shallow_nn.predict 的格式一样，可以直接和 Y 比较） """ L = class_likelihoods(name, X) return (L[1] > L[0]).astype(float)[None, :] def bayes_accuracy(name, m=20000, seed=999): """ 用一个很大的新样本估计 Bayes 最优分类器的准确率 —— 这个数据集"能做到的上限"。 【参数】name：数据集名字；m：样本数（越大越准，20000 个时随机误差约 ±0.2 个百分点）；seed：种子 【返回】float 【DL 联系】第 10 讲的 Bayes error（贝叶斯误差）：任何模型都消不掉的那部分错误。 真实任务里不知道它，常用"人类水平"来近似；我们的玩具世界里可以直接算出来。 """ X, Y = DATASETS[name](m, np.random.default_rng(seed)) return float(np.mean(bayes_predict(name, X) == Y)) if __name__ == "__main__": for name in DATASETS: print(f"{name:<8} Bayes 最优准确率（上限）≈ {bayes_accuracy(name):.3f}") for name in DATASETS: X_tr, Y_tr, X_te, Y_te = get_dataset(name) print(f"{name:<8} X_train {X_tr.shape} Y_train {Y_tr.shape} X_test {X_te.shape} " f"正类比例 {Y_tr.mean():.2f} x₁ 范围 [{X_tr[0].min():.2f}, {X_tr[0].max():.2f}] " f"x₂ 范围 [{X_tr[1].min():.2f}, {X_tr[1].max():.2f}]")