Import project sources

.gitignore

@@ -0,0 +1,15 @@
+
+# Local files and caches
+.DS_Store
+__pycache__/
+*.pyc
+.env
+.env.*
+!.env.example
+*.zip
+/analysis_data.json
+/*_results.json
+/*.log
+/*.png
+/内容统计分析报告.md
+*.xlsx

analyze.py

@@ -0,0 +1,136 @@
+import json
+import math
+import warnings
+from pathlib import Path
+
+import numpy as np
+import openpyxl
+import pandas as pd
+
+warnings.filterwarnings('ignore', message='Workbook contains no default style')
+OUT = Path(__file__).parent
+SOURCE = '/Users/zzyang/Downloads/笔记列表明细表.xlsx'
+
+
+def beta_cf(a, b, x):
+    qab, qap, qam = a + b, a + 1, a - 1
+    c, d = 1.0, 1 - qab * x / qap
+    d = 1 / (d if abs(d) > 1e-300 else 1e-300)
+    h = d
+    for m in range(1, 301):
+        m2 = 2 * m
+        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
+        d = 1 + aa * d
+        c = 1 + aa / c
+        d = 1 / (d if abs(d) > 1e-300 else 1e-300)
+        c = c if abs(c) > 1e-300 else 1e-300
+        h *= d * c
+        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
+        d = 1 + aa * d
+        c = 1 + aa / c
+        d = 1 / (d if abs(d) > 1e-300 else 1e-300)
+        c = c if abs(c) > 1e-300 else 1e-300
+        change = d * c
+        h *= change
+        if abs(change - 1) < 3e-14:
+            return h
+    raise ArithmeticError('Incomplete beta did not converge')
+
+
+def ibeta(a, b, x):
+    if x <= 0:
+        return 0.0
+    if x >= 1:
+        return 1.0
+    bt = math.exp(math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
+                  + a * math.log(x) + b * math.log1p(-x))
+    if x < (a + 1) / (a + b + 2):
+        return bt * beta_cf(a, b, x) / a
+    return 1 - bt * beta_cf(b, a, 1 - x) / b
+
+
+def t_pvalue(t, df):
+    return ibeta(df / 2, 0.5, df / (df + t * t))
+
+
+def t_critical(df):
+    lo, hi = 0.0, 100.0
+    for _ in range(80):
+        mid = (lo + hi) / 2
+        if t_pvalue(mid, df) > 0.05:
+            lo = mid
+        else:
+            hi = mid
+    return (lo + hi) / 2
+
+
+def ols(y, x, names):
+    y, x = np.asarray(y, dtype=float), np.asarray(x, dtype=float)
+    if x.ndim == 1:
+        x = x[:, None]
+    X = np.column_stack([np.ones(len(y)), x])
+    beta = np.linalg.lstsq(X, y, rcond=None)[0]
+    residual = y - X @ beta
+    n, k = X.shape
+    df = n - k
+    inv = np.linalg.inv(X.T @ X)
+    sse = residual @ residual
+    s2 = sse / df
+    se = np.sqrt(np.diag(inv) * s2)
+    leverage = np.sum((X @ inv) * X, axis=1)
+    scaled = residual**2 / (1 - leverage)**2
+    hc3_cov = inv @ (X.T @ (scaled[:, None] * X)) @ inv
+    hc3_se = np.sqrt(np.diag(hc3_cov))
+    critical = t_critical(df)
+    r2 = 1 - sse / np.sum((y - y.mean())**2)
+    result = {'n': n, 'df': df, 'r2': float(r2),
+              'adj_r2': float(1 - (1 - r2) * (n - 1) / df),
+              'residual_se': float(np.sqrt(s2)), 'terms': {}}
+    for i, name in enumerate(['intercept'] + names):
+        result['terms'][name] = {
+            'coef': float(beta[i]), 'se': float(se[i]),
+            'p': t_pvalue(beta[i] / se[i], df),
+            'ci95': [float(beta[i] - critical * se[i]), float(beta[i] + critical * se[i])],
+            'hc3_se': float(hc3_se[i]),
+            'hc3_p_t_approx': t_pvalue(beta[i] / hc3_se[i], df),
+            'hc3_ci95_t_approx': [float(beta[i] - critical * hc3_se[i]), float(beta[i] + critical * hc3_se[i])]
+        }
+    return result
+
+
+assert abs(t_pvalue(1, 1) - .5) < 1e-10
+assert abs(t_critical(22) - 2.0738730679) < 1e-8
+w = openpyxl.load_workbook(SOURCE, data_only=True)
+rows = list(w.active.values)
+d = pd.DataFrame(rows[2:], columns=rows[1])
+d['source_row'] = np.arange(3, 3 + len(d))
+v = d[d['曝光'] > 0].copy()
+v['ctr_pp'] = v['封面点击率'] * 100
+v['video'] = (v['体裁'] == '视频').astype(int)
+v['published'] = pd.to_datetime(v['首次发布时间'], format='%Y年%m月%d日%H时%M分%S秒')
+v['publish_day'] = (v['published'] - v['published'].min()).dt.total_seconds() / 86400
+v['ln_imp'] = np.log(v['曝光'])
+top = v['观看量'].idxmax()
+models = {}
+for subset_name, sub in [('all', v), ('without_top', v.drop(index=top))]:
+    for outcome, label in [('曝光', 'impressions'), ('观看量', 'views')]:
+        for transform in ['raw', 'ln']:
+            y = sub[outcome] if transform == 'raw' else np.log(sub[outcome])
+            models[f'{subset_name}_{label}_{transform}'] = ols(y, sub['ctr_pp'], ['ctr_pp'])
+        models[f'{subset_name}_{label}_ln_controls'] = ols(
+            np.log(sub[outcome]), sub[['ctr_pp', 'video', 'publish_day']],
+            ['ctr_pp', 'video', 'publish_day'])
+    models[f'{subset_name}_views_ln_imp_controls'] = ols(
+        np.log(sub['观看量']), sub[['ctr_pp', 'ln_imp', 'video']],
+        ['ctr_pp', 'ln_imp', 'video'])
+
+result = {'source': SOURCE, 'total_n': len(d), 'excluded_rows': d[d['曝光'] == 0]['source_row'].tolist(),
+          'top_title': v.loc[top, '笔记标题'], 'top_source_row': int(v.loc[top, 'source_row']),
+          'models': models}
+(OUT / 'regression_results.json').write_text(json.dumps(result, ensure_ascii=False, indent=2))
+for key, model in models.items():
+    t = model['terms']['ctr_pp']
+    print(key, json.dumps({'n': model['n'], 'intercept': model['terms']['intercept']['coef'],
+          'r2': model['r2'], **t}, ensure_ascii=False))
+print('INDEPENDENT_MODEL', json.dumps(models['all_views_ln_imp_controls'], ensure_ascii=False))
+v.to_json(OUT / 'analysis_data.json', orient='records', force_ascii=False, date_format='iso')

content_statistics.py

@@ -0,0 +1,181 @@
+import itertools
+import json
+import math
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+from analyze import ols
+
+ROOT = Path(__file__).parent
+SEED = 20261017
+rng = np.random.default_rng(SEED)
+d = pd.read_json(ROOT / 'analysis_data.json')
+categories = {
+    '技术作品/工具效果展示': [7, 9, 20, 24, 25, 27],
+    'Vibe Coding指南': [11, 15, 16],
+    '模型动态/观点': [5, 6, 8, 17, 18],
+    '安全/隐私事件': [14, 19, 21],
+    '其他': [10, 12, 13, 22, 23, 26, 28],
+}
+mapping = {row: group for group, rows in categories.items() for row in rows}
+d['category'] = d.source_row.map(mapping)
+d['october'] = d['首次发布时间'].str.contains('10月').astype(int)
+assert d.category.notna().all() and len(d) == 24
+
+
+def holm(ps):
+    ps = np.asarray(ps)
+    order = np.argsort(ps)
+    corrected = np.empty_like(ps)
+    running = 0.0
+    for i, index in enumerate(order):
+        running = max(running, min(1.0, (len(ps) - i) * ps[index]))
+        corrected[index] = running
+    return corrected.tolist()
+
+
+def rank_test(a, b):
+    a, b = np.asarray(a), np.asarray(b)
+    na, nb = len(a), len(b)
+    rank2 = np.rint(pd.Series(np.r_[a, b]).rank().to_numpy() * 2).astype(int)
+    observed = int(rank2[:na].sum())
+    center = na * (na + nb + 1)
+    distance = abs(observed - center)
+    hits, total = 0, 0
+    # Every allocation of the observed ranks is equally likely under exchangeability.
+    for selection in itertools.combinations(rank2.tolist(), na):
+        hits += abs(sum(selection) - center) >= distance
+        total += 1
+    u = observed / 2 - na * (na + 1) / 2
+    superiority = u / (na * nb)
+    return {'n_a': na, 'n_b': nb, 'U': u, 'p_exact': hits / total,
+            'allocations': total, 'superiority_probability': superiority,
+            'cliff_delta': 2 * superiority - 1,
+            'median_a': float(np.median(a)), 'median_b': float(np.median(b))}
+
+
+def boot_effect(a, b, draws=50000):
+    a, b = np.asarray(a), np.asarray(b)
+    aa = a[rng.integers(0, len(a), size=(draws, len(a)))]
+    bb = b[rng.integers(0, len(b), size=(draws, len(b)))]
+    med_ratio = np.median(aa, axis=1) / np.median(bb, axis=1)
+    delta = (aa[:, :, None] > bb[:, None, :]).mean(axis=(1, 2)) - (aa[:, :, None] < bb[:, None, :]).mean(axis=(1, 2))
+    return {'median_ratio': float(np.median(a) / np.median(b)),
+            'median_ratio_ci95_percentile': np.quantile(med_ratio, [.025, .975]).tolist(),
+            'cliff_delta_ci95_percentile': np.quantile(delta, [.025, .975]).tolist(),
+            'bootstrap_draws': draws}
+
+
+def kw_permutation(sub, outcome, draws=200000):
+    groups = [g for _, g in sub.groupby('category', sort=False)]
+    sizes = np.array([len(g) for g in groups])
+    vals = np.concatenate([g[outcome].to_numpy() for g in groups])
+    ranks = pd.Series(vals).rank().to_numpy()
+    n = len(ranks)
+    _, counts = np.unique(vals, return_counts=True)
+    correction = 1 - np.sum(counts**3 - counts) / (n**3 - n)
+    starts = np.r_[0, np.cumsum(sizes)[:-1]]
+
+    def statistic(r):
+        sums = np.add.reduceat(r, starts, axis=-1)
+        return (12 / (n * (n + 1)) * np.sum(sums**2 / sizes, axis=-1) - 3 * (n + 1)) / correction
+
+    h = float(statistic(ranks))
+    hits = 0
+    for begin in range(0, draws, 5000):
+        batch = min(5000, draws - begin)
+        shuffled = ranks[np.argsort(rng.random((batch, n)), axis=1)]
+        hits += int(np.sum(statistic(shuffled) >= h - 1e-10))
+    p = (hits + 1) / (draws + 1)
+    return {'n': n, 'H': h, 'p_permutation': p, 'draws': draws,
+            'monte_carlo_se': math.sqrt(p * (1 - p) / (draws + 1)),
+            'epsilon_squared': (h - len(groups) + 1) / (n - len(groups))}
+
+
+def loo_r2(y, x):
+    y, x = np.asarray(y), np.asarray(x)
+    X = np.column_stack([np.ones(len(y)), x])
+    beta = np.linalg.lstsq(X, y, rcond=None)[0]
+    e = y - X @ beta
+    h = np.sum((X @ np.linalg.inv(X.T @ X)) * X, axis=1)
+    press = np.sum((e / (1 - h))**2)
+    return float(1 - press / np.sum((y - y.mean())**2))
+
+
+result = {'seed': SEED, 'n': len(d), 'genre': {}, 'topic_omnibus': {},
+          'topic_pairwise_views': [], 'regressions': {}, 'topic_summary': []}
+genre_tests = []
+for outcome in ['观看量', '曝光', '封面点击率']:
+    test = rank_test(d[d.video == 1][outcome], d[d.video == 0][outcome])
+    result['genre'][outcome] = test
+    genre_tests.append(test['p_exact'])
+for outcome, adjusted in zip(result['genre'], holm(genre_tests)):
+    result['genre'][outcome]['p_holm_3'] = adjusted
+result['genre']['观看量'].update(boot_effect(d[d.video == 1]['观看量'], d[d.video == 0]['观看量']))
+
+for subset_name, sub in [('all', d), ('without_top', d[d.source_row != 9]), ('september', d[d.october == 0])]:
+    if subset_name != 'all':
+        result['genre'][subset_name + '_views'] = rank_test(sub[sub.video == 1]['观看量'], sub[sub.video == 0]['观看量'])
+    for outcome in ['观看量', '曝光']:
+        covariates = ['video', 'october'] if subset_name != 'september' else ['video']
+        fit = ols(np.log(sub[outcome]), sub[covariates], covariates)
+        fit['loo_r2_log_target'] = loo_r2(np.log(sub[outcome]), sub[covariates])
+        result['regressions'][f'{subset_name}_{outcome}_format_cohort'] = fit
+
+for outcome in ['观看量', '曝光', '封面点击率']:
+    result['topic_omnibus'][outcome] = kw_permutation(d, outcome)
+for outcome, adjusted in zip(result['topic_omnibus'], holm([x['p_permutation'] for x in result['topic_omnibus'].values()])):
+    result['topic_omnibus'][outcome]['p_holm_3'] = adjusted
+result['topic_omnibus']['without_top_views'] = kw_permutation(d[d.source_row != 9], '观看量')
+
+for (ka, ra), (kb, rb) in itertools.combinations(categories.items(), 2):
+    aa, bb = d[d.source_row.isin(ra)]['观看量'], d[d.source_row.isin(rb)]['观看量']
+    test = rank_test(aa, bb)
+    test.update({'a': ka, 'b': kb})
+    result['topic_pairwise_views'].append(test)
+for test, adjusted in zip(result['topic_pairwise_views'], holm([x['p_exact'] for x in result['topic_pairwise_views']])):
+    test['p_holm_10'] = adjusted
+
+trim_pairs = []
+for (ka, ra), (kb, rb) in itertools.combinations(categories.items(), 2):
+    aa = d[d.source_row.isin(ra) & (d.source_row != 9)]['观看量']
+    bb = d[d.source_row.isin(rb) & (d.source_row != 9)]['观看量']
+    test = rank_test(aa, bb)
+    test.update({'a': ka, 'b': kb})
+    trim_pairs.append(test)
+for test, adjusted in zip(trim_pairs, holm([x['p_exact'] for x in trim_pairs])):
+    test['p_holm_10'] = adjusted
+result['topic_pairwise_views_without_top'] = trim_pairs
+result['within_video_demo_vs_other'] = rank_test(
+    d[(d.video == 1) & (d.category == '技术作品/工具效果展示')]['观看量'],
+    d[(d.video == 1) & (d.category != '技术作品/工具效果展示')]['观看量'])
+
+for name, rows in categories.items():
+    sub = d[d.source_row.isin(rows)]
+    result['topic_summary'].append({'group': name, 'n': len(sub), 'video_n': int(sub.video.sum()),
+        'median_views': float(sub['观看量'].median()), 'median_impressions': float(sub['曝光'].median()),
+        'median_ctr': float(sub['封面点击率'].median())})
+
+# Conditional exact permutation within month cohorts, using the adjusted regression coefficient.
+y = np.log(d['观看量'].to_numpy())
+video = d.video.to_numpy()
+month = d.october.to_numpy()
+indices = [np.flatnonzero(month == value) for value in [0, 1]]
+counts = [int(video[ix].sum()) for ix in indices]
+centered_y = y.copy()
+for ix in indices:
+    centered_y[ix] -= y[ix].mean()
+observed = float(centered_y @ video)
+hits, total = 0, 0
+for a in itertools.combinations(indices[0].tolist(), counts[0]):
+    sa = centered_y[list(a)].sum()
+    for b in itertools.combinations(indices[1].tolist(), counts[1]):
+        value = sa + centered_y[list(b)].sum()
+        hits += abs(value) >= abs(observed) - 1e-12
+        total += 1
+result['genre']['views_month_stratified_permutation'] = {'p_exact': hits / total, 'allocations': total}
+
+(ROOT / 'content_statistics_results.json').write_text(json.dumps(result, ensure_ascii=False, indent=2))
+print('CONTENT_STATISTICS', json.dumps(result, ensure_ascii=False, indent=2))

correlations.py

@@ -0,0 +1,69 @@
+import json
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+ROOT = Path(__file__).parent
+d = pd.read_json(ROOT / 'analysis_data.json')
+base = ['封面点击率', '点赞', '评论', '收藏', '涨粉', '分享', '人均观看时长']
+rates = []
+for metric in ['点赞', '评论', '收藏', '涨粉', '分享']:
+    name = metric + '/千次观看'
+    d[name] = d[metric] / d['观看量'] * 1000
+    rates.append(name)
+features = base + rates
+rng = np.random.default_rng(20261017)
+results = []
+B = 200000
+for outcome in ['观看量', '曝光']:
+    X = d[features].rank().to_numpy()
+    X -= X.mean(axis=0)
+    X /= np.sqrt(np.sum(X**2, axis=0))
+    y = d[outcome].rank().to_numpy()
+    y -= y.mean()
+    y /= np.sqrt(np.sum(y**2))
+    observed = y @ X
+    hits = np.zeros(len(features), dtype=int)
+    for start in range(0, B, 5000):
+        batch = min(5000, B-start)
+        shuffled = y[np.argsort(rng.random((batch, len(y))), axis=1)]
+        null = shuffled @ X
+        hits += (np.abs(null) >= np.abs(observed) - 1e-12).sum(axis=0)
+    trim = d[d.source_row != 9]
+    for j, metric in enumerate(features):
+        results.append({'outcome': outcome, 'feature': metric,
+                        'spearman': float(observed[j]),
+                        'without_top': float(trim[metric].rank().corr(trim[outcome].rank())),
+                        'p_permutation': float((hits[j]+1)/(B+1))})
+# One family of 24 exploratory tests, BH-FDR correction.
+ps = np.array([row['p_permutation'] for row in results])
+order = np.argsort(ps)
+q = np.empty_like(ps)
+running = 1.0
+for i in range(len(ps)-1, -1, -1):
+    idx = order[i]
+    running = min(running, ps[idx]*len(ps)/(i+1))
+    q[idx] = running
+for i, row in enumerate(results):
+    row['q_bh24'] = float(q[i])
+
+# Bootstrap pair resampling. Re-rank within every resample, retaining ties.
+for row in results:
+    if row['outcome'] != '观看量' or row['feature'] not in ['封面点击率', '评论/千次观看', '分享/千次观看']:
+        continue
+    n = len(d)
+    samples = rng.integers(0, n, (20000, n))
+    a = d[row['feature']].to_numpy()[samples]
+    b = d[row['outcome']].to_numpy()[samples]
+    ar = pd.DataFrame(a).rank(axis=1).to_numpy()
+    br = pd.DataFrame(b).rank(axis=1).to_numpy()
+    ar -= ar.mean(axis=1, keepdims=True)
+    br -= br.mean(axis=1, keepdims=True)
+    denom = np.sqrt((ar**2).sum(axis=1)*(br**2).sum(axis=1))
+    valid = denom > 0
+    boot = (ar[valid]*br[valid]).sum(axis=1)/denom[valid]
+    row['bootstrap_ci95'] = np.quantile(boot,[.025,.975]).tolist()
+
+(ROOT / 'correlation_results.json').write_text(json.dumps(results, ensure_ascii=False, indent=2))
+print(json.dumps(results, ensure_ascii=False, indent=2))

plot.py

@@ -0,0 +1,79 @@
+import json
+import math
+from pathlib import Path
+
+from PIL import Image, ImageDraw, ImageFont
+
+ROOT = Path(__file__).parent
+data = json.loads((ROOT / 'analysis_data.json').read_text())
+results = json.loads((ROOT / 'regression_results.json').read_text())
+models = results['models']
+image = Image.new('RGB', (1600, 1010), '#F8FAFC')
+draw = ImageDraw.Draw(image)
+FONT = '/System/Library/Fonts/Supplemental/Arial Unicode.ttf'
+
+
+def text(x, y, s, size=25, color='#334155', anchor=None):
+    draw.text((x, y), s, font=ImageFont.truetype(FONT, size), fill=color, anchor=anchor)
+
+
+text(65, 40, '点击率与曝光、观看的回归关系', 42, '#0F172A')
+text(65, 107, '24篇有效笔记；数量轴取对数，每个点代表一篇笔记', 25, '#64748B')
+colors = {'图文': '#2563EB', '视频': '#D97706'}
+for i, (genre, color) in enumerate(colors.items()):
+    xx = 1080 + i * 200
+    draw.ellipse((xx, 110, xx + 17, 127), fill=color)
+    text(xx + 30, 101, genre, 25)
+
+for panel, (outcome, label, y_min, y_max) in enumerate([
+    ('曝光', 'impressions', 1, 5.25), ('观看量', 'views', 0, 5.25)
+]):
+    left, right = 135 + panel * 780, 735 + panel * 780
+    top, bottom = 245, 725
+
+    def xy(ctr, log_quantity):
+        return (left + ctr / 21 * (right - left),
+                bottom - (log_quantity - y_min) / (y_max - y_min) * (bottom - top))
+
+    text(left - 50, 182, f'{outcome}与点击率', 30, '#0F172A')
+    for power in range(int(y_min), 6):
+        if power > y_max:
+            continue
+        yy = xy(0, power)[1]
+        draw.line((left, yy, right, yy), fill='#E2E8F0', width=2)
+        text(left - 15, yy, f'{10**power:,}', 21, '#64748B', anchor='rm')
+    for tick in [0, 5, 10, 15, 20]:
+        xx = xy(tick, y_min)[0]
+        draw.line((xx, bottom, xx, bottom + 7), fill='#94A3B8', width=2)
+        text(xx, bottom + 15, f'{tick}%', 22, '#64748B', anchor='mt')
+    draw.line((left, top, left, bottom, right, bottom), fill='#94A3B8', width=2)
+    text((left + right) / 2, bottom + 62, '封面点击率', 24, anchor='mt')
+
+    for subset, color in [('all', '#0F172A'), ('without_top', '#94A3B8')]:
+        m = models[f'{subset}_{label}_ln']
+        b0, b1 = m['terms']['intercept']['coef'], m['terms']['ctr_pp']['coef']
+        points = [xy(c / 10, (b0 + b1 * c / 10) / math.log(10)) for c in range(20, 195)]
+        for j in range(len(points) - 1):
+            if subset == 'all' or (j // 7) % 2 == 0:
+                draw.line((*points[j], *points[j + 1]), fill=color, width=4)
+
+    for row in data:
+        xx, yy = xy(row['ctr_pp'], math.log10(row[outcome]))
+        draw.ellipse((xx - 7, yy - 7, xx + 7, yy + 7), fill=colors[row['体裁']], outline='white', width=2)
+        if row['source_row'] == results['top_source_row']:
+            draw.ellipse((xx - 12, yy - 12, xx + 12, yy + 12), outline='#D97706', width=2)
+            text(xx - 15, yy - 22, '最高流量笔记', 20, '#92400E', anchor='rb')
+
+    model = models[f'all_{label}_ln']
+    beta = model['terms']['ctr_pp']['coef']
+    p = model['terms']['ctr_pp']['hc3_p_t_approx']
+    text(left - 50, 840, f'ln({outcome}) = {model["terms"]["intercept"]["coef"]:.3f} + {beta:.3f} × 点击率(%)', 23)
+    text(left - 50, 883, f'R² = {model["r2"]:.3f}     HC3稳健 p = {p:.3f}', 23)
+
+draw.line((70, 962, 120, 962), fill='#0F172A', width=4)
+text(134, 947, '全部24篇拟合', 22)
+for x in range(410, 470, 18):
+    draw.line((x, 962, x + 10, 962), fill='#94A3B8', width=4)
+text(490, 947, '去掉最高流量笔记后拟合（其余23篇）', 22)
+text(1535, 947, '数据来源：用户提供的笔记明细表', 20, '#64748B', anchor='ra')
+image.save(ROOT / '点击率回归分析.png')

render_content_statistics.py

@@ -0,0 +1,137 @@
+import json
+import math
+from pathlib import Path
+
+from PIL import Image, ImageDraw, ImageFont
+
+ROOT = Path(__file__).parent
+r = json.loads((ROOT / 'content_statistics_results.json').read_text())
+data = json.loads((ROOT / 'analysis_data.json').read_text())
+FONT = '/System/Library/Fonts/Supplemental/Arial Unicode.ttf'
+im = Image.new('RGB', (1660, 1000), '#F8FAFC')
+dr = ImageDraw.Draw(im)
+
+
+def txt(x, y, value, size=24, color='#334155', anchor=None):
+    dr.text((x, y), value, font=ImageFont.truetype(FONT, size), fill=color, anchor=anchor)
+
+
+txt(65, 38, '视频与图文的观看差异：统计检验与稳健性', 40, '#0F172A')
+txt(65, 105, '24篇有效笔记，7篇视频、17篇图文；观察性比较，不代表形式的因果效应', 24, '#64748B')
+txt(95, 185, '观看量分布', 29, '#0F172A')
+left, right, top, bottom = 150, 690, 265, 755
+
+
+def yview(v):
+    return bottom - math.log10(v) / 5.2 * (bottom - top)
+
+
+for exponent in range(6):
+    yy = yview(10**exponent)
+    dr.line((left, yy, right, yy), fill='#E2E8F0', width=2)
+    txt(left - 15, yy, f'{10**exponent:,}', 21, '#64748B', 'rm')
+dr.line((left, top, left, bottom, right, bottom), fill='#94A3B8', width=2)
+for name, flag, x, color in [('图文', 0, 300, '#2563EB'), ('视频', 1, 570, '#D97706')]:
+    values = sorted([p for p in data if p['video'] == flag], key=lambda p: p['观看量'])
+    for i, p in enumerate(values):
+        xx = x + (((i * 7) % 11) - 5) * 7
+        yy = yview(p['观看量'])
+        dr.ellipse((xx - 7, yy - 7, xx + 7, yy + 7), fill=color, outline='white', width=2)
+        if p['source_row'] == 9:
+            txt(xx, yy - 18, '77,822', 21, '#92400E', 'mb')
+    median = r['genre']['观看量']['median_a' if flag else 'median_b']
+    yy = yview(median)
+    dr.line((x - 55, yy, x + 55, yy), fill=color, width=5)
+    txt(x, bottom + 23, f'{name} n={len(values)}', 24, '#334155', 'mt')
+    txt(x, bottom + 63, f'中位数 {median:,.0f}', 22, color, 'mt')
+txt(95, 860, '双侧精确秩和检验；3项指标Holm校正', 22)
+txt(95, 900, '观看量 p=0.00113；Cliff’s δ=0.866', 24, '#0F172A')
+
+txt(820, 185, '控制发布月份后的观看量比', 29, '#0F172A')
+txt(820, 230, '点为几何均值比，横线为HC3稳健95%置信区间', 21, '#64748B')
+fleft, fright = 1035, 1570
+
+
+def xratio(v):
+    return fleft + math.log2(v) / 8 * (fright - fleft)
+
+
+for tick in [1, 4, 16, 64, 256]:
+    xx = xratio(tick)
+    dr.line((xx, 300, xx, 720), fill='#CBD5E1' if tick == 1 else '#E2E8F0', width=2)
+    txt(xx, 742, f'{tick}倍', 21, '#64748B', 'mt')
+for i, (key, label) in enumerate([
+    ('all_观看量_format_cohort', '全部24篇'),
+    ('without_top_观看量_format_cohort', '去掉最高流量'),
+    ('september_观看量_format_cohort', '仅9月18篇')
+]):
+    yy = 355 + i * 145
+    term = r['regressions'][key]['terms']['video']
+    ratio = math.exp(term['coef'])
+    lo, hi = [math.exp(v) for v in term['hc3_ci95_t_approx']]
+    txt(820, yy - 14, label, 22)
+    dr.line((xratio(lo), yy, xratio(hi), yy), fill='#0F766E', width=4)
+    for end in [lo, hi]:
+        dr.line((xratio(end), yy - 8, xratio(end), yy + 8), fill='#0F766E', width=3)
+    xx = xratio(ratio)
+    dr.ellipse((xx - 8, yy - 8, xx + 8, yy + 8), fill='#0F766E')
+    txt(xx, yy + 23, f'{ratio:.1f} [{lo:.1f}, {hi:.1f}]', 22, '#0F766E', 'mt')
+txt(820, 860, '题材与形式高度重叠，不能稳定分离各自作用。', 22)
+txt(820, 900, '视频内仅1篇非展示内容，题材比较 p=0.286。', 22, '#0F172A')
+txt(65, 962, '数据来源：笔记列表明细表.xlsx；观看量轴、倍数轴均为对数刻度。', 20, '#64748B')
+im.save(ROOT / '内容统计检验.png')
+
+lines = [
+    '# 内容与流量的统计分析', '',
+    '源文件：/Users/zzyang/Downloads/笔记列表明细表.xlsx，Sheet1。原始26篇，排除第3、4行曝光为0的最新笔记，分析24篇。', '',
+    '主要结果变量为观看量，曝光与封面点击率为辅助结果。所有结果均为单账号观察性、探索性分析。题材标签为根据标题进行的事后人工归类，未读取正文或封面。', '',
+    '## 方法', '',
+    '- 形式比较：双侧精确置换秩和检验，以秩和偏离零假设期望的绝对值为统计量，枚举346,104种分组；观看、曝光、点击率3项p值作Holm校正。检验反映分布/秩差异，并非在任意分布下只检验中位数。',
+    '- 效应量：中位数之比、Cliff’s delta、随机一篇视频观看高于随机一篇图文的样本配对比例。中位数比及delta的区间为组内重抽样50,000次的百分位Bootstrap区间。',
+    '- 题材总体比较：含并列秩修正的Kruskal–Wallis H统计量，200,000次标签置换估计p值，p=(极端次数+1)/(置换次数+1)。3项结果作Holm校正；观看量的10组两两比较另作Holm校正。',
+    '- 回归：ln(观看量)=截距+视频指示变量+10月发布指示变量。9月为参照。HC3异方差稳健标准误，使用残差自由度t分布近似计算p值及95%区间。系数指数化报告条件几何均值比。回归p值未额外作多重校正。',
+    '- 稳健性：去掉第9行最高流量笔记、仅分析9月18篇、在9月与10月内分别固定视频篇数作分层精确置换（61,200种分组）。',
+    '- 预测评估：以帽子矩阵计算留一交叉验证PRESS及log尺度R²，区别于样本内R²。',
+    '- 随机种子20261017。检验将笔记近似视为独立、可交换的观察；同账号连续发布可能不满足严格独立性。', '',
+    '## 图文和视频', '',
+    '|结果|视频中位数|图文中位数|U|原始p|Holm校正p|Cliff delta|',
+    '|---|---:|---:|---:|---:|---:|---:|'
+]
+for outcome in ['观看量', '曝光', '封面点击率']:
+    s = r['genre'][outcome]
+    a = f'{s["median_a"]:.1%}' if outcome == '封面点击率' else f'{s["median_a"]:,.0f}'
+    b = f'{s["median_b"]:.1%}' if outcome == '封面点击率' else f'{s["median_b"]:,.0f}'
+    lines.append(f'|{outcome}|{a}|{b}|{s["U"]:.1f}|{s["p_exact"]:.6f}|{s["p_holm_3"]:.6f}|{s["cliff_delta"]:.3f}|')
+s = r['genre']['观看量']
+lines += ['', f'观看量中位数比={s["median_ratio"]:.2f}，95% Bootstrap区间={s["median_ratio_ci95_percentile"]}。Cliff delta区间={s["cliff_delta_ci95_percentile"]}。视频高于图文的样本配对比例={s["superiority_probability"]:.2%}，该比例不是未来成功概率。', '',
+          '## 内容题材', '', '|题材|篇数|视频篇数|观看量中位数|曝光中位数|点击率中位数|', '|---|---:|---:|---:|---:|---:|']
+for s in r['topic_summary']:
+    lines.append(f'|{s["group"]}|{s["n"]}|{s["video_n"]}|{s["median_views"]:,.0f}|{s["median_impressions"]:,.1f}|{s["median_ctr"]:.1%}|')
+lines += ['', '|结果|H|置换p|Holm校正p|epsilon²|', '|---|---:|---:|---:|---:|']
+for outcome in ['观看量', '曝光', '封面点击率']:
+    s = r['topic_omnibus'][outcome]
+    lines.append(f'|{outcome}|{s["H"]:.3f}|{s["p_permutation"]:.6f}|{s["p_holm_3"]:.6f}|{s["epsilon_squared"]:.3f}|')
+lines += ['', 'epsilon²为秩统计量效应量，不能解读为解释了原始观看量的相同比例方差。', '',
+          '### 观看量的全部10组两两比较', '', '|A|B|原始p|Holm校正p|Cliff delta|去掉最高流量后校正p|', '|---|---|---:|---:|---:|---:|']
+for s, t in zip(r['topic_pairwise_views'], r['topic_pairwise_views_without_top']):
+    lines.append(f'|{s["a"]}|{s["b"]}|{s["p_exact"]:.6f}|{s["p_holm_10"]:.6f}|{s["cliff_delta"]:.3f}|{t["p_holm_10"]:.6f}|')
+lines += ['', '## 回归及稳健性', '', '|样本|结果|视频系数|几何均值比|HC3 95%区间|稳健p|样本内R²|留一R²(log尺度)|', '|---|---|---:|---:|---|---:|---:|---:|']
+for key, m in r['regressions'].items():
+    s = m['terms']['video']
+    lo, hi = [math.exp(a) for a in s['hc3_ci95_t_approx']]
+    lines.append(f'|{key}|log结果|{s["coef"]:.4f}|{math.exp(s["coef"]):.2f}|[{lo:.2f}, {hi:.2f}]|{s["hc3_p_t_approx"]:.6f}|{m["r2"]:.3f}|{m["loo_r2_log_target"]:.3f}|')
+lines += ['', '全样本观看量方程：ln(观看量)=4.2505+3.0640×视频−0.1126×10月。月份控制只粗略区分两个发布批次，不能替代统一发布后观察窗口。', '',
+          f'按月份分层的精确置换检验，观看量p={r["genre"]["views_month_stratified_permutation"]["p_exact"]:.6f}。', '',
+          '## 能下的结论与限制', '',
+          '1. 视频与图文的观看分布存在较大差异，去掉最高流量笔记及仅保留9月笔记后仍存在。CTR差异未检出，不能据此证明两种形式CTR等同。',
+          '2. 题材总体存在差异。展示类对模型动态类的校正p从0.039变为去掉最高流量后的0.071，具体题材结论比形式比较更脆弱。',
+          '3. 展示类6篇全部为视频，7篇视频中仅1篇非展示类。视频内的题材秩和检验p=0.286，只有1篇参照，统计功效不足。不能稳定分离题材和形式各自的作用。',
+          '4. 24篇不足以支持包含大量内容标签及交互项的可靠选题模型。观看量发布前信息模型的留一R²为0.230，预测精度有限，且留一法不能充分检验未来时间上的泛化。',
+          '5. 所有题材分组及后续敏感性分析为探索性分析，未覆盖所有可能尝试的分组方式。Holm校正仅对应明确列出的比较族，不会消除事后选标签带来的偏差。',
+          '6. 未知统一统计窗口、投流、账号粉丝基数、受众变化及平台分发状态。回归倍数为条件关联，不是把图文换成视频的因果增幅。', '',
+          '## 题材标签与源行号', '']
+for summary in r['topic_summary']:
+    matching = [x for x in data if (x['source_row'] in {'技术作品/工具效果展示':[7,9,20,24,25,27], 'Vibe Coding指南':[11,15,16], '模型动态/观点':[5,6,8,17,18], '安全/隐私事件':[14,19,21], '其他':[10,12,13,22,23,26,28]}[summary['group']])]
+    lines.append(f'### {summary["group"]}')
+    lines += ['', *[f'- 第{x["source_row"]}行：{x["笔记标题"]}' for x in matching], '']
+(ROOT / '内容统计分析报告.md').write_text('\n'.join(lines))