*! version 3.0 26jul2026

Title

mmd_2s -- Two-sample Maximum Mean Discrepancy (MMD) test, weight optional

Syntax

mmd_2s varname [if] [in] [weight],
by(varname)
[options]

Required
by(varname) grouping variable, must have exactly 2 distinct values

Optional
boot(#) number of bootstrap replications; default boot(200)

reps(#) number of independent draws averaged for the observed statistic; default reps(20)

seed(#) sets the random-number seed via Mata rseed(); default uses current Stata seed

boxplot draws an (unweighted) boxplot of varname by by() annotated with p-boot, effect size, and neff

kdensity draws a (weighted, if applicable) kernel density plot of both groups

bw(#) bandwidth for kdensity; default is the test's own sigma

npoints(#) number of grid points for kdensity; default 200

aweights, pweights, and iweights are allowed; see weight.
Weighting is entirely optional -- omit [weight] for the unweighted test.

Description

mmd_2s tests the null hypothesis that two independent samples were
drawn from the same distribution, using a Maximum Mean Discrepancy (MMD)
statistic with an RBF (Gaussian) kernel and a bootstrap p-value.

Version 3.0 merges the former mmd_2s (unweighted) and
mmd_2s_pond (weighted) into a single command, following the same
design as kstest.ado: weighting is specified with Stata's
native [pweight/aweight/iweight] syntax, not a separate
required option. When no weight is specified, an internal weight of 1
is used for every observation, and the command reproduces the old
unweighted mmd_2s v1.0 exactly (to machine precision, same
seed()/boot()/reps()). There is no longer a separate
mmd_2s_pond file to keep in sync -- one codebase, one set of Mata
functions, weight defaults to 1.

The underlying estimator is the linear-time / paired-subsample MMD
(Gretton-style): random pairs are drawn within each group, a four-term
RBF kernel combination is computed per pair, and the (optionally
weighted) average across pairs is the statistic, clipped at 0 and
averaged over three bandwidths (sigma/2, sigma, 2*sigma). When weights
are supplied, which observations are paired is still chosen
uniformly at random; the weight only affects how much a given pair
counts in the average (see Methods and formulas).
This is a deliberate design choice suited to the paired-subsample
estimator, and is not the same normalization pattern used by
kstest (which weights each group's full empirical CDF by
its own total weight) -- the two commands weight in structurally
different ways appropriate to their respective statistics.

Options

by(varname) (required) specifies the grouping
variable. It must take exactly two distinct values in the estimation
sample (after excluding missing by(), which is done automatically);
the command exits with an error otherwise.

weight (aweight/pweight/iweight) is
optional. If specified, non-positive or missing weight values are
excluded from the sample automatically (same treatment as missing
by()). If omitted, every observation gets weight 1.

boot(#) number of bootstrap replications used to
build the null distribution (resampling with replacement from the
pooled sample; weights travel with each row). Default 200; increase for
finer p-value resolution (minimum detectable p-value is 1/(boot+1)).

reps(#) number of independent draws averaged to form
the observed MMD statistic. Default reps(20). Because the
estimator subsamples the data, a single draw (reps(1)) can be
severely noisy, especially with small or moderate effective sample
size -- confirmed empirically: with n=52/22 (a common real-world sample
size), a single draw gave a non-significant p-value (p=0.22) for a
difference that both a full Kolmogorov-Smirnov test and this same MMD
statistic averaged over many draws found highly significant (p<0.01).
reps(1) disables the stability diagnostic entirely (there is
nothing to compare a single draw against) and prints an explicit warning
for that reason. Do not lower reps() below the default to save
time without first checking r(mmd_stat_sd) / the printed cv_draw
and cv_se at a higher reps().

seed(#) sets the Mata seed via rseed() before
any resampling, for reproducibility. If omitted, the current Mata random
state is used.

boxplot draws a graph box of varname over by(),
annotated with the (weighted, if applicable) MMD statistic, p-boot,
effect size, and neff of both groups. The boxplot itself is always
unweighted, regardless of any [weight] specified -- this is a
limitation of Stata's native graph box command, which has no
weighting option at all, not a choice made by this package. The
statistic and p-value in the annotation ARE correctly weighted; only the
box-and-whisker drawing is not. If you need the plot itself to reflect
the weighting, use kdensity instead, which computes a weighted
kernel density directly (see below) and is not subject to this
limitation.

kdensity draws kernel density curves for both groups, weighted if
a weight was specified, with bandwidth anchored to the test's own sigma
unless bw() overrides it.

Stored results

Scalars

r(mmd_stat) observed MMD statistic (average of reps() draws)

r(mmd_stat_sd) SD across reps() draws (missing only if reps(1) is explicitly forced -- not the default)

r(mmd_boot_mean) mean of the bootstrap null distribution

r(p_boot) bootstrap p-value, (#(boot>=stat)+1)/(boot()+1) (Davison-Hinkley "+1" correction -- never exactly 0)

r(effect_size) mmd_stat / mmd_boot_mean

r(nA) / r(nB) raw sample size, each group

r(neff_A) / r(neff_B) Kish effective sample size (= raw n if unweighted)

r(sigma) RBF kernel bandwidth (median heuristic)

r(N_boot) / r(N_reps) number of replications/draws used

Macros

r(weighted) "yes" or "no"

Methods and formulas

Effective sample size: n_eff = (sum w)^2 / sum(w^2) (Kish, via
Monahan 2011, eq. 12.4.5). Reduces to raw n when w=1 for all.

MMD statistic (single draw): let m = floor(min(neff_A,neff_B)/2).
Draw 2m indices uniformly at random (without replacement) from each
group, split into two halves, form h = k(x1,x2)+k(y1,y2)-k(x1,y2)-k(x2,y1)
with an RBF kernel. The statistic is the weighted average of h
using per-pair weights w_par=(wx1*wx2+wy1*wy2)/2, floored at 0. Note
this weights the contribution of each already-uniformly-selected
pair -- it does not make high-weight observations more likely to be
selected into a pair. Repeated at bandwidths sigma/2, sigma, 2*sigma and
averaged, then averaged again over reps().

Bandwidth (sigma): median of pairwise absolute differences on the
pooled variable; if pooled n>2000, a subsample of 2000 is drawn first
(weighted, without replacement, Efraimidis-Spirakis) -- with n<=2000 the
result does not depend on weights at all.

Bootstrap: pooled sample (values and weights together) resampled
with replacement boot() times; same statistic recomputed each
time. p-value = (# replicates >= observed + 1) / (boot() + 1).

Examples

Setup

. sysuse auto, clear

Unweighted (no [weight] specified): does fuel economy (mpg) have the same distribution in imported vs. domestic cars?

. mmd_2s mpg, by(foreign) boot(200) reps(20) seed(12345)

Same example, with kernel density plot

. mmd_2s mpg, by(foreign) boot(200) reps(20) seed(12345) kdensity

Stabilizing the statistic with reps() and the diagnostic boxplot

. mmd_2s price, by(foreign) boot(200) reps(20) seed(12345) boxplot

Setup, second dataset

. sysuse nlsw88, clear

Weighted, with aweight: does wage have the same distribution by union affiliation?

. mmd_2s wage [aweight=hours], by(union) boot(500) reps(20) seed(12345) boxplot

Check that it reduces to the unweighted case (constant weight = 1)

. gen peso1 = 1

. mmd_2s wage [aweight=peso1], by(union) boot(500) reps(20) seed(12345)

. * compare r(mmd_stat), r(p_boot), and r(effect_size) against the run without [weight], same seed()

. mmd_2s wage, by(union) boot(500) reps(20) seed(12345)

Real-data example, weighted, with public microdata: harvested
surface area, transitory vs. permanent crops, using Peru's national
agricultural survey (ENA). Data downloaded directly from the official
INEI microdata portal with the companion command sriinei
(also available on SSC) -- no data files are distributed with this
package.

. sriinei, codigo(1036) modulo(1895) tipo(stata) destino("C:\BD_INEI\mic")

. cd "C:\BD_INEI\mic\1036-Modulo1895"

. use 03_CAP200AB, clear

. keep if CODIGO==1 & inlist(P204_TIPO,1,2)

. gen SUP_COSECHADA=P217_SUP_ha

. drop if missing(SUP_COSECHADA, FACTOR_PRODUCTOR)

. gen double SUP_COSECHADA_log = log(SUP_COSECHADA + 1)

. mmd_2s SUP_COSECHADA_log [aweight=FACTOR_PRODUCTOR], by(P204_TIPO) boot(30) reps(20) kdensity

Author / development notes

mmd_2s v3.0 merges the former mmd_2s/mmd_2s_pond pair
into a single command, adopting the weight-optional design of
kstest.ado (Ariel Linden). See also
sriinei, a companion command for downloading public
microdata directly from an official statistics portal (used in the
real-data example above).

Andres Talavera Cuya

Direccion Nacional de Censos y Encuestas -- INEI Peru

Email: atalaveracuya@gmail.com

Junio 2026

License

This module is made available under the terms of the GPL v3
(https://www.gnu.org/licenses/gpl-3.0.txt).

Version

mmd_2s v3.0 -- 26jul2026