commit 68383b90a9ab066f70163a743d2dc2a7b8f91cfd Author: Codex Bot Date: Thu Nov 5 01:04:53 2026 +0800 Import project sources diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a654501 --- /dev/null +++ b/.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 diff --git a/analyze.py b/analyze.py new file mode 100644 index 0000000..9650d2d --- /dev/null +++ b/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') diff --git a/content_statistics.py b/content_statistics.py new file mode 100644 index 0000000..a99c900 --- /dev/null +++ b/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)) diff --git a/correlations.py b/correlations.py new file mode 100644 index 0000000..5022af8 --- /dev/null +++ b/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)) diff --git a/plot.py b/plot.py new file mode 100644 index 0000000..4e60662 --- /dev/null +++ b/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') diff --git a/render_content_statistics.py b/render_content_statistics.py new file mode 100644 index 0000000..a82a80e --- /dev/null +++ b/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))