¿Cómo simular datos que satisfagan restricciones específicas, como tener una media específica y una desviación estándar?

56

Esta pregunta está motivada por mi pregunta sobre el metanálisis . Pero imagino que también sería útil para enseñar contextos en los que desea crear un conjunto de datos que refleje exactamente un conjunto de datos publicado existente.

Sé cómo generar datos aleatorios a partir de una distribución dada. Entonces, por ejemplo, si leo sobre los resultados de un estudio que tuvo:

una media de 102,
una desviación estándar de 5.2, y
un tamaño de muestra de 72.

Podría generar datos similares usando rnormen R. Por ejemplo,

set.seed(1234)
x <- rnorm(n=72, mean=102, sd=5.2)

Por supuesto, la media y la desviación estándar no serían exactamente iguales a 102 y 5.2 respectivamente:

round(c(n=length(x), mean=mean(x), sd=sd(x)), 2)
##     n   mean     sd 
## 72.00 100.58   5.25

En general, estoy interesado en cómo simular datos que satisfagan un conjunto de restricciones. En el caso anterior, las constantes son el tamaño de la muestra, la media y la desviación estándar. En otros casos, puede haber restricciones adicionales. Por ejemplo,

puede conocerse un mínimo y un máximo en los datos o en la variable subyacente.
se sabe que la variable toma solo valores enteros o solo valores no negativos.
los datos pueden incluir múltiples variables con interrelaciones conocidas.

Preguntas

En general, ¿cómo puedo simular datos que satisfagan exactamente un conjunto de restricciones?
¿Hay artículos escritos sobre esto? ¿Hay algún programa en R que haga esto?
Por el bien de ejemplo, ¿cómo podría y debería simular una variable para que tenga una media y un SD específicos?

r dataset simulation random-generation Jeromy Anglim
fuente

1

¿Por qué quieres que sean exactamente como los resultados publicados? ¿No son estas estimaciones de la media de la población y la desviación estándar dada su muestra de datos? Dada la incertidumbre en esas estimaciones, ¿quién puede decir que la muestra que muestra arriba no es consistente con sus observaciones?

Restablece a Monica - G. Simpson el

44

Debido a que esta pregunta parece estar recopilando respuestas que no dan en el blanco (en mi humilde opinión), me gustaría señalar que conceptualmente la respuesta es sencilla: las restricciones de igualdad se tratan como distribuciones marginales y las restricciones de desigualdad son análogos multivariados de truncamiento. El truncamiento es relativamente fácil de manejar (a menudo con muestreo de rechazo); El problema más difícil es encontrar una manera de muestrear estas distribuciones marginales. Esto significa tomar muestras marginales dada la distribución y la restricción, o integrarse para encontrar la distribución marginal y el muestreo a partir de ellas.

whuber

44

Por cierto, la última pregunta es trivial para las familias de distribución a escala de ubicación. Por ejemplo, x<-rnorm(72);x<-5.2*(x-mean(x))/sd(x)+102hace el truco.

whuber

1

@whuber, como alude el cardenal en un comentario a mi respuesta (que menciona este "truco") y un comentario a otra respuesta: este método, en general, no mantendrá las variables dentro de la misma familia de distribución, ya que está dividiendo por la desviación estándar de la muestra.

Macro

55

@Macro Este es un buen punto, pero quizás la mejor respuesta es "¡por supuesto que no tendrán la misma distribución"! La distribución que desea es la distribución condicional a las restricciones. En general, eso no será de la misma familia que la distribución de los padres. Por ejemplo, cada elemento de una muestra de tamaño 4 con media 0 y DE 1 extraída de una distribución normal tendrá una probabilidad casi uniforme en [-1.5, 1.5], porque las condiciones colocan límites superior e inferior en los valores posibles.

whuber

26

En general, para que la media y la varianza de su muestra sean exactamente iguales a un valor especificado previamente, puede desplazar y escalar adecuadamente la variable. En concreto, si es una muestra, entonces las nuevas variables $X_1, X_2, ..., X_n$

Z_{i} = \sqrt{c_{1}} (\frac{X_{i} - \bar{X}}{s_{X}}) + c_{2}

$Z_i = \sqrt{c_{1}} \left( \frac{X_i-\overline{X}}{s_{X}} \right) + c_{2}$

