Estimation and maximum likelihood

import numpy as np, pandas as pd

df = pd.read_parquet("lumen_orders.parquet")
ov = df["order_value"].to_numpy()

true_var = ov.var(ddof=0)          # variance of the whole population ~ 38^2 ~ 1450
rng = np.random.default_rng(7)

biased, unbiased = [], []
for _ in range(5000):
    s = rng.choice(ov, size=30, replace=False)
    biased.append(np.var(s, ddof=0))    # divide by N   - the tempting one
    unbiased.append(np.var(s, ddof=1))  # divide by N-1 - Bessel

print("true variance      :", round(true_var, 1))
print("mean of /N   (ddof=0):", round(np.mean(biased), 1))    # LOW - underestimates
print("mean of /N-1 (ddof=1):", round(np.mean(unbiased), 1))  # matches true_var
# ddof=0 lands ~ (n-1)/n = 29/30 of the truth; ddof=1 is unbiased.

rng = np.random.default_rng(42)
converted = rng.random(ov.size) < 0.032   # the Lumen conversion flag, ~3.2% true

k = converted.sum()          # number of conversions
n = converted.size           # number of sessions
print("observed rate:", round(k / n, 4))   # ~ 0.032

# log-likelihood of Binomial(n, p) over a grid of candidate p
p_grid = np.linspace(0.001, 0.08, 800)
loglik = k * np.log(p_grid) + (n - k) * np.log(1 - p_grid)

p_hat = p_grid[np.argmax(loglik)]
print("MLE p-hat    :", round(p_hat, 4))    # ~ 0.032 - the peak
print("sample prop  :", round(k / n, 4))    # ~ 0.032 - identical
# The MLE of a proportion IS the sample proportion. The curve peaks at 3.2%.

from scipy import stats
from scipy.optimize import minimize

# (a) let scipy fit a lognormal to the right-skewed order values
shape, loc, scale = stats.lognorm.fit(ov, floc=0)   # floc=0 pins location
print("lognorm sigma-hat:", round(shape, 3))   # ~ 0.55 - the value b1 generated with
print("lognorm scale    :", round(scale, 1))   # ~ exp(mu), near the median

# (b) roll your own: minimize the NEGATIVE log-likelihood of a Normal on log(ov)
logs = np.log(ov)
def neg_loglik(params, data):
    mu, sigma = params
    if sigma <= 0:
        return np.inf
    return -np.sum(stats.norm.logpdf(data, loc=mu, scale=sigma))

res = minimize(neg_loglik, x0=[logs.mean(), 1.0], args=(logs,))
mu_hat, sigma_hat = res.x
print("mu-hat   :", round(mu_hat, 3))     # ~ 4.0  (matches log-mean)
print("sigma-hat:", round(sigma_hat, 3))  # ~ 0.55 (matches log-std, ddof=0)
# MLE of a Normal mean IS the sample mean; MLE of its variance divides by N.

rng = np.random.default_rng(42)
B = 2000
boot = np.empty(B)
for b in range(B):
    boot[b] = rng.choice(ov, size=ov.size, replace=True).mean()

ci = np.percentile(boot, [2.5, 97.5])
print("estimate      :", round(ov.mean(), 2))          # ~ 74
print("95% bootstrap CI:", [round(ci[0], 2), round(ci[1], 2)])  # ~ [73.7, 74.3]
# Same tool works for the MLE of any parameter - resample, refit, read the spread.