AI Scientist Bench
Results overview / Analyze data

Run a statistical analysis

Write and run code for a specified statistical analysis.

Sonnet 5·8 examples·One attempt with each setup·2026-09-15

One more completed analysis

Baseline passed seven of eight tasks; checklist and skill passed eight. The extra pass reflects a baseline tool-call failure, not a wrong statistical result.

AI alone
7/8
analyses passed
$0.0290 · 20.3s per case
With checklist
8/8
analyses passed
$0.0166 · 13.1s per case
With skill
8/8
analyses passed
$0.0277 · 17.4s per case

Skill tested: K-Dense: statistical-analysis ↗. “AI alone” uses the task instructions and required tools; the other setups add a checklist or skill guidance. Everything else is held the same.

How this was scored, and what it does not establish

A case passes only when the submitted program passes numerical checks on both the original data and a changed version. Failed tool calls remain failed cases.

  • Eight tasks share four public datasets and one synthetic fixture; the statistical methods are prescribed.
  • Passing programs ran on original and hidden changed data, but this does not establish general analysis competence.
  • The fixed environment supplies NumPy, SciPy and pandas; the full upstream dependency stack is unavailable.
  • Multiplicity checks do not change the rejection set; a test-only outcome perturbation would strengthen the leakage check.

SKILL.md supplied; no explicit reference-file reads recorded.

Next useful test: First validate one complete reproduction in a native agent, then test independent datasets that require choosing a statistical method.

Examples and actual answers

Open an example to see all three attempts.

Penguins: unequal-variance comparisonpenguins.csvAI alonePassedWith checklistPassedWith skillPassed

Estimate male minus female Adelie body mass with Welch uncertainty and an effect size.

Source dataset ↗ · Source cluster: penguins

AI alone

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_independent",
  "n_left": 73,
  "n_right": 73,
  "difference": 674.6575342465753,
  "se": 51.39744312894034,
  "df": 135.69357893530383,
  "p_value": 6.402319748031793e-26,
  "ci_low": 573.0139063407537,
  "ci_high": 776.301162152397,
  "hedges_g": 2.1613425922546154
}
Submitted code excerpt
import pandas as pd, numpy as np, json
from scipy import stats

df = pd.read_csv("data.csv")
df = df[df['sex'].isin(['male','female'])]
df = df.dropna(subset=['body_mass_g'])
df['body_mass_g'] = pd.to_numeric(df['body_mass_g'], errors='coerce')
df = df.dropna(subset=['body_mass_g'])

m = df.loc[df['sex']=='male','body_mass_g'].values
f = df.loc[df['sex']=='female','body_mass_g'].values

n1, n2 = len(m), len(f)
mean1, mean2 = m.mean(), f.mean()
var1, var2 = m.var(ddof=1), f.var(ddof=1)

diff = mean1 - mean2
se = np.sqrt(var1/n1 + var2/n2)
df_welch = (var1/n1 + var2/n2)**2 / ((var1/n1)**2/(n1-1) + (var2/n2)**2/(n2-1))

t = diff / se
p = 2*stats.t.sf(abs(t), df_welch)

tcrit = stats.t.ppf(0.975, df_welch)
ci_low = diff - tcrit*se
ci_high = diff + tcrit*se

sp = np.sqrt(((n1-1)*var1 + (n2-1)*var2)/(n1+n2-2))
J = 1 - 3/(4*(n1+n2)-9)
g = J * diff / sp

result = {
    "method": "welch_independent",
    "n_left": n1,
    "n_right": n2,
    "difference": float(diff),
    "se": float(se),
    "df": float(df_welch),
    "p_value": float(p),
    "ci_low": float(ci_low),
    "ci_high": float(ci_high),
    "hedges_g": float(g)
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0114 · 9.6s

With checklist

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_independent",
  "n_left": 73,
  "n_right": 73,
  "difference": 674.6575342465753,
  "se": 51.39744312894034,
  "df": 135.69357893530383,
  "p_value": 6.402319748031793e-26,
  "ci_low": 573.0139063407537,
  "ci_high": 776.301162152397,
  "hedges_g": 2.1613425922546154
}
Submitted code excerpt
import pandas as pd, numpy as np, json
from scipy import stats

df = pd.read_csv("data.csv")
df = df[df['sex'].isin(['male','female'])]
df = df.dropna(subset=['body_mass_g'])
df['body_mass_g'] = pd.to_numeric(df['body_mass_g'], errors='coerce')
df = df.dropna(subset=['body_mass_g'])

m = df.loc[df['sex']=='male','body_mass_g'].astype(float)
f = df.loc[df['sex']=='female','body_mass_g'].astype(float)

n1, n2 = len(m), len(f)
mean1, mean2 = m.mean(), f.mean()
var1, var2 = m.var(ddof=1), f.var(ddof=1)

diff = mean1 - mean2
se = np.sqrt(var1/n1 + var2/n2)
df_welch = (var1/n1 + var2/n2)**2 / ((var1/n1)**2/(n1-1) + (var2/n2)**2/(n2-1))

t_stat = diff/se
p_value = 2*stats.t.sf(np.abs(t_stat), df_welch)

ci_low = diff - stats.t.ppf(0.975, df_welch)*se
ci_high = diff + stats.t.ppf(0.975, df_welch)*se

pooled_sd = np.sqrt(((n1-1)*var1 + (n2-1)*var2)/(n1+n2-2))
J = 1 - 3/(4*(n1+n2)-9)
hedges_g = J * diff / pooled_sd