donde es la media muestral yes la varianza muestral de tal manera que la media muestral de loses exactamentey su varianza muestral es exactamente. Un ejemplo construido de manera similar puede restringir el rango: $\overline{X} = \frac{1}{n} \sum_{i=1}^{n} X_i$ $s^{2}_{X} = \frac{1}{n-1} \sum_{i=1}^{n} (X_i - \overline{X})^2$ $Z_{i}$ $c_2$ $c_1$

B_{i} = a + (b - a) (\frac{X_{i} - min ({X_{1}, . . ., X_{n}})}{max ({X_{1}, . . ., X_{n}}) - min ({X_{1}, . . ., X_{n}})})

$B_i = a + (b-a) \left( \frac{ X_i - \min (\{X_1, ..., X_n\}) }{\max (\{X_1, ..., X_n\}) - \min (\{X_1, ..., X_n\}) } \right)$

producirá un conjunto de datos que está restringido al intervalo . $B_1, ..., B_n$ $(a,b)$

Nota: Estos tipos de desplazamiento / escalado, en general, cambiarán la familia de distribución de los datos, incluso si los datos originales provienen de una familia de escala de ubicación.

Dentro del contexto de la distribución normal, la mvrnormfunción en le R permite simular datos normales (o multivariados normales) con una media / covarianza muestra especificada previamente mediante la configuración empirical=TRUE. Específicamente, esta función simula datos de la distribución condicional de una variable normalmente distribuida, dada la media muestral y la (co) varianza es igual a un valor preespecificado . Tenga en cuenta que las distribuciones marginales resultantes no son normales, como señaló @whuber en un comentario a la pregunta principal.

Aquí hay un ejemplo univariado simple donde la media muestral (de una muestra de ) está restringida a 0 y la desviación estándar muestral es 1. Podemos ver que el primer elemento es mucho más similar a una distribución uniforme que una normal distribución: $n=4$

library(MASS)
 z = rep(0,10000)
for(i in 1:10000)
{
    x = mvrnorm(n = 4, rep(0,1), 1, tol = 1e-6, empirical = TRUE)
    z[i] = x[1]
}
hist(z, col="blue")

$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$ ingrese la descripción de la imagen aquí

Macro
fuente

1

El no se distribuirá normalmente, aunque puede ser aproximadamente así si el tamaño de la muestra es grande. El primer comentario a la respuesta de @ Sean alude a esto.

Z_{i}

$Z_i$

cardenal

1

Bueno, es algo muy natural querer hacer ... y muchas veces no causa demasiados problemas.

cardenal

1

+1. En el ejemplo, el uniforme es la respuesta exacta , por cierto. (La aparente caída al final de la trama es un artefacto de cómo R dibuja los histogramas.)

whuber

1

@whuber, gracias por motivar este ejemplo. Dado que las distribuciones marginales cambian una vez que se condiciona la media / varianza de la muestra, parece que la mejor "respuesta" en el espíritu de la pregunta del OP es simplemente simular datos con media / varianza de la población igual a la informada como muestra cantidades (como lo sugirió el propio OP), ¿no? De esa manera, obtienes cantidades de muestra "similares" a las deseadas, y las distribuciones marginales son lo que querías que fueran.

Macro

1

@whuber, si su muestra es normal, entonces tiene una distribución , ¿sí? La "nueva" variable en cuestión será solo una combinación lineal de .

T_{i} = (X_{i} - \bar{X}) / s

$T_i = (X_i-\overline{X})/s$

t

$t$

T_{i}

$T_i$

Macro

22

En cuanto a su solicitud de documentos, hay:

Chatterjee, S. y Firat, A. (2007). Generación de datos con estadísticas idénticas pero gráficos diferentes: un seguimiento del conjunto de datos Anscombe. The American Statistician, 61 , 3, págs. 248-254.

Esto no es exactamente lo que está buscando, pero podría servir como grano para el molino.

Hay otra estrategia que nadie parece haber mencionado. Es posible generar datos (pseudo) aleatorios a partir de un conjunto de tamaño modo que el conjunto completo cumpla con las restricciones siempre que los datos restantes se fijen en los valores apropiados. Los valores requeridos deben poder resolverse con un sistema de ecuaciones, álgebra y algo de grasa en el codo. $N-k$ $N$ $k$ $k$ $k$

