标签：#statistics

Short version of a post on gex.live/research ; the full write-up, definitions and reproduce block live there. Most published dealer-gamma numbers are built from open interest : yesterday's outstanding contracts, multiplied by a convention about who holds which side. Whether the convention is right is a separate question. The prior question is simpler: how much of what trades today was already in that book this morning — and how much of tomorrow's book is being created today? Open interest and volume are enough to answer it, with no assumption about who bought. Sample: SPX and SPXW, 2022-04-14 to 2026-08-14, 1,081 trading days, every expiry within about a month (0DTE plus the 21 nearest), 8.6 million contract-days, 4.3 million with volume. Definitions Per contract (expiry, strike, right) and session D: OI(D) is open interest at the start of D, OI(D+1) at the start of the next session, ΔOI = OI(D+1) − OI(D) , vol the day's volume in that contract. |ΔOI| / vol is a lower bound on how one-sided the day's trading in that contract was — 1.0 means every lot opened (or every lot closed), 0 means opens and closes cancelled. Contracts expiring on D have no next-day OI and drop out of the ΔOI statistics; 4.1% of rows (3.8% of volume) show |ΔOI| > vol, which is impossible (OI snapshot timing) and are excluded. The book grows by 40% of what trades, every day days to expiry on D net ΔOI / volume |ΔOI| / volume (lower bound on one-sidedness) share of volume in contracts whose OI rose contract-days 1–5 37.8% 41.8% 90.3% 813,013 6–21 42.7% 53.3% 81.1% 2,206,446 22+ 42.8% 57.6% 76.5% 831,896 Across the whole book, net ΔOI is 39.9% of the day's volume on the median session (IQR 36.2–44.0%), positive in every year and every expiry bucket: the SPX book is always being built faster than it is unwound, until expiry does the unwinding. Far expiries are open-and-hold (a day's trading in a 22+ DTE contract is at least 58% one-sided); the nearest expiries churn (42% at 1–5 DTE). Per contract-

Every A/B testing tutorial ends the same way: run the test, wait for significance, ship the winner. Then you run a real test and variant B converts 12% better on newsletter signups, brings in 4% less revenue per visitor, and bounce is flat. Nothing is significant except the signups. Ship it? I spent an embarrassing amount of time on this question while building an A/B engine, and most of what I read online didn't help, because most of it assumes one metric. This post is what I ended up with. It's not novel — the statistics are decades old — but I couldn't find it written down in one place with working code, so here it is. Why the p-value doesn't answer the question you're asking Two problems, and the second one is the bad one. Multiple comparisons. Three metrics at α = 0.05 means roughly a 14% chance of at least one false positive if nothing is actually different. Bonferroni fixes this, but now you need α = 0.017 per metric and your test needs to run three times as long. On a site doing 300 conversions a month that's not a fix, it's a refusal. The p-value is answering a different question. It tells you the probability of your data assuming no difference exists. What you actually want to know is: if I ship B, how much do I expect to lose if I'm wrong? Those are not the same question and no amount of Bonferroni turns one into the other. There's also the peeking problem — everyone checks the dashboard daily and stops when it goes green, which quietly inflates the false positive rate well past whatever α you wrote down. I'll come back to that, because Bayesian methods do not magically solve it, whatever you may have read. Posterior first, decision second For a conversion rate, the Beta-Binomial conjugate pair gives you the posterior in one line. With a uniform prior, after c conversions out of n visitors: p | data ~ Beta ( 1 + c , 1 + n - c ) That's it. No closed-form comparison between two Betas that's worth implementing, so sample. PHP has no Beta sampler in core, and

Non-Gaussian Distributions Explained from First Principles (Beginner Friendly) As we all know, the real-world dataset is not normalized , but most of us thought every dataset followed the famous bell curve . After all, everyone talks about the Normal Distribution . But then I looked at real-world datasets like: Income of people Stock market returns Website traffic YouTube views Population of cities None of them looked like a bell curve. That's when I realized something important. Not every dataset in the real world is normally distributed. In this article, we'll understand Non-Gaussian (Non-Normal) Distributions from first principles using simple language, intuition, and real-world examples. First, What Does "Non-Gaussian" Mean? The Normal Distribution (also called the Gaussian Distribution) has a very specific shape. It is: Bell-shaped Symmetrical Mean = Median = Mode Most observations lie near the average But what if our data doesn't look like that? Then it is called a Non-Gaussian Distribution . In simple words, Any probability distribution that does not follow the Normal Distribution is called a Non-Gaussian Distribution. Why Should We Care? Imagine you are analyzing the salaries of employees. Most employees earn between ₹25,000 and ₹1,00,000. But a few CEOs earn ₹50 lakh or even ₹2 crore. Will this data form a perfect bell curve? No. The extremely high salaries pull the distribution toward one side. If we wrongly assume the data is normal, our analysis can become misleading. That's why understanding Non-Gaussian Distributions is extremely important in Data Science. Before Learning Other Distributions... Let's understand two important ideas. These help us decide whether our data is normally distributed or not. Kurtosis — How Heavy Are the Tails? When beginners hear the word Kurtosis , they usually think it measures how tall the peak of a graph is. That's actually a common misconception. A better way to think about Kurtosis is this: How likely is the distribution