result = {
    "method": "welch_independent",
    "n_left": int(n1),
    "n_right": int(n2),
    "difference": float(diff),
    "se": float(se),
    "df": float(df_welch),
    "p_value": float(p_value),
    "ci_low": float(ci_low),
    "ci_high": float(ci_high),
    "hedges_g": …

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0127 · 10.1s

With skill

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_independent",
  "n_left": 73,
  "n_right": 73,
  "difference": 674.6575342465753,
  "se": 51.39744312894034,
  "df": 135.69357893530383,
  "p_value": 6.402319748031793e-26,
  "ci_low": 573.0139063407537,
  "ci_high": 776.301162152397,
  "hedges_g": 2.1613425922546154
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv('data.csv')
df = df[df['sex'].isin(['male','female'])]
df = df.dropna(subset=['body_mass_g'])
male = df[df['sex']=='male']['body_mass_g'].astype(float)
female = df[df['sex']=='female']['body_mass_g'].astype(float)

n1,n2 = len(male), len(female)
m1,m2 = male.mean(), female.mean()
v1,v2 = male.var(ddof=1), female.var(ddof=1)

diff = m1-m2
se = np.sqrt(v1/n1+v2/n2)
df_welch = (v1/n1+v2/n2)**2 / ((v1/n1)**2/(n1-1)+(v2/n2)**2/(n2-1))
t = diff/se
p = 2*stats.t.sf(abs(t), df_welch)
ci_low = diff - stats.t.ppf(0.975, df_welch)*se
ci_high = diff + stats.t.ppf(0.975, df_welch)*se

sp = np.sqrt(((n1-1)*v1+(n2-1)*v2)/(n1+n2-2))
d = diff/sp
J = 1-3/(4*(n1+n2)-9)
g = d*J

res = {
"method":"welch_independent",
"n_left": n1,
"n_right": n2,
"difference": diff,
"se": se,
"df": df_welch,
"p_value": p,
"ci_low": ci_low,
"ci_high": ci_high,
"hedges_g": g
}
print(json.dumps(res))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0448 · 16.3s

Recorded case analyze-penguin-welch · Run four-track-canary-v3. Output excerpts are bounded; the full submitted answer was scored.

Norris: calibration regressionNorris.datAI alonePassedWith checklistPassedWith skillPassed

Fit an intercept and slope to the calibration observations and report classical uncertainty.

Source dataset ↗ · Source cluster: nist_norris

AI alone

Both original and changed-input executions passed.

Numerical answer
{
  "method": "ols_with_intercept",
  "n": 36,
  "intercept": -0.2623230737740414,
  "slope": 1.0021168180204545,
  "slope_ci_low": 1.0012433657355782,
  "slope_ci_high": 1.0029902703053308,
  "residual_sd": 0.884796396144381,
  "r_squared": 0.9999937458837117
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv("data.csv")
x = df['x'].values.astype(float)
y = df['y'].values.astype(float)
n = len(x)

xm = x.mean(); ym = y.mean()
Sxx = ((x-xm)**2).sum()
Sxy = ((x-xm)*(y-ym)).sum()
slope = Sxy/Sxx
intercept = ym - slope*xm
yhat = intercept + slope*x
resid = y - yhat
ss_res = (resid**2).sum()
dof = n-2
resid_var = ss_res/dof
resid_sd = np.sqrt(resid_var)
se_slope = np.sqrt(resid_var/Sxx)
tcrit = stats.t.ppf(0.975, dof)
ci_low = slope - tcrit*se_slope
ci_high = slope + tcrit*se_slope
ss_tot = ((y-ym)**2).sum()
r2 = 1 - ss_res/ss_tot

result = {
    "method": "ols_with_intercept",
    "n": n,
    "intercept": intercept,
    "slope": slope,
    "slope_ci_low": ci_low,
    "slope_ci_high": ci_high,
    "residual_sd": resid_sd,
    "r_squared": r2
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0187 · 29.3s

With checklist

Both original and changed-input executions passed.

Numerical answer
{
  "method": "ols_with_intercept",
  "n": 36,
  "intercept": -0.2623230737737886,
  "slope": 1.0021168180204536,
  "slope_ci_low": 1.0012433657355773,
  "slope_ci_high": 1.00299027030533,
  "residual_sd": 0.8847963961443822,
  "r_squared": 0.9999937458837117
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv("data.csv")
df = df.dropna(subset=["y","x"])
x = df["x"].values.astype(float)
y = df["y"].values.astype(float)
n = len(x)

X = np.column_stack([np.ones(n), x])
beta, res, rank, sv = np.linalg.lstsq(X, y, rcond=None)
intercept, slope = beta
yhat = X @ beta
resid = y - yhat
dof = n - 2
rss = np.sum(resid**2)
residual_sd = np.sqrt(rss/dof)
XtX_inv = np.linalg.inv(X.T @ X)
se_slope = residual_sd * np.sqrt(XtX_inv[1,1])
tcrit = stats.t.ppf(0.975, dof)
ci_low = slope - tcrit*se_slope
ci_high = slope + tcrit*se_slope

ss_tot = np.sum((y-np.mean(y))**2)
r_squared = 1 - rss/ss_tot

result = {
    "method": "ols_with_intercept",
    "n": n,
    "intercept": float(intercept),
    "slope": float(slope),
    "slope_ci_low": float(ci_low),
    "slope_ci_high": float(ci_high),
    "residual_sd": float(residual_sd),
    "r_squared": float(r_squared)
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0168 · 15.9s

With skill

Both original and changed-input executions passed.

Numerical answer
{
  "method": "ols_with_intercept",
  "n": 36,
  "intercept": -0.2623230737737886,
  "slope": 1.0021168180204536,
  "slope_ci_low": 1.0012433657355773,
  "slope_ci_high": 1.00299027030533,
  "residual_sd": 0.8847963961443822,
  "r_squared": 0.9999937458837117
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv("data.csv")
x = df['x'].values
y = df['y'].values
n = len(x)

X = np.column_stack([np.ones(n), x])
beta, res, rank, sv = np.linalg.lstsq(X, y, rcond=None)
intercept, slope = beta

yhat = X @ beta
resid = y - yhat
dof = n - 2
rss = np.sum(resid**2)
residual_sd = np.sqrt(rss/dof)

XtX_inv = np.linalg.inv(X.T @ X)
se = np.sqrt(residual_sd**2 * np.diag(XtX_inv))
se_slope = se[1]

tcrit = stats.t.ppf(0.975, dof)
ci_low = slope - tcrit*se_slope
ci_high = slope + tcrit*se_slope

ss_tot = np.sum((y-np.mean(y))**2)
r2 = 1 - rss/ss_tot

result = {
    "method": "ols_with_intercept",
    "n": n,
    "intercept": float(intercept),
    "slope": float(slope),
    "slope_ci_low": float(ci_low),
    "slope_ci_high": float(ci_high),
    "residual_sd": float(residual_sd),
    "r_squared": float(r2)
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0195 · 26.1s

Recorded case analyze-norris-calibration · Run four-track-canary-v3. Output excerpts are bounded; the full submitted answer was scored.

Mavro: ordered measurementsMavro.datAI alonePassedWith checklistPassedWith skillPassed

Compare independent-observation uncertainty with five consecutive block means.

Source dataset ↗ · Source cluster: nist_mavro

AI alone

Both original and changed-input executions passed.

Numerical answer
{
  "method": "ordered_block_mean_t",
  "n": 50,
  "mean": 2.0018559999999996,
  "sample_sd": 0.0004291234540030854,
  "lag1_rho": 0.9379891834382406,
  "iid_se": 6.068722085835504e-05,
  "block_se": 0.0001895415521726317,
  "block_ci_low": 2.0013297482852392,
  "block_ci_high": 2.002382251714762
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv('data.csv')
df = df.sort_values('sequence')
y = df['transmittance'].values
n = len(y)
mean = y.mean()
sample_sd = y.std(ddof=1)
iid_se = sample_sd/np.sqrt(n)
num = np.sum((y[1:]-mean)*(y[:-1]-mean))
den = np.sum((y-mean)**2)
lag1_rho = num/den

blocks = y[:50].reshape(5,10)
block_means = blocks.mean(axis=1)
block_sd = block_means.std(ddof=1)
block_se = block_sd/np.sqrt(5)
tcrit = stats.t.ppf(0.975,4)
bm = block_means.mean()
ci_low = bm - tcrit*block_se
ci_high = bm + tcrit*block_se

res = {
 "method":"ordered_block_mean_t",
 "n": int(n),
 "mean": float(mean),
 "sample_sd": float(sample_sd),
 "lag1_rho": float(lag1_rho),
 "iid_se": float(iid_se),
 "block_se": float(block_se),
 "block_ci_low": float(ci_low),
 "block_ci_high": float(ci_high)
}
print(json.dumps(res))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0193 · 12.5s

With checklist

Both original and changed-input executions passed.

Numerical answer
{
  "method": "ordered_block_mean_t",
  "n": 50,
  "mean": 2.0018559999999996,
  "sample_sd": 0.0004291234540030854,
  "lag1_rho": 0.9379891834382406,
  "iid_se": 6.068722085835504e-05,
  "block_se": 0.0001895415521726317,
  "block_ci_low": 2.0013297482852392,
  "block_ci_high": 2.002382251714762
}
Submitted code excerpt
import pandas as pd, numpy as np, json
from scipy import stats

df = pd.read_csv("data.csv").sort_values("sequence")
y = df["transmittance"].to_numpy(dtype=float)
n = len(y)
mean = y.mean()
sample_sd = y.std(ddof=1)
iid_se = sample_sd/np.sqrt(n)

num = np.sum((y[1:]-mean)*(y[:-1]-mean))
den = np.sum((y-mean)**2)
lag1_rho = num/den

nblocks = 5
bsize = n//nblocks
blocks = y[:nblocks*bsize].reshape(nblocks,bsize)
block_means = blocks.mean(axis=1)
bm_mean = block_means.mean()
bm_sd = block_means.std(ddof=1)
block_se = bm_sd/np.sqrt(nblocks)
tcrit = stats.t.ppf(0.975, 4)
block_ci_low = bm_mean - tcrit*block_se
block_ci_high = bm_mean + tcrit*block_se

result = {
    "method": "ordered_block_mean_t",
    "n": int(n),
    "mean": float(mean),
    "sample_sd": float(sample_sd),
    "lag1_rho": float(lag1_rho),
    "iid_se": float(iid_se),
    "block_se": float(block_se),
    "block_ci_low": float(block_ci_low),
    "block_ci_high": float(block_ci_high)
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0080 · 8.3s

With skill

Both original and changed-input executions passed.

Numerical answer
{
  "method": "ordered_block_mean_t",
  "n": 50,
  "mean": 2.0018559999999996,
  "sample_sd": 0.0004291234540030854,
  "lag1_rho": 0.9379891834382406,
  "iid_se": 6.068722085835504e-05,
  "block_se": 0.0001895415521726317,
  "block_ci_low": 2.0013297482852392,
  "block_ci_high": 2.002382251714762
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv("data.csv").sort_values("sequence")
y = df["transmittance"].values
n = len(y)
mean = y.mean()
sample_sd = y.std(ddof=1)
iid_se = sample_sd/np.sqrt(n)

num = np.sum((y[1:]-mean)*(y[:-1]-mean))
den = np.sum((y-mean)**2)
lag1_rho = num/den

blocks = y.reshape(5,10)
block_means = blocks.mean(axis=1)
block_sd = block_means.std(ddof=1)
block_se = block_sd/np.sqrt(5)
tcrit = stats.t.ppf(0.975,4)
block_ci_low = block_means.mean() - tcrit*block_se
block_ci_high = block_means.mean() + tcrit*block_se

res = {
    "method": "ordered_block_mean_t",
    "n": int(n),
    "mean": float(mean),
    "sample_sd": float(sample_sd),
    "lag1_rho": float(lag1_rho),
    "iid_se": float(iid_se),
    "block_se": float(block_se),
    "block_ci_low": float(block_ci_low),
    "block_ci_high": float(block_ci_high),
}
print(json.dumps(res, indent=2))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0198 · 14.6s

Recorded case analyze-mavro-blocks · Run four-track-canary-v3. Output excerpts are bounded; the full submitted answer was scored.

Silver: instrument comparisonAtmWtAg.datAI alonePassedWith checklistPassedWith skillPassed

Compare two independent instruments while preserving very small atomic-weight differences.

Source dataset ↗ · Source cluster: nist_silver

AI alone

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_independent",
  "n_left": 24,
  "n_right": 24,
  "difference": -1.7412500028513023e-05,
  "se": 4.360389250340865e-06,
  "df": 43.25183428311259,
  "p_value": 0.00024855676056512177,
  "ci_low": -2.620458259250799e-05,
  "ci_high": -8.620417464518057e-06,
  "hedges_g": -1.1338788698859152
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv("data.csv")
g1 = df[df.instrument==2]["atomic_weight"].values  # left
g2 = df[df.instrument==1]["atomic_weight"].values  # right

n1, n2 = len(g1), len(g2)
m1, m2 = g1.mean(), g2.mean()
v1, v2 = g1.var(ddof=1), g2.var(ddof=1)

diff = m1 - m2
se = np.sqrt(v1/n1 + v2/n2)
df_welch = (v1/n1 + v2/n2)**2 / ((v1/n1)**2/(n1-1) + (v2/n2)**2/(n2-1))

t_stat = diff/se
p_value = 2*stats.t.sf(abs(t_stat), df_welch)

ci_low = diff - stats.t.ppf(0.975, df_welch)*se
ci_high = diff + stats.t.ppf(0.975, df_welch)*se

sp = np.sqrt(((n1-1)*v1+(n2-1)*v2)/(n1+n2-2))
J = 1 - 3/(4*(n1+n2)-9)
hedges_g = J*diff/sp

result = {
    "method": "welch_independent",
    "n_left": n1,
    "n_right": n2,
    "difference": float(diff),
    "se": float(se),
    "df": float(df_welch),
    "p_value": float(p_value),
    "ci_low": float(ci_low),
    "ci_high": float(ci_high),
    "hedges_g": float(hedges_g)
}

print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0099 · 8.8s

With checklist

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_independent",
  "n_left": 24,
  "n_right": 24,
  "difference": -1.7412500028513023e-05,
  "se": 4.360389250340865e-06,
  "df": 43.25183428311259,
  "p_value": 0.00024855676056512177,
  "ci_low": -2.620458259250799e-05,
  "ci_high": -8.620417464518057e-06,
  "hedges_g": -1.133878869885915
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv("data.csv")
df = df.dropna(subset=["instrument","atomic_weight"])

g1 = df.loc[df.instrument==1, "atomic_weight"].astype(float)
g2 = df.loc[df.instrument==2, "atomic_weight"].astype(float)

n_left = len(g2)
n_right = len(g1)

mean1, mean2 = g1.mean(), g2.mean()
var1, var2 = g1.var(ddof=1), g2.var(ddof=1)

difference = mean2 - mean1

se = np.sqrt(var1/n_right + var2/n_left)

df_welch = (var1/n_right + var2/n_left)**2 / (
    (var1/n_right)**2/(n_right-1) + (var2/n_left)**2/(n_left-1)
)

t_stat = difference / se
p_value = 2*stats.t.sf(np.abs(t_stat), df_welch)

tcrit = stats.t.ppf(0.975, df_welch)
ci_low = difference - tcrit*se
ci_high = difference + tcrit*se

# pooled SD for hedges g
pooled_sd = np.sqrt(((n_right-1)*var1 + (n_left-1)*var2) / (n_right + n_left - 2))
J = 1 - 3/(4*(n_left+n_right)-9)
hedges_g = J * (difference / pooled_sd)

result = {
    "method": "welch_independent",
    "n_left": int(n_left),
    "n_right": int(n_right),
    "difference": float(difference),
    "se": float(se),
    "df": float(df_welch),
    "p_value": float(p_value),
    "ci_low": float(ci_low),
    "ci_high": f…

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0092 · 9.7s

With skill

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_independent",
  "n_left": 24,
  "n_right": 24,
  "difference": -1.7412500028513023e-05,
  "se": 4.360389250340865e-06,
  "df": 43.25183428311259,
  "p_value": 0.00024855676056512177,
  "ci_low": -2.620458259250799e-05,
  "ci_high": -8.620417464518057e-06,
  "hedges_g": -1.133878869885915
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv('data.csv')
g1 = df[df.instrument==2]['atomic_weight'].values  # left = instrument 2
g2 = df[df.instrument==1]['atomic_weight'].values  # right = instrument 1

n_left = len(g1)
n_right = len(g2)
m1, m2 = g1.mean(), g2.mean()
v1, v2 = g1.var(ddof=1), g2.var(ddof=1)

diff = m1 - m2
se = np.sqrt(v1/n_left + v2/n_right)
df_welch = (v1/n_left + v2/n_right)**2 / ((v1/n_left)**2/(n_left-1) + (v2/n_right)**2/(n_right-1))
t = diff/se
p = 2*stats.t.sf(np.abs(t), df_welch)
tcrit = stats.t.ppf(0.975, df_welch)
ci_low = diff - tcrit*se
ci_high = diff + tcrit*se

sp = np.sqrt(((n_left-1)*v1 + (n_right-1)*v2)/(n_left+n_right-2))
d = diff/sp
J = 1 - 3/(4*(n_left+n_right)-9)
g = d*J

result = {
    "method": "welch_independent",
    "n_left": n_left,
    "n_right": n_right,
    "difference": float(diff),
    "se": float(se),
    "df": float(df_welch),
    "p_value": float(p),
    "ci_low": float(ci_low),
    "ci_high": float(ci_high),
    "hedges_g": float(g)
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0219 · 17.3s

Recorded case analyze-silver-instruments · Run four-track-canary-v3. Output excerpts are bounded; the full submitted answer was scored.

Synthetic data: incomplete pairspaired_synthetic.csvAI alonePassedWith checklistPassedWith skillPassed

Compute after-minus-before changes using complete matched subjects only.

Synthetic fixture · Source cluster: synthetic_paired

AI alone

Both original and changed-input executions passed.

Numerical answer
{
  "method": "paired_complete_cases",
  "n_pairs": 12,
  "excluded_subjects": 2,
  "mean_change": 1.508333333333334,
  "sd_change": 1.1843358962522066,
  "p_value": 0.001042637825800151,
  "ci_low": 0.7558422052400875,
  "ci_high": 2.2608244614265804,
  "cohen_dz": 1.2735688735825765
}
Submitted code excerpt
import pandas as pd, numpy as np, json
from scipy import stats

df = pd.read_csv('data.csv')
total = len(df)
complete = df.dropna(subset=['before','after'])
excluded = total - len(complete)

diff = complete['after'] - complete['before']
n = len(diff)
mean_change = diff.mean()
sd_change = diff.std(ddof=1)
se = sd_change/np.sqrt(n)
tstat = mean_change/se
p_value = 2*stats.t.sf(abs(tstat), df=n-1)
tcrit = stats.t.ppf(0.975, df=n-1)
ci_low = mean_change - tcrit*se
ci_high = mean_change + tcrit*se
cohen_dz = mean_change/sd_change

result = {
    "method": "paired_complete_cases",
    "n_pairs": int(n),
    "excluded_subjects": int(excluded),
    "mean_change": float(mean_change),
    "sd_change": float(sd_change),
    "p_value": float(p_value),
    "ci_low": float(ci_low),
    "ci_high": float(ci_high),
    "cohen_dz": float(cohen_dz)
}

print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0079 · 7.1s

With checklist

Both original and changed-input executions passed.

Numerical answer
{
  "method": "paired_complete_cases",
  "n_pairs": 12,
  "excluded_subjects": 2,
  "mean_change": 1.508333333333334,
  "sd_change": 1.1843358962522066,
  "p_value": 0.001042637825800151,
  "ci_low": 0.7558422052400875,
  "ci_high": 2.2608244614265804,
  "cohen_dz": 1.2735688735825765
}
Submitted code excerpt
import pandas as pd, numpy as np, json
from scipy import stats

df = pd.read_csv("data.csv")
total_subjects = len(df)
complete = df.dropna(subset=["before","after"])
excluded_subjects = total_subjects - len(complete)

diffs = complete["after"] - complete["before"]
n = len(diffs)
mean_change = diffs.mean()
sd_change = diffs.std(ddof=1)
se = sd_change / np.sqrt(n)
tstat = mean_change / se
p_value = 2*stats.t.sf(abs(tstat), df=n-1)
tcrit = stats.t.ppf(0.975, df=n-1)
ci_low = mean_change - tcrit*se
ci_high = mean_change + tcrit*se
cohen_dz = mean_change / sd_change

result = {
    "method": "paired_complete_cases",
    "n_pairs": int(n),
    "excluded_subjects": int(excluded_subjects),
    "mean_change": float(mean_change),
    "sd_change": float(sd_change),
    "p_value": float(p_value),
    "ci_low": float(ci_low),
    "ci_high": float(ci_high),
    "cohen_dz": float(cohen_dz)
}

print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0065 · 9.1s

With skill

Both original and changed-input executions passed.

Numerical answer
{
  "method": "paired_complete_cases",
  "n_pairs": 12,
  "excluded_subjects": 2,
  "mean_change": 1.508333333333334,
  "sd_change": 1.1843358962522066,
  "p_value": 0.001042637825800151,
  "ci_low": 0.7558422052400875,
  "ci_high": 2.2608244614265804,
  "cohen_dz": 1.2735688735825765
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv('data.csv')
total = len(df)
complete = df.dropna(subset=['before','after'])
excluded = total - len(complete)
diff = complete['after'] - complete['before']
n = len(diff)
mean_change = diff.mean()
sd_change = diff.std(ddof=1)
se = sd_change/np.sqrt(n)
t_stat = mean_change/se
p_value = 2*stats.t.sf(abs(t_stat), df=n-1)
ci = stats.t.interval(0.95, n-1, loc=mean_change, scale=se)
cohen_dz = mean_change/sd_change

result = {
    "method": "paired_complete_cases",
    "n_pairs": int(n),
    "excluded_subjects": int(excluded),
    "mean_change": float(mean_change),
    "sd_change": float(sd_change),
    "p_value": float(p_value),
    "ci_low": float(ci[0]),
    "ci_high": float(ci[1]),
    "cohen_dz": float(cohen_dz)
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0168 · 16.3s

Recorded case analyze-paired-missing · Run new-sources-and-analysis-v3. Output excerpts are bounded; the full submitted answer was scored.

Penguins: species-adjusted regressionpenguins.csvAI aloneFailedWith checklistPassedWith skillPassed

Regress body mass on flipper length and species, retaining rows with missing non-model variables.

Source dataset ↗ · Source cluster: penguins

AI alone

A tool call exhausted its output allowance and omitted the code argument; no numerical answer was graded.

Recorded score: 0/1 · $0.1353 · 71.1s

With checklist

Both original and changed-input executions passed.

Numerical answer
{
  "method": "ols_species_adjusted",
  "n": 342,
  "intercept": -4031.476890693633,
  "flipper_slope": 40.705400777284765,
  "chinstrap_offset": -206.51012033971995,
  "gentoo_offset": 266.80960317921796,
  "slope_ci_low": 34.664682514400155,
  "slope_ci_high": 46.746119040169376,
  "residual_sd": 375.53507478218387
}
Submitted code excerpt
import pandas as pd
import numpy as np
from scipy import stats
import json

df = pd.read_csv('data.csv', na_values=['NA'])

needed = ['body_mass_g','flipper_length_mm','species']
d = df.dropna(subset=needed).copy()

y = d['body_mass_g'].astype(float).values
flipper = d['flipper_length_mm'].astype(float).values
species = d['species'].astype(str)

n = len(d)
chin = (species == 'Chinstrap').astype(float).values
gent = (species == 'Gentoo').astype(float).values

X = np.column_stack([np.ones(n), flipper, chin, gent])

beta, residuals, rank, sv = np.linalg.lstsq(X, y, rcond=None)

fitted = X @ beta
resid = y - fitted
p = X.shape[1]
dof = n - p
sigma2 = np.sum(resid**2) / dof
residual_sd = np.sqrt(sigma2)

XtX_inv = np.linalg.inv(X.T @ X)
cov_beta = sigma2 * XtX_inv
se = np.sqrt(np.diag(cov_beta))

tcrit = stats.t.ppf(0.975, dof)

intercept = beta[0]
flipper_slope = beta[1]
chinstrap_offset = beta[2]
gentoo_offset = beta[3]

slope_se = se[1]
slope_ci_low = flipper_slope - tcrit*slope_se
slope_ci_high = flipper_slope + tcrit*slope_se

result = {
    "method": "ols_species_adjusted",
    "n": int(n),
    "intercept": float(intercept),
    "flipper_slope": float(flipper_slope),
    "chinstr…

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0287 · 15.2s

With skill

Both original and changed-input executions passed.

Numerical answer
{
  "method": "ols_species_adjusted",
  "n": 342,
  "intercept": -4031.476890693633,
  "flipper_slope": 40.705400777284765,
  "chinstrap_offset": -206.51012033971995,
  "gentoo_offset": 266.80960317921796,
  "slope_ci_low": 34.664682514400155,
  "slope_ci_high": 46.746119040169376,
  "residual_sd": 375.5350747821838
}
Submitted code excerpt
import pandas as pd, numpy as np, json
from scipy import stats

df = pd.read_csv("data.csv")
df2 = df.dropna(subset=["body_mass_g","flipper_length_mm","species"])
df2 = df2[df2['species'].notna()]

y = df2['body_mass_g'].astype(float).values
sp = df2['species']
X = pd.DataFrame({
    'const':1.0,
    'flipper': df2['flipper_length_mm'].astype(float).values,
    'chinstrap': (sp=='Chinstrap').astype(float).values,
    'gentoo': (sp=='Gentoo').astype(float).values,
})
Xm = X.values
n,k = Xm.shape
beta, _, _, _ = np.linalg.lstsq(Xm, y, rcond=None)
resid = y - Xm@beta
dof = n-k
sigma2 = (resid@resid)/dof
XtX_inv = np.linalg.inv(Xm.T@Xm)
cov = sigma2*XtX_inv
se = np.sqrt(np.diag(cov))

tcrit = stats.t.ppf(0.975, dof)

slope = beta[1]
slope_se = se[1]
ci_low = slope - tcrit*slope_se
ci_high = slope + tcrit*slope_se

residual_sd = np.sqrt(sigma2)

result = {
    "method": "ols_species_adjusted",
    "n": int(n),
    "intercept": float(beta[0]),
    "flipper_slope": float(slope),
    "chinstrap_offset": float(beta[2]),
    "gentoo_offset": float(beta[3]),
    "slope_ci_low": float(ci_low),
    "slope_ci_high": float(ci_high),
    "residual_sd": float(residual_sd),
}
print(json.dumps(resul…

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0622 · 18.6s

Recorded case analyze-penguin-adjusted · Run new-sources-and-analysis-v3. Output excerpts are bounded; the full submitted answer was scored.

Norris: held-out predictionNorris.datAI alonePassedWith checklistPassedWith skillPassed

Fit only the training rows, then calculate prediction error on held-out rows.

Source dataset ↗ · Source cluster: nist_norris

AI alone

Both original and changed-input executions passed.

Numerical answer
{
  "method": "train_only_ols_holdout",
  "n_train": 24,
  "n_test": 12,
  "intercept": -0.25917168496622917,
  "slope": 1.0029666153104575,
  "test_rmse": 1.4837478210602169,
  "test_mae": 1.1711495352409773
}
Submitted code excerpt
import pandas as pd, numpy as np, json

df = pd.read_csv("data.csv")
train = df[df.split=="train"]
test = df[df.split=="test"]

x = train.x.values
y = train.y.values
n = len(x)
slope, intercept = np.polyfit(x, y, 1)

xt = test.x.values
yt = test.y.values
pred = intercept + slope*xt
resid = yt - pred
rmse = np.sqrt(np.mean(resid**2))
mae = np.mean(np.abs(resid))

result = {
    "method": "train_only_ols_holdout",
    "n_train": int(n),
    "n_test": int(len(xt)),
    "intercept": float(intercept),
    "slope": float(slope),
    "test_rmse": float(rmse),
    "test_mae": float(mae)
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0197 · 15.1s

With checklist

Both original and changed-input executions passed.

Numerical answer
{
  "method": "train_only_ols_holdout",
  "n_train": 24,
  "n_test": 12,
  "intercept": -0.25917168496622917,
  "slope": 1.0029666153104575,
  "test_rmse": 1.4837478210602169,
  "test_mae": 1.1711495352409773
}
Submitted code excerpt
import pandas as pd, numpy as np, json

df = pd.read_csv("data.csv")
train = df[df.split=="train"].dropna(subset=["x","y"])
test = df[df.split=="test"].dropna(subset=["x","y"])

x_tr = train["x"].values
y_tr = train["y"].values
x_te = test["x"].values
y_te = test["y"].values

n_train = len(x_tr)
n_test = len(x_te)

slope, intercept = np.polyfit(x_tr, y_tr, 1)

pred = intercept + slope * x_te
resid = y_te - pred
test_rmse = float(np.sqrt(np.mean(resid**2)))
test_mae = float(np.mean(np.abs(resid)))

result = {
    "method": "train_only_ols_holdout",
    "n_train": int(n_train),
    "n_test": int(n_test),
    "intercept": float(intercept),
    "slope": float(slope),
    "test_rmse": test_rmse,
    "test_mae": test_mae
}

print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0061 · 9.1s

With skill

Both original and changed-input executions passed.

Numerical answer
{
  "method": "train_only_ols_holdout",
  "n_train": 24,
  "n_test": 12,
  "intercept": -0.25917168496622917,
  "slope": 1.0029666153104575,
  "test_rmse": 1.4837478210602169,
  "test_mae": 1.1711495352409773
}
Submitted code excerpt
import pandas as pd, numpy as np, json

df = pd.read_csv('data.csv')
train = df[df.split=='train']
test = df[df.split=='test']

x = train.x.values
y = train.y.values
slope, intercept = np.polyfit(x, y, 1)

pred = intercept + slope*test.x.values
resid = test.y.values - pred
rmse = np.sqrt(np.mean(resid**2))
mae = np.mean(np.abs(resid))

result = {
    "method": "train_only_ols_holdout",
    "n_train": int(len(train)),
    "n_test": int(len(test)),
    "intercept": float(intercept),
    "slope": float(slope),
    "test_rmse": float(rmse),
    "test_mae": float(mae)
}
print(json.dumps(result))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0148 · 15.5s

Recorded case analyze-norris-holdout · Run new-sources-and-analysis-v3. Output excerpts are bounded; the full submitted answer was scored.

Penguins: four simultaneous testspenguins.csvAI alonePassedWith checklistPassedWith skillPassed

Compare 2009 with 2007 for four measurements and apply Benjamini-Hochberg adjustment.

Source dataset ↗ · Source cluster: penguins

AI alone

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_bh_four_tests",
  "bill_length_mm": {
    "difference": -0.5150349650349781,
    "p_value": 0.3731192317432984,
    "q_value": 0.4657145087383828,
    "reject_fdr_05": false
  },
  "bill_depth_mm": {
    "difference": -0.7772727272727238,
    "p_value": 0.004087705017456997,
    "q_value": 0.01635082006982799,
    "reject_fdr_05": true
  },
  "flipper_length_mm": {
    "difference": 4.11888111888112,
    "p_value": 0.018363579336678298,
    "q_value": 0.036727158673356596,
    "reject_fdr_05": true
  },
  "body_mass_g": {
    "difference": -55.157342657342724,
    "p_value": 0.4657145087383828,
    "q_value": 0.4657145087383828,
    "reject_fdr_05": false
  }
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv('data.csv')
vars_ = ['bill_length_mm','bill_depth_mm','flipper_length_mm','body_mass_g']

results = {}
pvals = []
diffs = []
for v in vars_:
    g2009 = df.loc[df['year']==2009, v].dropna()
    g2007 = df.loc[df['year']==2007, v].dropna()
    t, p = stats.ttest_ind(g2009, g2007, equal_var=False)
    diff = g2009.mean() - g2007.mean()
    diffs.append(diff)
    pvals.append(p)

pvals = np.array(pvals)
n = len(pvals)
order = np.argsort(pvals)
ranked = pvals[order]
q = ranked * n / (np.arange(1, n+1))
# enforce monotonicity from the end
q_mono = np.minimum.accumulate(q[::-1])[::-1]
q_final = np.empty(n)
q_final[order] = np.minimum(q_mono, 1.0)

out = {"method": "welch_bh_four_tests"}
for i, v in enumerate(vars_):
    out[v] = {
        "difference": float(diffs[i]),
        "p_value": float(pvals[i]),
        "q_value": float(q_final[i]),
        "reject_fdr_05": bool(q_final[i] <= 0.05)
    }

print(json.dumps(out))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0100 · 8.9s

With checklist

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_bh_four_tests",
  "bill_length_mm": {
    "difference": -0.5150349650349781,
    "p_value": 0.3731192317432984,
    "q_value": 0.4657145087383828,
    "reject_fdr_05": false
  },
  "bill_depth_mm": {
    "difference": -0.7772727272727238,
    "p_value": 0.004087705017456997,
    "q_value": 0.01635082006982799,
    "reject_fdr_05": true
  },
  "flipper_length_mm": {
    "difference": 4.11888111888112,
    "p_value": 0.018363579336678298,
    "q_value": 0.036727158673356596,
    "reject_fdr_05": true
  },
  "body_mass_g": {
    "difference": -55.157342657342724,
    "p_value": 0.4657145087383828,
    "q_value": 0.4657145087383828,
    "reject_fdr_05": false
  }
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv("data.csv")

vars_ = ["bill_length_mm","bill_depth_mm","flipper_length_mm","body_mass_g"]
pvals = []
diffs = {}
for v in vars_:
    sub = df[["year", v]].dropna()
    g2009 = sub[sub.year==2009][v].values
    g2007 = sub[sub.year==2007][v].values
    diff = float(np.mean(g2009) - np.mean(g2007))
    t, p = stats.ttest_ind(g2009, g2007, equal_var=False)
    diffs[v] = diff
    pvals.append(float(p))

m = len(pvals)
order = np.argsort(pvals)
ranked = np.array(pvals)[order]
q = np.empty(m)
prev = 1.0
for i in range(m-1, -1, -1):
    rank = i+1
    val = ranked[i]*m/rank
    prev = min(prev, val)
    q[i] = prev
q_final = np.empty(m)
q_final[order] = q

out = {"method": "welch_bh_four_tests"}
for idx, v in enumerate(vars_):
    out[v] = {
        "difference": diffs[v],
        "p_value": pvals[idx],
        "q_value": float(q_final[idx]),
        "reject_fdr_05": bool(q_final[idx] <= 0.05)
    }

print(json.dumps(out))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0452 · 27.6s

With skill

Both original and changed-input executions passed.

Numerical answer
{
  "method": "welch_bh_four_tests",
  "bill_length_mm": {
    "difference": -0.5150349650349781,
    "p_value": 0.3731192317432984,
    "q_value": 0.4657145087383828,
    "reject_fdr_05": false
  },
  "bill_depth_mm": {
    "difference": -0.7772727272727238,
    "p_value": 0.004087705017456997,
    "q_value": 0.01635082006982799,
    "reject_fdr_05": true
  },
  "flipper_length_mm": {
    "difference": 4.11888111888112,
    "p_value": 0.018363579336678298,
    "q_value": 0.036727158673356596,
    "reject_fdr_05": true
  },
  "body_mass_g": {
    "difference": -55.157342657342724,
    "p_value": 0.4657145087383828,
    "q_value": 0.4657145087383828,
    "reject_fdr_05": false
  }
}
Submitted code excerpt
import pandas as pd, numpy as np
from scipy import stats
import json

df = pd.read_csv('data.csv')
vars_ = ['bill_length_mm','bill_depth_mm','flipper_length_mm','body_mass_g']
pvals = []
diffs = []
for v in vars_:
    g09 = df[df.year==2009][v].dropna()
    g07 = df[df.year==2007][v].dropna()
    t,p = stats.ttest_ind(g09, g07, equal_var=False)
    diff = g09.mean()-g07.mean()
    diffs.append(diff)
    pvals.append(p)

pvals = np.array(pvals)
n = len(pvals)
order = np.argsort(pvals)
ranked = pvals[order]
q = ranked * n / (np.arange(n)+1)
q = np.minimum.accumulate(q[::-1])[::-1]
qvals = np.empty(n)
qvals[order] = np.minimum(q,1)

out = {"method":"welch_bh_four_tests"}
for i,v in enumerate(vars_):
    out[v] = {
        "difference": float(diffs[i]),
        "p_value": float(pvals[i]),
        "q_value": float(qvals[i]),
        "reject_fdr_05": bool(qvals[i]<=0.05)
    }

print(json.dumps(out))

Original input: passed · Changed input: passed

Recorded score: 1/1 · $0.0218 · 14.1s

Recorded case analyze-penguin-multiplicity · Run new-sources-and-analysis-v3. Output excerpts are bounded; the full submitted answer was scored.

← Results overviewDownload experiment data