Por ejemplo, para generar un conjunto de datos a partir de una distribución normal que tendrá una media de muestra dada, , y una varianza, , deberá fijar los valores de dos puntos: y . Como la media muestral es: debe ser: La varianza muestral es: tanto (después de sustituir lo anterior por , frustrar / distribuir y reorganizar ... ) obtenemos: $N$ $\bar x$ $s^2$ $y$ $z$

\bar{x} = \frac{\sum_{i = 1}^{N - 2} x_{i} + y + z}{N}

$\bar x = \frac{\sum_{i=1}^{N-2}x_i\; + \;y\!+\!z}{N}$

y

$y$

y = N \bar{x} - (\sum_{i = 1}^{N - 2} x_{i} + z)

$y = N\bar x\; - \;\left(\sum_{i=1}^{N-2}x_i\!+\!z\right)$

s^{2} = \frac{\sum_{i = 1}^{N - 2} (x_{i} - \bar{x})^{2} + (y - \bar{x})^{2} + (z - \bar{x})^{2}}{N - 1}

$s^2 = \frac{\sum_{i=1}^{N-2}(x_i - \bar x)^2\; + \;(y - \bar x)^2\!+\!(z - \bar x)^2}{N-1}$

y

$y$

2 (N \bar{x} - \sum_{i = 1}^{N - 2} x_{i}) z - 2 z^{2} = N {\bar{x}}^{2} (N - 1) + \sum_{i = 1}^{N - 2} x_{i}^{2} + {[\sum_{i = 1}^{N - 2} x_{i}]}^{2} - 2 N \bar{x} \sum_{i = 1}^{N - 2} x_{i} - (N - 1) s^{2}

$2(N\bar{x}\! - \!\sum_{i=1}^{N-2}x_i)z - 2z^2 = N\bar{x}^2(N\!-\!1) + \sum_{i=1}^{N-2}x_i^2 + \left[\sum_{i=1}^{N-2}x_i\right]^2 - 2N\bar{x}\sum_{i=1}^{N-2}x_i - (N\!-\!1)s^2$ Si tomamos , , como negación del RHS, podemos resolver usando la fórmula cuadrática . Por ejemplo, en , se podría usar el siguiente código:

a = - 2

$a=-2$

b = 2 (N \bar{x} - \sum_{i = 1}^{N - 2} x_{i})

$b=2(N\bar{x} - \sum_{i=1}^{N-2}x_i)$

c

$c$

z

$z$ R

find.yz = function(x, xbar, s2){
  N    = length(x) + 2
  sumx = sum(x)
  sx2  = as.numeric(x%*%x)          # this is the sum of x^2
  a    = -2
  b    = 2*(N*xbar - sumx)
  c    = -N*xbar^2*(N-1) - sx2 - sumx^2 + 2*N*xbar*sumx + (N-1)*s2
  rt   = sqrt(b^2 - 4*a*c)

  z    = (-b + rt)/(2*a)
  y    = N*xbar - (sumx + z)
  newx = c(x, y, z)
  return(newx)
}

set.seed(62)
x    = rnorm(2)
newx = find.yz(x, xbar=0, s2=1)
newx                                # [1] 0.8012701  0.2844567  0.3757358 -1.4614627
mean(newx)                          # [1] 0
var(newx)                           # [1] 1

Hay algunas cosas que debes entender sobre este enfoque. Primero, no se garantiza que funcione. Por ejemplo, es posible que sus datos iniciales de sean tales que no existan valores y que hagan que la varianza del conjunto resultante sea igual a . Considerar: $N-2$ $y$ $z$ $s^2$

set.seed(22)    
x    = rnorm(2)
newx = find.yz(x, xbar=0, s2=1)
Warning message:
In sqrt(b^2 - 4 * a * c) : NaNs produced
newx                                # [1] -0.5121391  2.4851837        NaN        NaN
var(c(x, mean(x), mean(x)))         # [1] 1.497324

En segundo lugar, mientras que la estandarización hace que las distribuciones marginales de todas sus variantes sean más uniformes, este enfoque solo afecta los dos últimos valores, pero hace que sus distribuciones marginales estén sesgadas:

