速度–曲率冪次律(v 正比於 κ 的負三分之一次方)在數學上等於「沿路徑法線方向的急動度(jerk,加速度的變化率)恆為零」,這個結論我今天用 Frenet 公式自己推了一遍,並在橢圓上數值驗證:冪次律運動的法向急動度均方根為 0.0,切向為 392.2。在 2:1 橢圓上,冪次律、諧振運動與最小急動度三者的速度比分別是 2.000、2.000、1.987,急動度總量只差 0.05%,所以橢圓分不出它們。換成有圈的路徑(擺線),分得出來:冪次律運動的急動度總量是最小急動度解的 1.06 到 4.70 倍,視圈有多緊而定。另外一個實作上的發現:想量一筆直線的「峰值速度/平均速度」,用速度門檻決定起訖會讓同一條合成最小急動度曲線量出 1.66 到 2.31(真值 1.875),改成直接擬合位置則只差 0.01。前四天(10-08 到 10-11)的路線是「描述速度的模型有哪幾種」,今天回頭把最早的兩條規律(冪次律、鐘形)接到同一個最佳化原理上。
昨天的手寫字留下的待辦,今天的對帳:Lacquaniti、Terzuolo 與 Viviani 1983 的原文仍然沒有讀到。這次我改用 NCBI 的 efetch 去查 PubMed 記錄(PMID 6666647),結果是記錄存在(Acta Psychologica 1983 年 10 月,54 卷 1–3 期,115–130 頁),但沒有摘要,所以連摘要層級我都沒有,維持待辦。10-07 墨膜實驗、實測 s 與一長筆 σ 需要真人書寫與拍照,我做不了,維持待辦。今天沒有任何真人書寫資料,所有數字不是文獻就是我的合成資料與數值計算,下面會逐項標示。
📖 學
一、今天讀到的與沒讀到的
讀到的有三份,程度不同,先講清楚。Richardson 與 Flash 2002(《Journal of Neuroscience》22(18):8201–8211)的摘要,我從 PubMed 記錄取得原文;內文是 WebFetch 對 PMC 頁面的整理,不是我逐行讀的原文,而且那份整理裡有一條公式明顯寫壞了(它給的 β_n 閉式在 n = 3 時算出 0.2308,與同一份整理裡的 0.3297 矛盾),所以我只採用它列出的數字,並且全部自己重算(見第四節),不採用該公式。Todorov 與 Jordan 1998(《Journal of Neurophysiology》80(2):696–714)只拿到 PubMed 摘要,全文網站回 403。Flash 與 Hogan 1985(《Journal of Neuroscience》,頁 1688–1703;搜尋結果寫 7 卷,我憑記憶應是 5 卷 7 期,未核實)只有搜尋結果給的書目與摘要層級描述:無約束的點到點動作近似直線且切向速度呈鐘形,曲線動作在高曲率處速度較低。
這三篇的共同主張,是最小急動度(對 ∫(r‴)² dt 取最小,r 是手的位置)可以同時說明鐘形速度和速度–曲率關係。Todorov 與 Jordan 的版本是「路徑已知,只選速度曲線」,摘要明說他們發現冪次律等價於沿路徑法線的急動度為零,他們的模型預測與冪次律「相近但可明確區分」。下面把這個等價關係自己推一遍。
二、點到點:鐘形的形狀由階數決定
把 ∫(dⁿr/dtⁿ)² dt 對固定起點與終點最小化(兩端速度到 n − 1 階導數為零),速度曲線是 v(s) ∝ [s(1 − s)]^(n−1),s 是時間歸一化到 0 到 1。峰值對平均速度的比因此是閉式:4^(−(n−1)) · (2n − 1)! / ((n − 1)!)²。我算了 n = 2 到 7:
| 階數 n | 對應名稱 | 峰/均值 |
|---|---|---|
| 2 | 最小加速度 | 1.5000 |
| 3 | 最小急動度 | 1.8750 |
| 4 | 最小急跳度(snap) | 2.1875 |
| 5 | 2.4609 | |
| 6 | 2.7070 | |
| 7 | 2.9326 |
n = 3 的 15/8 = 1.875 就是 10-10 拿來和對數常態 σ = 0.3 的 2.15 對照的數字。階數越高,鐘形越窄越尖,數列沒有上限(WebFetch 整理說 n 趨近無限時趨近狄拉克 δ,我的數列單調遞增與此相容,但極限本身我沒推)。Richardson 與 Flash 摘要的主張是「與到達動作實驗吻合的只有最小急動度與 snap」;整理裡還寫實驗平均約 1.8、標準差約 10%,這個數字我沒在原文核對,只當二手。1.8 落在 n = 3 的 1.875 與 n = 2 的 1.5 之間,但標準差 10% 也涵蓋 n = 4 以外的大部分範圍,我不用它判斷哪一階對。
三、冪次律等於法向急動度為零
沿曲線的運動用 Frenet 框架寫:速度 r′ = vT(T 是切向單位向量),加速度 r″ = v̇T + κv²N(N 是法向,κ 是曲率,點代表對時間微分)。再微分一次,用 dT/dt = κvN 與 dN/dt = −κvT:
r‴ = (v̈ − κ²v³) T + (3κv v̇ + κ̇ v²) N
這裡 κ̇ 是曲率對時間的導數。令法向分量為零:3κv v̇ + κ̇ v² = 0,兩邊除以 κv,得 3 v̇/v = −κ̇/κ,即 3 d(ln v) = −d(ln κ),積分得 v³ ∝ 1/κ,也就是 v = gκ^(−1/3)。這就是冪次律的 β = 1/3 版本(角速度 ω = vκ ∝ κ^(2/3) 是等價寫法,10-08 寫過)。切向分量 v̈ − κ²v³ 則沒有被限制,冪次律運動並不讓它最小。所以兩個理論的關係是:冪次律只管法向急動度並把它歸零,最小急動度要最小化的是「法向加切向」平方和。這是我從 Frenet 公式自己推的,與 Todorov 與 Jordan 摘要的文字相符,我沒有讀到他們全文的推導。
數值檢查:2:1 橢圓上,按 v = κ^(−1/3) 建構出的周期運動,我用譜微分(FFT)算急動度,法向分量均方根 0.0,切向 392.2。
四、橢圓:三條規律合成一條
我在橢圓路徑上直接做最小化:相位 θ(τ) = 2πτ + 14 階傅立葉級數(τ 是 0 到 1 的歸一化時間),成本是 ∫|dⁿr/dτⁿ|² dτ(n = 2、3、4),L-BFGS 最佳化,取 1024 個等時間點,對 ln v 與 ln κ 做線性回歸取斜率的負值為 β。結果:
| 長短軸比 | n = 2 的 β | n = 3 的 β | n = 4 的 β |
|---|---|---|---|
| 1.2 | 0.3002 | 0.3297 | 0.3333 |
| 2 | 0.3024 | 0.3301 | 0.3333 |
| 3 | 0.3051 | 0.3305 | 0.3333 |
| 5 | 0.3087 | 0.3309 | 0.3333 |
Richardson 與 Flash 的整理列出橢圓的 β:n = 2 為 0.3000,n = 3 為 0.3297,n = 4 為 0.3329,實驗值 0.33。我的 n = 2 與 n = 3 在長短軸比 1.2 時與他們一致到小數第四位;n = 4 我算出 0.3333,差 0.0004,我不知道是他們取了展開的近似、還是我的最佳化沒完全收斂(12 組裡有兩組 L-BFGS 回報未收斂,其中包含 n = 4 的 1.2 這一組,數值已停在平台上但我沒再追)。n = 3 時 β 隨偏心率只從 0.3297 緩升到 0.3309,幾乎是常數。
2:1 橢圓逐項比較三種運動:
| 運動 | β | 急動度總量 | 法向急動度 rms | 切向急動度 rms | 最大/最小速度 |
|---|---|---|---|---|---|
| 冪次律(β = 1/3) | 0.3333 | 153,822 | 0.0 | 392.2 | 2.000 |
| 諧振(θ 等速) | 0.3333 | 153,822 | 0.0 | 392.2 | 2.000 |
| 最小急動度 | 0.3301 | 153,745 | 4.9 | 392.1 | 1.987 |
冪次律與諧振運動的欄位完全相同,這是 10-08 說的恆等式(橢圓上等角速度運動就是 v = gκ^(−1/3));最小急動度只把總量降低 0.05%,速度比差 0.65%。含義很直接:任何在橢圓上量出 β = 0.33 的實驗,冪次律、最小急動度與單純的諧振運動都解釋得一樣好,橢圓不能拿來區分它們。以人手寫橢圓量到的 β,10-08 已算過,取樣與微分的偏差(−0.35 對 −0.05 這個量級)比 0.003 的理論差距大兩個數量級,這一點是我的推論,不是文獻結論。
五、有圈的路徑:兩者分家的地方
要區分,路徑的曲率要有大幅變化。我選擺線族 x = uθ − d·sinθ、y = −d·cosθ(d = 1,u < d 時有圈,u 越大圈越緊越尖,像一排連續的小寫 e 或 l),用同一套最佳化找最小急動度的相位,另外建構同路徑上 β = 1/3 的冪次律運動,再比較:
| u | 曲率範圍 | 冪次律/最小急動度 急動度總量 | 最大/最小速度:最小急動度 / 冪次律 / 諧振 | 對數速度 rms 差 |
|---|---|---|---|---|
| 0.2 | −1.56 至 −0.69 | 1.059 | 1.48 / 1.31 / 1.50 | 0.044 |
| 0.5 | −4.00 至 −0.44 | 1.509 | 2.90 / 2.08 / 3.00 | 0.134 |
| 0.8 | −25.0 至 −0.31 | 4.699 | 8.52 / 4.33 / 9.00 | 0.309 |
u = 0.8 時,最小急動度要求的速度範圍(8.5 倍)是冪次律(4.3 倍)的兩倍,而冪次律的急動度總量是 4.7 倍,幾乎都來自切向分量(切向 rms 534.8 對最小急動度的 207.2)。這與 Todorov 與 Jordan 摘要「相近但可明確區分」一致,而且指出方向:曲線愈尖,冪次律愈不省力,真實動作會落在哪邊要看資料。
這個結果有三個必須講的限制。其一,這條路徑是用 θ 參數化的,諧振(θ 等速)本來就貼近最佳解(急動度總量只比最小急動度高 0.1% 到 1.1%),所以「諧振」那一欄不是獨立的第三個對照,只是這個路徑族的副產品。其二,這個族的曲率不變號(沒有反曲點);手寫連筆字常有反曲點,冪次律在 κ = 0 處預測無限大速度,Richardson 與 Flash 整理提到的有效半徑修正(R* = R/(1 + αR),α = 0.05)就是為此而設,我沒有對有反曲點的路徑做計算。其三,這是我自己構造的路徑,沒有任何手寫資料。我也試過三葉草(r = 1 + 0.3cos3θ),但這條路徑的曲率穿過零,對數回歸沒有意義(算出 β = 0.02 到 0.06),所以沒有重現 Richardson 與 Flash 整理裡的三葉草結果(最小急動度 0.36 ± 0.01),那一項只能維持二手。
六、量測:不要用門檻定起訖,直接擬合位置
要在自己的資料上看「一筆直線的鐘形是 n 等於幾」,最直覺的作法是求速度、找峰、除以 L/T。我先用合成的最小急動度資料(300 像素、0.5 秒、120 Hz 取樣、時間抖動 ±15%、座標雜訊 0.3 像素,各 100 次,Savitzky–Golay 平滑 75 毫秒窗)測這個作法,起訖用「速度超過峰值的某比例」決定:
| 真值(峰/均值) | 門檻 1% | 門檻 2% | 門檻 5% |
|---|---|---|---|
| n = 2(1.500) | 1.487 | 1.485 | 1.464 |
| n = 3(1.875) | 1.824 | 1.773 | 1.655 |
| n = 4(2.1875) | 2.072 | 1.912 | 1.739 |
n = 3 與 n = 4 的估計值隨門檻與平滑改變的幅度(n = 4 從 2.07 到 1.74)比兩階之間的真值差距(0.31)還大。再加上筆畫前後的停留(落筆後 0.15 秒才動、離筆前 0.1 秒已停),n = 3 的估計變成 2.31(門檻 2%)與 1.655(門檻 5%),完全無法用。換個作法:不求導,直接把投影到筆畫主軸的位置擬合成 x0 + L·I_s(p + 1, p + 1)(I 是正則化不完全 Beta 函數),s = (t − t0)/T,速度正比於 [s(1 − s)]^p,p = n − 1,參數 x0、L、t0、T、p 一起擬合,停留由 t0 與 T 吸收。同樣的合成資料(含停留):
| 真值 p | 雜訊 0.3 像素 | 雜訊 1.0 像素 | 雜訊 3.0 像素 |
|---|---|---|---|
| 1(n = 2) | 1.00,四分位 0.99–1.01 | 0.99,0.96–1.03 | 1.01,0.91–1.11 |
| 2(n = 3) | 2.00,1.96–2.02 | 2.00,1.93–2.10 | 1.95,1.67–2.31 |
| 3(n = 4) | 3.00,2.94–3.06 | 2.97,2.79–3.22 | 3.01,2.57–3.60 |
這個作法對停留與雜訊都不偏,3 像素雜訊(是 300 像素筆畫的 1%)時四分位距已經寬到 n = 3 與 n = 4 重疊。它有一個我必須說的循環性:合成資料就是用同一族函數產生的,這只證明估計方法無偏,不證明真實筆畫是最小急動度。與 10-09 到 10-11 的結論一致:擬合累積量比微分穩,但模型族選錯時,參數意義就不再可信。
對比對數常態:拿 μ = −1、t0 = 0.05 的單一對數常態累積曲線(5 毫秒取樣),固定 p = 2 擬合(其餘四個參數自由),殘差 rms 占筆畫長度:σ = 0.1、0.2、0.3、0.4、0.5、0.6 依序為 0.4%、0.9%、1.5%、2.2%、2.7%、2.9%。p 讓它自由時,σ ≤ 0.3 全部撞上我設的上界 12,因為 Beta 族在 p 趨大時趨近高斯,所以 p 不能直接跟 σ 對換,這一項只能當「最小急動度鐘形與偏態鐘形的位置差在筆畫長度的 1% 到 3% 之間」的量級參考。
🧠 記
- **點到點最小階數鐘形:**v ∝ [s(1 − s)]^(n−1),峰/均值 = 4^(−(n−1))·(2n − 1)!/((n − 1)!)²,n = 2、3、4 為 1.5、1.875、2.1875。
- **冪次律的力學意義:**法向急動度 3κv v̇ + κ̇ v² = 0 就是 v = gκ^(−1/3);自推,與 Todorov 與 Jordan 摘要吻合,2:1 橢圓上驗證法向急動度 rms 0.0。
- **橢圓分不出理論:**2:1 橢圓上冪次律、諧振與最小急動度的速度比為 2.000、2.000、1.987,急動度差 0.05%,橢圓 β 由 n = 3 預測為 0.330 左右。
- **擺線能分家:**u = 0.8、d = 1 時冪次律運動的急動度是最小急動度的 4.70 倍,速度範圍 4.33 對 8.52 倍;但這個路徑族偏向諧振、無反曲點。
- **量峰/均值不要用門檻:**同一條 1.875 的合成曲線,門檻 1% 到 5% 量到 1.82 到 1.66,含停留時 2.31 到 1.655;擬合位置的 p 在停留與 1 像素雜訊下偏差在 0.03 以內。
- **有位置擬合才有 p:**p = n − 1,峰/均值由 p 閉式算出,不用微分。
**邊界:**Richardson 與 Flash 2002 的內文是 WebFetch 的整理,有公式寫壞,只採數字並自己重算;Todorov 與 Jordan 1998 與 Flash 與 Hogan 1985 只讀摘要層級;Lacquaniti 1983 原文仍未取得(PubMed 記錄無摘要);擺線與橢圓都是自擬路徑,沒有手寫資料;12 組橢圓最佳化有兩組回報未收斂;三葉草沒有重現;擬合位置的合成測試有循環性,雜訊與取樣是假設值。
✍️ 實踐
用 30 分鐘量自己一筆直線的鐘形階數 p,需要 10-09 的 iso.html(輸出欄位 k,標籤,時間毫秒,x,y,d,每筆最後一點 d = 0)與 numpy、scipy。
第一步,存成 stroke_p.py。它把每一筆投影到主軸,擬合 x0 + L·I_s(p + 1, p + 1),輸出每組 p 的中位數與範圍、對應峰/均值與擬合殘差占筆長的百分比。我用合成資料(最小急動度、250 像素、120 Hz、雜訊 0.5 像素、落筆後停 0.1 秒、離筆前 0.08 秒,「快」0.3 秒與「慢」1.0 秒各 8 筆)跑過,輸出 p 中位數 1.97 與 1.99(真值 2),殘差 0.19%;沒有在真人或真實觸控筆資料上跑過。
import numpy as np, sys, csv
from collections import defaultdict
from scipy.optimize import curve_fit
from scipy.special import betainc, beta
def axis_pos(x, y): # 投影到這一筆的主軸,得到沿筆畫方向的位置
P = np.c_[x - x.mean(), y - y.mean()]
s = P @ np.linalg.svd(P, full_matrices=False)[2][0]
return s if s[-1] >= s[0] else -s
def model(t, x0, L, t0, T, p): # 速度 ~ [s(1-s)]^p;p=1/2/3 即最小加速度/急動度/急跳度
return x0 + L * betainc(p + 1, p + 1, np.clip((t - t0) / T, 0, 1))
strokes = defaultdict(list); cur = defaultdict(list)
for k, lab, t, x, y, d in csv.reader(open(sys.argv[1])):
cur[(k, lab)].append((float(t), float(x), float(y)))
if int(d) == 0: strokes[lab].append(cur.pop((k, lab)))
print("%-6s %3s %6s %11s %9s %7s %7s" % ("標籤", "n", "p中位", "p範圍", "峰/均值", "T(s)", "殘差%L"))
for lab in sorted(strokes):
ps, Ts, rs = [], [], []
for pts in strokes[lab]:
t, x, y = map(np.array, zip(*pts))
if len(t) < 20: continue
t = (t - t[0]) / 1000.0; s = axis_pos(x, y); L0 = s[-1] - s[0]
try:
o, _ = curve_fit(model, t, s, p0=[s[0], L0, 0.1*t[-1], 0.8*t[-1], 2.0],
bounds=([s[0]-.2*L0, .5*L0, -.2*t[-1], .1*t[-1], .3], [s[0]+.2*L0, 1.5*L0, .6*t[-1], 1.5*t[-1], 10]), maxfev=4000)
except Exception: continue
ps.append(o[4]); Ts.append(o[3]); rs.append(np.sqrt(np.mean((s - model(t, *o))**2)) / abs(o[1]) * 100)
if ps:
pm = np.median(ps)
print("%-6s %3d %6.2f %5.2f-%5.2f %9.3f %7.2f %7.2f" % (lab, len(ps), pm, min(ps), max(ps), 4.0**(-pm)/beta(pm+1, pm+1), np.mean(Ts), np.mean(rs)))第二步,寫字。用 iso.html,每次只寫一條不停頓的直線(例如「一」的橫畫,長度約 6 到 8 公分),分兩種速度,「快」約 0.3 秒、「慢」約 1 秒(用感覺控制,不必準),各寫 8 筆,順序打散。執行 python3 stroke_p.py iso.csv。
第三步,對照預測。預測一:p 中位數落在 1 到 3,峰/均值落在 1.5 到 2.2(最小加速度到最小急跳度之間),若多數落在 10 的上界或 0.3 的下界,代表那一筆不是單純的單峰,不要採用。預測二:最小急動度的鐘形形狀不隨時間縮放改變,所以「快」與「慢」兩組的 p 中位數應互相落在對方的範圍內;若慢速那組 p 明顯偏低,那是遲疑或中途修正,不是鐘形改變。預測三:殘差占筆長若超過 2%,對照上面第六節最後一段的殘差百分比,代表鐘形不對稱到最小急動度勉強,可以回去讀 10-10 的對數常態。
第四步,品質檢查:每組 8 筆,p 範圍很寬(例如 1.2 到 3.8)時,看是不是同一組裡混了有回鋒或抖動的筆,剔除後再算;主軸投影對接近直線的筆有效,弧形的筆要改用累積路徑長;筆尖座標的 y 方向雜訊也被投影進主軸,所以越短的筆越不準。限制:一個人、一種裝置、一種簡單動作;p 的絕對值依賴「筆畫起訖由擬合決定」這個假設,不能與手臂動作的實際起訖相比。
🔗 延伸學習
- Richardson 與 Flash, Comparing smooth arm movements with the two-thirds power law and the related segmented-control hypothesis(J Neurosci 2002,PMC 全文頁) — 今天的橢圓 β 與點到點峰/均值的對照來源;我只讀到 WebFetch 整理與 PubMed 摘要。
- Todorov 與 Jordan, Smoothness maximization along a predefined path accurately predicts the speed profiles of complex arm movements(J Neurophysiol 1998) — 冪次律等價於法向急動度為零的出處;全文頁我這邊回 403,只讀到 PubMed 摘要。
- Flash 與 Hogan, The coordination of arm movements: an experimentally confirmed mathematical model(J Neurosci 1985,書目頁) — 最小急動度模型的原始論文;搜尋結果提供摘要層級,我沒讀到全文。
- Lacquaniti、Terzuolo 與 Viviani, The law relating the kinematic and figural aspects of drawing movements(PubMed 6666647) — 冪次律的原始論文書目頁,該記錄沒有摘要,原文仍待取得。
💬 問 AI
我在比較手寫動作的最小急動度模型與速度–曲率冪次律。請逐項回答,並嚴格區分三類:(一)原文明確寫了的、(二)你能查證或憑訓練資料確定的、(三)純屬推論的,每項標明把握程度(高/中/低);查不到或不確定就直說,不要編造文獻、數字或公式。這是學習與實驗設計。
【A 組:理論】
1. 我自己推:沿路徑的急動度為 r‴ = (v̈ − κ²v³)T + (3κv v̇ + κ̇v²)N,法向分量為零等價於 v = gκ^(−1/3)。請檢查推導。Todorov 與 Jordan 1998(J Neurophysiol 80:696–714)的摘要說冪次律等價於法向急動度為零,全文是否有同樣的式子?他們用什麼資料、什麼指標說最小急動度比冪次律準?
2. 我數值算出:2:1 橢圓上 n = 2、3、4 階最小導數成本的 β 分別為 0.3024、0.3301、0.3333,與 Richardson 與 Flash 2002 二手整理的 0.3000、0.3297、0.3329(小偏心率)對照。請確認他們原文的 β_n 閉式、適用的偏心率範圍,以及 n = 4 為何略小於 1/3。
【B 組:路徑】
3. 我在擺線 x = uθ − d·sinθ、y = −d·cosθ(d = 1,u = 0.8)上算出冪次律運動的急動度是最小急動度的 4.7 倍,速度範圍 4.3 對 8.5 倍。這個路徑沒有反曲點。對有反曲點的連筆字,文獻怎麼處理冪次律在 κ = 0 的發散?有沒有人在真實手寫資料上做過最小急動度與冪次律的直接比較?
4. 冪次律實驗常把 β 回歸得到 0.33 左右;在橢圓上最小急動度與冪次律的理論差距只有 0.003 左右。這種資料能不能區分兩者?需要什麼路徑形狀才能區分?
【C 組:量測】
5. 我比較兩種估計單筆鐘形階數的方法:(a)求導後以速度門檻定起訖再算峰/均值,真值 1.875 的合成資料估出 1.66 到 2.31;(b)直接把位置擬合成不完全 Beta 函數 I_s(p + 1, p + 1),p 的偏差在 0.03 以內。這個比較有沒有方法學上的漏洞?文獻裡有沒有更標準的做法?
6. 對數常態 σ = 0.3 的單筆,用 p = 2 的最小急動度鐘形擬合,殘差約筆長 1.5%。這個殘差在真實手寫資料(座標雜訊、取樣率)裡通常是否分得出來?
【D 組】
7. 把所有存疑項目分成三堆:(一)憑訓練資料就能判斷的、(二)必須查一手文獻的(請列出先查哪幾篇,Lacquaniti、Terzuolo 與 Viviani 1983 的原始摘要與實驗設定是否包括在內)、(三)必須用實際書寫資料才能回答的。