137 lines
5.0 KiB
Python
137 lines
5.0 KiB
Python
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')
|