set.seed(82)
xScaled = matrix(NA, ncol=4, nrow=10000)
for(i in 1:10000){
  x           = rnorm(4)
  xScaled[i,] = scale(x)
}

(insertar gráfico)

set.seed(82)
xDf = matrix(NA, ncol=4, nrow=10000)
i   = 1
while(i<10001){
  x       = rnorm(2)
  xDf[i,] = try(find.yz(x, xbar=0, s2=2), silent=TRUE)  # keeps the code from crashing
  if(!is.nan(xDf[i,4])){ i = i+1 }                      # increments if worked
}

(insertar gráfico)

Tercero, la muestra resultante puede no parecer muy normal; puede parecer que tiene "valores atípicos" (es decir, puntos que provienen de un proceso de generación de datos diferente al resto), ya que ese es esencialmente el caso. Es menos probable que esto sea un problema con tamaños de muestra más grandes, ya que las estadísticas de muestra de los datos generados deben converger a los valores requeridos y, por lo tanto, necesitan menos ajustes. Con muestras más pequeñas, siempre puede combinar este enfoque con un algoritmo de aceptación / rechazo que lo intenta nuevamente si la muestra generada tiene estadísticas de forma (por ejemplo, asimetría y curtosis) que están fuera de los límites aceptables (cf. comentario de @ cardinal ), o ampliar este enfoque para generar una muestra con una media fija, varianza, asimetría ycurtosis (aunque te dejaré el álgebra) Alternativamente, podría generar una pequeña cantidad de muestras y utilizar la que tenga la estadística más pequeña (digamos) de Kolmogorov-Smirnov.

library(moments)
set.seed(7900)  
x = rnorm(18)
newx.ss7900 = find.yz(x, xbar=0, s2=1)
skewness(newx.ss7900)                       # [1] 1.832733
kurtosis(newx.ss7900) - 3                   # [1] 4.334414
ks.test(newx.ss7900, "pnorm")$statistic     # 0.1934226

set.seed(200)  
x = rnorm(18)
newx.ss200 = find.yz(x, xbar=0, s2=1)
skewness(newx.ss200)                        # [1] 0.137446
kurtosis(newx.ss200) - 3                    # [1] 0.1148834
ks.test(newx.ss200, "pnorm")$statistic      # 0.1326304 

set.seed(4700)  
x = rnorm(18)
newx.ss4700 = find.yz(x, xbar=0, s2=1)
skewness(newx.ss4700)                       # [1]  0.3258491
kurtosis(newx.ss4700) - 3                   # [1] -0.02997377
ks.test(newx.ss4700, "pnorm")$statistic     # 0.07707929S

(agregar trama)

gung - Restablece a Monica
fuente

10

La técnica general es el 'Método de rechazo', donde simplemente rechaza los resultados que no cumplen con sus restricciones. A menos que tenga algún tipo de orientación (como MCMC), ¡podría generar muchos casos (dependiendo de su escenario) que se rechazan!

Cuando esté buscando algo así como una desviación media y estándar y pueda crear una métrica de distancia de algún tipo para decir qué tan lejos está de su objetivo, puede usar la optimización para buscar las variables de entrada que le dan la salida deseada valores.

Como un ejemplo feo donde buscaremos un vector uniforme aleatorio con longitud 100 que tenga media = 0 y desviación estándar = 1.

# simplistic optimisation example
# I am looking for a mean of zero and a standard deviation of one
# but starting from a plain uniform(0,1) distribution :-)
# create a function to optimise
fun <- function(xvec, N=100) {
  xmin <- xvec[1]
  xmax <- xvec[2]
  x <- runif(N, xmin, xmax)
  xdist <- (mean(x) - 0)^2 + (sd(x) - 1)^2
  xdist
}
xr <- optim(c(0,1), fun)

# now lets test those results
X <- runif(100, xr$par[1], xr$par[2])
mean(X) # approx 0
sd(X)   # approx 1

Sean
fuente

77

Las restricciones que ocurren con probabilidad cero son difíciles de satisfacer. ;-) Para el ejemplo específico en cuestión, un cambio y una dilatación apropiados logran fácilmente los objetivos establecidos , aunque uno podría analizar un poco más a fondo para ver cómo la distribución de los datos se ve perturbada por dicha operación.