Every binomial confidence interval you have ever computed on an eval pass rate, Wald, Wilson, Clopper-Pearson, all of them, rests on one assumption: each example is an independent draw. Most eval sets violate it. You have 40 questions generated from the same 8 documents, or 200 turns from the same 30 conversations, or 150 examples that are really 50 cases with 3 paraphrases each. Those are not 200 independent observations. And when you feed a correlated set into a formula that assumes independence, the interval comes out too narrow, which means you declare differences significant that aren't. I want to walk through why, put a number on how much it matters, and show the fix, because this one is invisible: the code runs, the interval prints, and it is quietly wrong. Why clustering shrinks your real sample size Independent examples each carry their own information. Correlated examples carry overlapping information. If five questions come from the same document, and the model either understands that document or doesn't, those five outcomes move together. You did not learn five independent things about the model. You learned something closer to one and a half. The survey-statistics name for this is the design effect (Kish, "Survey Sampling," 1965). For clustered data it is approximately: Deff = 1 + (m̄ - 1) · ICC where m̄ is the average cluster size and ICC is the intra-cluster correlation, the fraction of total variance that lives between clusters rather than within them. Your effective sample size is: n_eff = n / Deff That is the number of independent examples your clustered set is actually worth. The number Take a realistic eval set: n = 200 examples, drawn from 40 source documents, so average cluster size m̄ = 5. Suppose the ICC is 0.3, which is unremarkable for "questions from the same document" (I have measured higher). Deff = 1 + (5 - 1) · 0.3 = 2.2 n_eff = 200 / 2.2 ≈ 91 Your 200-example eval is worth about 91 independent examples. The correct confidence interval

Most experiments that come back "no clear winner" were unwinnable on the day they launched. The data could not resolve an effect that size, and no amount of extra runtime was going to change that. You can find this out in about two minutes, before you spend anything, with one formula and a resampling pass over your own data. Here is the check, in three steps. Step 1. Compute the smallest lift your data can see For a two-arm test on a conversion rate, the smallest lift detectable at 95% confidence and 80% power is a one-liner: from math import sqrt Z_ALPHA = 1.96 # two-sided 95% Z_BETA = 0.84 # 80% power def mde ( baseline_cvr : float , n_per_arm : int ) -> tuple [ float , float ]: """ Minimum detectable effect: absolute (pp) and relative (%). """ se = sqrt ( 2 * baseline_cvr * ( 1 - baseline_cvr ) / n_per_arm ) abs_lift = ( Z_ALPHA + Z_BETA ) * se return abs_lift * 100 , abs_lift / baseline_cvr * 100 At a 3% conversion rate: clicks per arm smallest lift you can detect 5,000 +32% relative 20,000 +16% relative 100,000 +7% relative Read the middle row twice. Twenty thousand clicks per arm is a serious amount of traffic for a mid-market account, and a real 15% improvement still lands inside the confidence interval. The report will say "inconclusive," and the team will read that as a verdict on the idea. It is a verdict on the instrument. Invert the same formula and the planning question gets easier: at 3% baseline, detecting a 10% lift needs about 51,000 clicks per arm, and detecting a 5% lift needs about 203,000. If your account produces 8,000 clicks a month, you now know the honest answer to "how long should we run this." Step 2. Stop assuming your conversions are independent The formula above treats every click as an independent coin flip with the same probability. Account data does not behave that way, and the gap is not small. In a corpus of 31 advertiser accounts I maintain for diagnostic work (9.46 million search term rows, roughly $133M of spend, September 2024

