*! version 3.0-es 26jul2026

Titulo

mmd_2s -- Test de dos muestras Maximum Mean Discrepancy (MMD), ponderacion opcional

Sintaxis

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

Requerido
by(varname) variable de agrupamiento; debe tener exactamente 2 valores distintos

Opcional
boot(#) numero de replicas de bootstrap; por defecto boot(200)

reps(#) numero de corridas independientes que se promedian para el estadistico observado; por defecto reps(20)

seed(#) semilla del generador aleatorio (Mata rseed()); por defecto usa la semilla actual de Stata

boxplot dibuja un boxplot (NO ponderado) de varname por by(), anotado con p-boot, effect size y neff

kdensity dibuja la densidad kernel (ponderada, si corresponde) de ambos grupos

bw(#) ancho de banda para kdensity; por defecto usa el sigma del propio test

npoints(#) numero de puntos de grilla para kdensity; por defecto 200

Se permiten aweight, pweight e iweight; ver weight.
La ponderacion es totalmente opcional -- omiti [weight] para el test sin ponderar.

Descripcion

mmd_2s pone a prueba la hipotesis nula de que dos muestras
independientes provienen de la misma distribucion, usando un
estadistico Maximum Mean Discrepancy (MMD) con kernel RBF (gaussiano) y
un p-valor por bootstrap.

La version 3.0 fusiona los anteriores mmd_2s (sin ponderar) y
mmd_2s_pond (ponderado) en un solo comando, siguiendo el mismo
diseño que kstest.ado: la ponderacion se especifica con la
sintaxis nativa de Stata ([pweight/aweight/iweight]), no con
una opcion propia obligatoria. Si no se especifica ningun peso, se usa
un peso interno de 1 para cada observacion, y el comando reproduce
exactamente (a precision de maquina) el comportamiento del antiguo
mmd_2s v1.0 sin ponderar, con los mismos seed()/boot()/
reps(). Ya no existe un archivo mmd_2s_pond separado que
mantener sincronizado -- un solo codigo, un solo set de funciones Mata,
el peso vale 1 por defecto.

El estimador de base es el MMD de tiempo lineal por pares aleatorios
(estilo Gretton): se sortean pares al azar dentro de cada grupo, se
calcula una combinacion de 4 terminos con kernel RBF por par, y el
promedio (ponderado, si corresponde) entre pares es el estadistico,
acotado en 0 y promediado sobre tres anchos de banda (sigma/2, sigma,
2*sigma). Cuando se especifican pesos, quienes entran a cada par se
sigue eligiendo de forma uniforme al azar; el peso solo afecta
cuanto cuenta un par ya elegido (ver Metodos
y formulas). Es una eleccion de diseño deliberada, adecuada al
estimador por pares -- no es el mismo patron de normalizacion que
usa kstest (que pondera la ECDF completa de cada grupo por
su propio peso total). Los dos comandos ponderan de formas
estructuralmente distintas, apropiadas a su propio estadistico.

Opciones

by(varname) (requerido) especifica la variable de
agrupamiento. Debe tomar exactamente dos valores distintos en la
muestra de estimacion (tras excluir automaticamente los missing de
by()); el comando termina con error si no es asi.

weight (aweight/pweight/iweight) es
opcional. Si se especifica, los valores de peso no positivos o
missing se excluyen automaticamente de la muestra (mismo tratamiento
que un by() missing). Si se omite, cada observacion recibe peso 1.

boot(#) numero de replicas de bootstrap para armar
la distribucion nula (remuestreo con reemplazo del pool combinado; los
pesos viajan junto con cada fila). Por defecto 200; subilo para mas
resolucion en el p-valor (el minimo detectable es 1/(boot+1)).

reps(#) numero de corridas independientes que se
promedian para formar el estadistico observado. Por defecto
reps(20). Como el estimador submuestrea los datos, una sola
corrida (reps(1)) puede ser severamente ruidosa, sobre todo
con tamaño de muestra efectivo chico o moderado -- confirmado en la
practica: con n=52/22 (un tamaño de muestra comun en casos reales), una
sola corrida dio un p-valor no significativo (p=0.22) para una
diferencia que tanto un test de Kolmogorov-Smirnov completo como este
mismo estadistico MMD promediado sobre muchas corridas encontraron
altamente significativa (p<0.01). reps(1) desactiva por completo
el diagnostico de estabilidad (no hay nada contra que comparar una
corrida unica) y por eso imprime un aviso explicito. No bajes
reps() por debajo del default para ahorrar tiempo sin antes
revisar r(mmd_stat_sd) / el cv_draw y cv_se impresos con un
reps() mas alto.

seed(#) fija la semilla de Mata via rseed()
antes de cualquier remuestreo, para reproducibilidad. Si se omite, se
usa el estado aleatorio actual de Mata.

boxplot dibuja un graph box de varname segun by(),
anotado con el estadistico MMD (ponderado, si corresponde), p-boot,
effect size y neff de ambos grupos. El boxplot en si SIEMPRE sale
sin ponderar, sin importar que peso hayas especificado -- es una
limitacion del propio comando nativo graph box de Stata, que no
tiene ninguna opcion de ponderacion, no una decision de este paquete. El
estadistico y el p-valor de la anotacion SI estan correctamente
ponderados; solo el dibujo de cajas no. Si necesitas que el grafico en
si refleje la ponderacion, usa kdensity en su lugar, que calcula
la densidad kernel ponderada directamente (ver abajo) y no tiene esta
limitacion.

kdensity dibuja curvas de densidad kernel para ambos grupos,
ponderadas si se especifico un peso, con ancho de banda anclado al
sigma del propio test salvo que bw() lo sobreescriba.

Resultados guardados

Escalares

r(mmd_stat) estadistico MMD observado (promedio de reps() corridas)

r(mmd_stat_sd) desviacion estandar entre corridas de reps() (missing solo si se fuerza reps(1) explicitamente -- no es el default)

r(mmd_boot_mean) media de la distribucion nula de bootstrap

r(p_boot) p-valor de bootstrap, (#(boot>=stat)+1)/(boot()+1) (correccion "+1" de Davison-Hinkley -- nunca da exactamente 0)

r(effect_size) mmd_stat / mmd_boot_mean

r(nA) / r(nB) tamaño de muestra crudo, cada grupo

r(neff_A) / r(neff_B) tamaño de muestra efectivo de Kish (= n crudo si no se pondero)

r(sigma) ancho de banda del kernel RBF (heuristica de mediana)

r(N_boot) / r(N_reps) numero de replicas/corridas usadas

Macros

r(weighted) "yes" o "no"

Metodos y formulas

Tamaño de muestra efectivo: n_eff = (sum w)^2 / sum(w^2) (Kish, via
Monahan 2011, ec. 12.4.5). Se reduce al n crudo cuando w=1 para todos.

Estadistico MMD (una corrida): sea m = floor(min(neff_A,neff_B)/2).
Se sortean 2m indices de forma uniforme al azar (sin reemplazo) de
cada grupo, se parten en dos mitades, y se arma
h = k(x1,x2)+k(y1,y2)-k(x1,y2)-k(x2,y1) con kernel RBF. El estadistico
es el promedio ponderado de h usando pesos por par
w_par=(wx1*wx2+wy1*wy2)/2, acotado en 0. Ojo: esto pondera la
contribucion de cada par ya elegido de forma uniforme -- NO hace
que una observacion de peso alto tenga mas chance de ser sorteada para
un par. Se repite con anchos de banda sigma/2, sigma, 2*sigma y se
promedia, y eso se vuelve a promediar sobre reps().

Ancho de banda (sigma): mediana de las diferencias absolutas
pareadas sobre la variable combinada; si el pool tiene mas de 2000
observaciones, primero se sortea una submuestra de 2000 (ponderada, sin
reemplazo, Efraimidis-Spirakis) -- con 2000 o menos, el resultado no
depende de los pesos en absoluto.

Bootstrap: la muestra combinada (valores y pesos juntos) se
remuestrea con reemplazo boot() veces; se recalcula el mismo
estadistico cada vez. p-valor = (# replicas >= observado + 1) /
(boot() + 1).

Ejemplos

Preparacion

. sysuse auto, clear

Sin ponderar (no se especifica [weight]): ¿el rendimiento (mpg) tiene la misma distribucion en autos importados vs. nacionales?

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

Mismo ejemplo, con grafico de densidad kernel

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

Estabilizando el estadistico con reps() y el boxplot diagnostico

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

Preparacion, segunda base

. sysuse nlsw88, clear

Ponderado, con aweight: ¿el salario (wage) tiene la misma distribucion segun afiliacion sindical (union)?

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

Chequeo de reduccion al caso no ponderado (peso constante = 1)

. gen peso1 = 1

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

. * comparar r(mmd_stat), r(p_boot) y r(effect_size) contra la corrida sin [weight], misma seed()

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

Ejemplo con datos reales, microdatos publicos: superficie cosechada
entre cultivos transitorios y permanentes (P204_TIPO), usando la
Encuesta Nacional Agropecuaria (ENA) de Peru, ponderado por el factor de
expansion. Los datos se descargan directo del portal oficial de
microdatos de INEI con el comando complementario sriinei
(tambien disponible en SSC) -- este paquete no distribuye ningun archivo
de datos. Probado y validado por el autor con datos reales:

. 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

Autor / notas de desarrollo

La v3.0 de mmd_2s fusiona el antiguo par
mmd_2s/mmd_2s_pond en un solo comando, adoptando el diseño
de ponderacion opcional de kstest.ado (Ariel Linden). Ver
tambien sriinei, un comando complementario para descargar
microdatos publicos directo de un portal oficial de estadisticas
(usado en el ejemplo con datos reales de arriba).

Andres Talavera Cuya

Direccion Nacional de Censos y Encuestas -- INEI Peru

Email: atalaveracuya@gmail.com

Junio 2026

Licencia

Este modulo se distribuye bajo los terminos de la GPL v3
(https://www.gnu.org/licenses/gpl-3.0.txt).

Version

mmd_2s v3.0 -- 26jul2026