cardenal

Gracias. Ciertamente, sería fácil rechazar observaciones menores que el mínimo y mayores que el máximo. Y puedo ver cómo podría definirlo como un problema de optimización. Sería genial ver algunos ejemplos o tal vez tener algunas sugerencias sobre qué leer a continuación.

Jeromy Anglim

1

@cardinal - de acuerdo. ¡Uno debe mirar las distribuciones (es decir, un histograma) de los números simulados de entrada y salida, ya que a veces pueden parecer muy extraños!

Sean

9

¿Hay algún programa en R que haga esto?

El paquete Runuran R contiene muchos métodos para generar variaciones aleatorias. Utiliza bibliotecas C del proyecto UNU.RAN (generador universal de números aleatorios no uniformes). Mi propio conocimiento del campo de la generación de variables aleatorias es limitado, pero la viñeta de Runuran ofrece una buena visión general. A continuación se muestran los métodos disponibles en el paquete Runuran, tomados de la viñeta:

Distribuciones continuas:

Muestreo de rechazo adaptativo
Rechazo de densidad transformada inversa
Interpolación polinómica de CDF inversa
Método simple de relación de uniformes
Rechazo de densidad transformada

Distribuciones discretas:

Inversión discreta de rechazo automático
Método Alias-Urna
Método de tabla de guía para inversión discreta

Distribuciones multivariantes:

Algoritmo de golpe y carrera con el método de relación de uniformes
Método de relación de uniformes ingenuo multivariante

Ejemplo:

Para un ejemplo rápido, suponga que desea generar una distribución Normal limitada entre 0 y 100:

require("Runuran")

## Normal distribution bounded between 0 and 100
d1 <- urnorm(n = 1000, mean = 50, sd = 25, lb = 0, ub = 100)

summary(d1)
sd(d1)
hist(d1)

La urnorm()función es una conveniente función de envoltura. Creo que detrás de escena utiliza el método de interpolación polinómica de CDF inversa, pero no estoy seguro. Para algo más complejo, digamos, una distribución normal discreta limitada entre 0 y 100:

require("Runuran")

## Discrete normal distribution bounded between 0 and 100
# Create UNU.RAN discrete distribution object
discrete <- unuran.discr.new(pv = dnorm(0:100, mean = 50, sd = 25), lb = 0, ub = 100)

# Create UNU.RAN object using the Guide-Table Method for Discrete Inversion
unr <- unuran.new(distr = discrete, method = "dgt")

# Generate random variates from the UNU.RAN object
d2 <- ur(unr = unr, n = 1000)

summary(d2)
sd(d2)
head(d2)
hist(d2)

jthetzel
fuente

3

¡Parece que hay un paquete R que cumple con su requisito publicado ayer! simstudy Por Keith Goldfeld

Simula conjuntos de datos para explorar técnicas de modelado o comprender mejor los procesos de generación de datos. El usuario especifica un conjunto de relaciones entre covariables y genera datos basados en estas especificaciones. Los conjuntos de datos finales pueden representar datos de ensayos de control aleatorio, diseños de medidas repetidas (longitudinales) y ensayos aleatorios grupales. La falta puede generarse utilizando varios mecanismos (MCAR, MAR, NMAR).

Tyelcie
fuente

1

Ni en la viñeta ni en la página de inicio del programa se menciona la reunión exacta de restricciones. ¿Por qué cree que este paquete cumple con el requisito de extraer de distribuciones condicionales?

gg

2

Esta es una respuesta que llega tan tarde que presumiblemente no tiene sentido, pero siempre hay una solución MCMC para la pregunta. Es decir, para proyectar la densidad conjunta de la muestra en la variedad definida por las restricciones, por ejemplo El único problema es entonces simular valores sobre esa variedad, es decir, encontrar una parametrización de la dimensión correcta. Un artículo de 2015 de Bornn, Shephard y Solgi estudia este mismo problema (con una respuesta interesante, si no la última ).

\prod_{i = 1}^{n} f (x_{i})

$\prod_{i=1}^n f(x_i)$