Adapted from an appendix of my MS thesis. Markov Chain Monte Carlo Almost as soon as computers were invented, they were used for simulation. Markov chain Monte Carlo (MCMC) was invested as Los Alamos, Metropolis et al (1953) simulated a liquid in equilibrium with its gas phase. Their tour de force was the realization that they did not need to simulate the exact dynamics, they only needed to simulate some Markov chain with the same equilibrium distribution. The Metropolis algorithm was widely used by chemists and physicists, but was not widely known among statisticians until after 1990. Hastings (1970) generalized the Metropolis algorithm, and simulations following his scheme are said to use the Metropolis-Hastings (MH) algorithm [1]. A special case of the MH algorithm was introduced by Geman et al (1984) discussing optimization to find the posterior mode rather than simulation. Algorithms following their scheme are said to use the Gibbs sampler. It took some time for the spatial statistics community to understand that the Gibbs sampler simulated the posterior distribution, thus enabling full Bayesian inference of all kinds. Gelfand et al (1990) made the wider Bayesian community aware of the Gibbs sampler, and then it was rapidly realized that most Bayesian inference could be done using MCMC, whereas very little could be done without MCMC. Green (1995) generalized the MH algorithm as much as it could be generalized [1]. Theoretical Foundations A sequence X 1 ​ , X 2 ​ , … of random elements of some set is a Markov chain if the conditional distribution of X n + 1 ​ given X 1 ​ , … , X n ​ depends on X n ​ only. The set in which the X i ​ take values is called the state space of the Markov chain. A Markov chain has stationary transition probabilities if the conditional distribution of X n + 1 ​ given X n ​ does not depend on n . This is the main kind of Markov chain of interest in MCMC. The joint distribution of a Markov chain is determined by the following [1]. The ma

A friend messaged me a photo of a sealed booster box last month with one question: "worth it?" He'd already decided, really. The chase card in that set was all over his feed, so the box felt like a good deal. I asked him to send me the pull rates instead of the hype, and we spent twenty minutes turning "worth it?" into something we could actually compute. That exercise is a small, self-contained data problem. It's also a good example of how a clean-looking model can hand you a confident number that doesn't survive contact with reality. If you like building little estimators, this one is worth doing carefully, because the interesting part isn't the formula. It's everything the formula assumes. The formula is the easy part Expected value of a box is a weighted sum. Each card you can pull has a probability and a market value, and you multiply the two across every slot the box gives you. That's it. Undergrad probability. Here's a stripped-down version for a hypothetical set. I'm using made-up numbers so nobody mistakes this for real pull data — the point is the shape of the computation, not the specific set. # One "hit slot" in a box: probabilities cover the full outcome space. # Values are illustrative market estimates in USD. hit_table = [ { " name " : " Alt-art chase " , " p " : 0.0125 , " value " : 180.00 }, { " name " : " Secret rare " , " p " : 0.030 , " value " : 45.00 }, { " name " : " Full-art rare " , " p " : 0.100 , " value " : 8.00 }, { " name " : " Standard hit " , " p " : 0.400 , " value " : 0.55 }, { " name " : " No notable hit " , " p " : 0.4575 , " value " : 0.06 }, ] assert abs ( sum ( row [ " p " ] for row in hit_table ) - 1.0 ) < 1e-9 ev_per_slot = sum ( row [ " p " ] * row [ " value " ] for row in hit_table ) hit_slots_per_box = 36 # e.g. one meaningful slot per pack ev_box = ev_per_slot * hit_slots_per_box print ( f " EV per slot: $ { ev_per_slot : . 2 f } " ) # $4.65 print ( f " EV per box: $ { ev_box : . 2 f } " ) # $167.31 The box costs $150 sea

I recently calibrated a recovery-rate model that had only two weak features. Its point accuracy was almost nothing — R² basically zero. I expected its uncertainty estimates to be junk too. They weren't: the 90% conformal prediction intervals covered ~89% of held-out outcomes. Valid, just wide . That surprised me enough to nail it down, because it contradicts a belief a lot of us carry around: "my model isn't accurate, so I can't trust its uncertainty." For split conformal prediction, that's backwards. Here's the precise statement, a runnable demo, and the one caveat that actually bites. Coverage is a property of the procedure, not the model Split conformal prediction gives a distribution-free, finite-sample marginal coverage guarantee : P( Y ∈ Ĉ(X) ) ≥ 1 − α and it holds for any point model, as long as the calibration and test data are exchangeable. The model is a black box. You fit it however you like, then on a held-out calibration set you take the (1−α) quantile of the absolute residuals, and that quantile becomes the half-width of your intervals. Nowhere does that construction require the model to be good. A bad model just has large residuals, so the calibration quantile is large, so the intervals are wide — wide enough to still cover at the stated rate. Accuracy doesn't buy you validity ; it buys you efficiency (narrower intervals at the same coverage). The demo (numbers are reproducible, seed fixed) Same dataset and target, three models from strong to useless, target coverage 90%: model R² marginal coverage mean interval width gradient boosting 0.741 0.895 5.39 weak linear (1 noisy feature) 0.061 0.905 10.39 predict-the-mean −0.000 0.907 10.83 All three land at ~90% coverage. The only thing that changes is width: the good model's intervals are half as wide . That's the whole story in one table — validity is constant, efficiency tracks accuracy. import numpy as np from sklearn.linear_model import LinearRegression from sklearn.ensemble import GradientBoostingReg

In [1], the authors present a method for constructing a symbolic (nominal) representation for real-valued time series data. A symbolic representation is desirable because then it becomes possible to use many of the effective algorithms that require symbolic representation, like hashing and Markov models. The authors claim that one of the most useful time series operations is measuring the similarity between two time series data sets. To do this on the original time series, the Euclidean distance formula can be used. Therefore, for a time series transformation to be useful, distance measures applied to the corresponding transformations should provide some guaranteed lower bound on the true distance. This is a basic requirement for almost all time series algorithms in data mining. Non-symbolic transformations like Discrete Fourier Transform (DFT) and Piecewise Aggregate Approximation (PAA) models have this lower-bounding property. However, the authors claim no previously proposed symbolic representations do, which limits their usefulness. Additionally, the authors observe that most raw time series data sets have very high dimensionality. This is problematic because time series mining algorithms are $\mathcal{O}(cn)$, where n is the number of dimensions. Therefore, preferably any transformations on the original time series will reduce the dimensionality to a more manageable size. Unfortunately, the authors observe, previously proposed symbolic representations preserve the original time series dimensionality. Next, the authors present their symbolic representation, SAX (Symbolic Aggregate approXimation), which addresses each of the previously mentioned shortcomings of symbolic representations. SAX is unique in that it uses an intermediate transformation, PAA, and then nominalizes the PAA representation into a sequence of characters'a string. By using the intermediate PAA representation, SAX enjoys two benefits: It is able to exploit the dimensionality reducing propertie

I put my coursework from SIUe's Master's in Mathematics program up on the problem sets section of this site. Five courses from 2021-2022 that formed the core of my graduate statistics work. The Courses STAT 478 - Time Series Analysis Spring 2021, Dr. Beidi ARIMA models, forecasting, spectral analysis, state space methods. The problem sets mix theoretical derivations with R implementations for temporal data. STAT 482 - Regression Analysis Fall 2022, Dr. Andrew Neath Linear models, diagnostics, variable selection, model comparison. Neath's approach put the theoretical foundations front and center, with application following from understanding. STAT 575 - Computational Statistics Summer 2021, Dr. Qiang Beidi Probably the most hands-on course of the five. Topics included: Newton-Raphson and numerical optimization Monte Carlo simulation Sampling methods (inverse transform, acceptance-rejection) Hand-coded MLE for Poisson regression Implementing these algorithms from scratch instead of calling library functions teaches you what the methods are actually doing. You hit the edge cases. You debug convergence failures. That's where the understanding comes from. STAT 579 - Discrete Multivariate Analysis Spring 2021, Dr. Andrew Neath Categorical data analysis, log-linear models, contingency tables. Both the mathematical theory and R implementations for discrete multivariate data. STAT 581 - Statistical Methods Fall 2021, Dr. Neath Experimental design, ANOVA, general linear models with practical applications. Why Share This? A few reasons. Learning resource. Worked solutions for graduate statistics are surprisingly hard to find online. If someone studying this material stumbles across these and they help, good. Personal archive. I did much of this work during cancer treatment. Keeping it organized and accessible matters to me. Reference. I still look up my own derivations and implementations when something comes up in research. Easier to find them here than to dig through old dir

When can you safely use a simpler model for a series system? I ran extensive simulation studies with likelihood ratio tests to get a quantitative answer. The Problem In series system reliability, you estimate component parameters from masked failure data. For Weibull components, that means estimating (2m) parameters: shape (k_j) and scale (\lambda_j) for each of (m) components. But what if the components have similar failure characteristics? A reduced model with homogeneous shape parameters uses only (m+1) parameters (one common (k) plus (m) scales). This roughly halves the parameter count and has a nice property: the system itself becomes Weibull-distributed. The question is when this simplification is justified. Key Findings Robustness of the Reduced Model For well-designed series systems (components with similar failure characteristics), the result is striking: The reduced homogeneous-shape model cannot be rejected even with sample sizes approaching 30,000, far larger than anything typically available in practice. With realistic sample sizes (50 to 500), the likelihood ratio test shows no evidence against the reduced model when components truly have similar shapes. This is strong justification for using the simpler model. Sharp Boundaries The paper pins down exactly how much heterogeneity it takes to trigger rejection: Shape Deviation Sample Size LRT Decision 0.25 30,000 Fail to reject 0.50 1,000+ Reject 1.0 100+ Strong reject 3.0 50+ Very strong reject Even modest deviations in a single component's shape parameter provide evidence against the reduced model. The boundaries are clean. Practical Guidance Use the reduced model when: Components come from similar manufacturing processes Historical data suggests similar wear-out patterns Sample sizes are moderate ((n < 500)) You need a quick reliability assessment Use the full model when: Components have fundamentally different failure modes (infant mortality vs wear-out) Large samples are available ((n > 1000)) Precis