\sum_{i = 1}^{n} x_{i} = μ_{0} \sum_{i = 1}^{n} x_{i}^{2} = σ_{0}^{2}

$\sum_{i=1}^n x_i=\mu_0\qquad\sum_{i=1}^n x_i^2=\sigma_0^2$

Xi'an
fuente

2

Esta respuesta considera otro enfoque para el caso en el que desea forzar a las variaciones a ubicarse en un rango específico y, además, dictar la media y / o la varianza.

Restrinja nuestra atención al intervalo de la unidad . una media ponderada para la generalidad, así que arregle algunos pesos con , o configure si desea una ponderación estándar. Suponga que las cantidades y representan la media deseada (ponderada) y la varianza (ponderada), respectivamente. El límite superior en es necesario porque esa es la varianza máxima posible en un intervalo de unidades. Estamos interesados en dibujar algunas variantes de con estas restricciones de momento. $[0,1]$ $w_k\in[0,1]$ $\sum_{k=1}^Nw_k=1$ $w_k=1/N$ $\mu\in(0,1)$ $0<\sigma^2<\mu(1-\mu)$ $\sigma^2$ $x_1,...,x_N$ $[0,1]$

Primero dibujamos algunas variantes de cualquier distribución, como . Esta distribución afectará la forma de la distribución final. Luego los restringimos al intervalo unitario usando una función logística: $y_1,...,y_N$ $N(0,1)$ $[0,1]$

x_{k} = \frac{1}{1 + e^{- (y_{k} v - h)}}

$x_k=\frac{1}{1+e^{-(y_k v-h)}}$

Sin embargo, antes de hacer eso, como se ve en la ecuación anterior, transformamos los 's con la traducción y la escala . Esto es análogo a la primera ecuación en la respuesta de @ Macro. El truco ahora es elegir y para que las variables transformadas tengan los momentos deseados. Es decir, necesitamos uno o ambos de los siguientes elementos: $y_k$ $h$ $v$ $h$ $v$ $x_1,...,x_N$

μ = \sum_{k = 1}^{N} \frac{w_{k}}{1 + e^{- (y_{k} v - h)}} σ^{2} = \sum_{k = 1}^{N} \frac{w_{k}}{(1 + e^{- (y_{k} v - h)})^{2}} - {(\sum_{k = 1}^{N} \frac{w_{k}}{1 + e^{- (y_{k} v - h)}})}^{2}

$\mu=\sum_{k=1}^N \frac{w_k}{1+e^{-(y_k v-h)}} \\ \sigma^2=\sum_{k=1}^N \frac{w_k}{(1+e^{-(y_k v-h)})^2} - \left( \sum_{k=1}^N \frac{w_k}{1+e^{-(y_k v-h)}} \right)^2$

Invertir estas ecuaciones para y analíticamente no es factible, pero hacerlo numéricamente es sencillo, especialmente porque las derivadas con respecto a y son fáciles de calcular; solo toma unas pocas iteraciones del método de Newton. $v$ $h$ $v$ $h$

Como primer ejemplo, digamos que solo nos importa restringir la media ponderada y no la varianza. Fix , , , . Luego, para las distribuciones subyacentes , y terminamos con los siguientes histogramas, respectivamente, y de modo que la media de las variables es exactamente (incluso para pequeñas ): $\mu=0.8$ $v=1$ $w_k=1/N$ $N=200000$ $N(0,1)$ $N(0,0.1)$ $\text{Unif}(0,1)$ $0.8$ $N$

A continuación, limitemos la media y la varianza. Tome , , y considere las tres desviaciones estándar deseadas . Usando la misma distribución subyacente , aquí están los histogramas para cada uno: $\mu=0.2$ $w_k=1/N$ $N=2000$ $\sigma=0.1,0.05,0.01$ $N(0,1)$

Tenga en cuenta que estos pueden parecer un poco beta distribuidos, pero no lo son.

Ian Hincks
fuente

1

En mi respuesta aquí , enumeré tres paquetes R para hacer esto:

abalter
fuente

Debe haber algún formato para un enlace a referencias. ¿Debería ser un comentario en su lugar?

abalter

¿Cómo simular datos que satisfagan restricciones específicas, como tener una media específica y una desviación estándar?

Preguntas

Respuestas: