Comment simuler des données qui satisfont à des contraintes spécifiques, telles que la moyenne et l’écart-type spécifiques?

56

Cette question est motivée par ma question sur la méta-analyse . Mais j'imagine que cela serait également utile dans les contextes pédagogiques dans lesquels vous souhaitez créer un jeu de données qui reflète exactement un jeu de données publié existant.

Je sais comment générer des données aléatoires à partir d'une distribution donnée. Ainsi, par exemple, si je lisais les résultats d’une étude portant sur:

  • une moyenne de 102,
  • un écart type de 5,2, et
  • une taille d'échantillon de 72.

Je pourrais générer des données similaires en utilisant rnormdans R. Par exemple,

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

Bien sûr, la moyenne et l'écart-type ne seraient pas exactement égaux à 102 et 5,2 respectivement:

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

En général, je suis intéressé par la façon de simuler des données qui satisfont à un ensemble de contraintes. Dans le cas ci-dessus, les constantes sont la taille de l'échantillon, la moyenne et l'écart type. Dans d'autres cas, il peut y avoir des contraintes supplémentaires. Par exemple,

  • un minimum et un maximum dans les données ou dans la variable sous-jacente peuvent être connus.
  • la variable peut être connue pour ne prendre que des valeurs entières ou uniquement des valeurs non négatives.
  • les données peuvent inclure plusieurs variables avec des corrélations connues.

Des questions

  • En général, comment puis-je simuler des données qui répondent exactement à un ensemble de contraintes?
  • Y a-t-il des articles écrits à ce sujet? Y at-il des programmes dans R qui font cela?
  • Par exemple, comment puis-je et devrais-je simuler une variable de sorte qu'elle ait une moyenne et une dd spécifiques?
Jeromy Anglim
la source
1
Pourquoi voulez-vous qu'ils soient exactement comme les résultats publiés? Ces estimations de la moyenne et de l'écart type de la population ne sont-elles pas basées sur leur échantillon de données? Étant donné l’incertitude de ces estimations, qui peut dire que l’échantillon que vous montrez ci-dessus ne correspond pas à leurs observations?
Réintégrer Monica - G. Simpson le
4
Étant donné que cette question semble collecter des réponses qui manquent la note (IMHO), je voudrais souligner que, sur le plan conceptuel, la réponse est simple: les contraintes d' égalité sont traitées comme des distributions marginales et les contraintes d' inégalité sont des analogues multivariés de troncature. La troncature est relativement facile à gérer (souvent avec un échantillonnage de rejet); le problème le plus difficile consiste à trouver un moyen d'échantillonner ces distributions marginales. Cela signifie soit des marges d'échantillonnage compte tenu de la distribution et de la contrainte, soit une intégration pour trouver la distribution marginale et en tirer un échantillon.
whuber
4
BTW, la dernière question est triviale pour les familles de distribution échelle de localisation. Par exemple, x<-rnorm(72);x<-5.2*(x-mean(x))/sd(x)+102fait le tour.
whuber
1
@whuber, comme le dit cardinal dans un commentaire à ma réponse (qui mentionne cette "astuce") et un commentaire dans une autre réponse - cette méthode, en général, ne gardera pas les variables dans la même famille de distribution, car vous divisez par l’écart type de l’échantillon.
Macro
5
@ Macro C'est un bon point, mais peut-être que la meilleure réponse est "bien sûr, ils n'auront pas la même distribution"! La distribution souhaitée est la distribution conditionnelle aux contraintes. En général, cela ne sera pas de la même famille que la distribution parente. Par exemple, chaque élément d'un échantillon de taille 4 avec une moyenne de 0 et SD 1 tirés d'une distribution normale aura une probabilité presque uniforme sur [-1.5, 1.5], car les conditions placent des limites supérieure et inférieure sur les valeurs possibles.
whuber

Réponses:

26

En général, pour que la moyenne et la variance de votre échantillon soient exactement égales à une valeur prédéfinie, vous pouvez déplacer et mettre à l'échelle la variable de manière appropriée. Plus précisément, si est un échantillon, alors les nouvelles variablesX1,X2,...,Xn

Zi=c1(XiX¯sX)+c2

X¯=1ni=1nXisX2=1n1i=1n(XiX¯)2Zic2c1

Bi=a+(ba)(Ximin({X1,...,Xn})max({X1,...,Xn})min({X1,...,Xn}))

produira un ensemble de données limité à l'intervalle . B1,...,Bn(a,b)

Remarque: ces types de décalage / mise à l'échelle modifieront généralement la famille de distribution des données, même si les données d'origine proviennent d'une famille d'échelle de localisation.

Dans le contexte de la distribution normale, la mvrnormfonction R vous permet de simuler des données normales (ou normales à plusieurs variables) avec une moyenne / covariance d' échantillon prédéfinie par réglage empirical=TRUE. Spécifiquement, cette fonction simule les données de la distribution conditionnelle d'une variable normalement distribuée, étant donné que la moyenne de l'échantillon et la (co) variance sont égales à une valeur prédéfinie . Notez que les distributions marginales résultantes ne sont pas normales, comme l'a souligné @whuber dans un commentaire à la question principale.

Voici un exemple univarié simple où la moyenne de l'échantillon (d'un échantillon de ) est contrainte à 0 et l'écart type de l'échantillon est 1. Nous pouvons voir que le premier élément est beaucoup plus similaire à une distribution uniforme qu'à une normale Distribution: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")

                  entrez la description de l'image ici

Macro
la source
1
Le ne sera pas distribué normalement, bien qu'il puisse l'être approximativement si la taille de l'échantillon est grande. Le premier commentaire à la réponse de @ Sean y fait allusion. Zi
cardinal
1
Eh bien, c'est une chose assez naturelle que vous voulez faire ... et qui ne cause souvent pas trop de problèmes.
cardinal
1
+1 Dans l'exemple, l'uniforme est la réponse exacte , à propos. (La chute apparente aux extrémités de l'intrigue est un artefact de la façon dont R dessine les histogrammes.)
whuber
1
@ Whuber, merci d'avoir motivé cet exemple. Étant donné que les distributions marginales changent une fois que vous conditionnez la moyenne / variance de l'échantillon, il semble que la meilleure "réponse" dans l'esprit de la question du PO consiste simplement à simuler des données avec une moyenne / variance de la population égale à celle indiquée comme échantillon. quantités (comme suggéré par le PO lui-même), n'est-ce pas? De cette façon, vous obtenez des quantités d'échantillons "similaires" à celles souhaitées, et les distributions marginales sont ce que vous souhaitiez.
Macro
1
@whuber, si votre échantillon est normal, alors a une distribution , oui? La "nouvelle" variable en question sera simplement une combinaison linéaire de . Ti=(XiX¯)/stTi
Macro
22

En ce qui concerne votre demande de papiers, il y a:

Ce n’est pas tout à fait ce que vous cherchez, mais pourrait servir à abreuver le moulin.


Il y a une autre stratégie que personne ne semble avoir mentionnée. Il est possible de générer (pseudo) données aléatoires à partir d'un ensemble de taille sorte qu'un ensemble complet respecte contraintes, pour autant que les données restantes soient fixées à des valeurs appropriées. Les valeurs requises doivent pouvoir être résolues avec un système d' équations , d'algèbre et de graisse pour coudes. NkNkkk

Par exemple, pour générer un ensemble de données à partir d'une distribution normale qui aura une moyenne d'échantillon donnée, , et une variance, , vous devrez fixer les valeurs de deux points: et . Étant donné que la moyenne de l'échantillon est: doit être: L'écart exemple est: ainsi (après avoir remplacé ce qui précède par , déjouant / distribuant, et réarrangeant ... ) on a: Nx¯s2yz

x¯=i=1N2xi+y+zN
y
y=Nx¯(i=1N2xi+z)
s2=i=1N2(xix¯)2+(yx¯)2+(zx¯)2N1
y
2(Nx¯i=1N2xi)z2z2=Nx¯2(N1)+i=1N2xi2+[i=1N2xi]22Nx¯i=1N2xi(N1)s2
Si nous prenons , , et comme négation de RHS, nous pouvons résoudre pour utilisant la formule quadratique . Par exemple, dans , le code suivant pourrait être utilisé: a=2b=2(Nx¯i=1N2xi)czR
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

Il y a certaines choses à comprendre à propos de cette approche. Tout d'abord, il n'est pas garanti de fonctionner. Par exemple, il est possible que vos données initiales soient telles qu'il n'existe pas de valeurs et rendant la variance de l'ensemble résultant égale à . Considérer: y z s 2N2yzs2

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

Deuxièmement, alors que la standardisation uniformise les distributions marginales de toutes vos variables, cette approche n’affecte que les deux dernières valeurs, mais rend leurs distributions marginales faussées:

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

(insérer la parcelle)

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
}

(insérer la parcelle)

Troisièmement, l'échantillon résultant peut ne pas sembler très normal; il peut sembler qu'il a des «valeurs aberrantes» (c'est-à-dire des points qui proviennent d'un processus de génération de données différent du reste), car c'est essentiellement le cas. Ceci est moins susceptible de poser problème avec des tailles d'échantillon plus grandes, car les statistiques d'échantillon à partir des données générées doivent converger vers les valeurs requises et nécessitent par conséquent moins d'ajustement. Avec des échantillons plus petits, vous pouvez toujours combiner cette approche avec un algorithme d'acceptation / rejet qui tente à nouveau si l'échantillon généré a des statistiques de forme (par exemple, asymétrie et kurtosis) qui sont en dehors des limites acceptables (cf., commentaire de @ cardinal ), ou étendent cette approche pour générer un échantillon avec une moyenne, une variance, une asymétrie etkurtosis (je laisserai cependant l'algèbre à vous). Alternativement, vous pouvez générer un petit nombre d’échantillons et utiliser celui qui contient la statistique la plus petite (disons) 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

(ajouter un complot)

gung - Rétablir Monica
la source
10

La technique générale est la «méthode de rejet», dans laquelle vous rejetez simplement des résultats qui ne respectent pas vos contraintes. À moins que vous n'ayez une sorte de guide (comme MCMC), vous pourriez générer beaucoup de cas (selon votre scénario) qui sont rejetés!

Lorsque vous recherchez quelque chose comme un écart moyen et standard et que vous pouvez créer une métrique de distance pour indiquer votre distance par rapport à votre objectif, vous pouvez utiliser l'optimisation pour rechercher les variables d'entrée qui vous donnent le résultat souhaité. valeurs.

Comme un vilain exemple où nous allons chercher un vecteur aléatoire uniforme de longueur 100 qui a une moyenne = 0 et l' écart type = 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
la source
7
Les contraintes qui surviennent avec une probabilité nulle sont difficiles à satisfaire. ;-) En ce qui concerne l’exemple spécifique, un décalage et une dilatation appropriés permettent d’atteindre facilement les objectifs énoncés , bien que l’on veuille analyser un peu plus en profondeur pour examiner comment la distribution des données est perturbée par une telle opération.
cardinal
Merci. Certes, il serait facile de rejeter les observations inférieures au minimum et supérieures au maximum. Et je peux voir comment vous pourriez le définir comme un problème d'optimisation. Il serait bon de voir quelques exemples ou peut-être quelques suggestions de choses à lire par la suite.
Jeromy Anglim
1
@ cardinal - d'accord. Il faut regarder les distributions (c.-à-d. Un histogramme) des nombres simulés et des résultats en entrée, car ils peuvent parfois sembler très étranges!
Sean
9

Y at-il des programmes dans R qui font cela?

Le paquet Runuran R contient de nombreuses méthodes pour générer des variables aléatoires. Il utilise les bibliothèques C du projet UNU.RAN ( générateur universel de nombres non uniformes de RAndom). Ma connaissance du domaine de la génération de variables aléatoires est limitée, mais la vignette Runuran fournit un bon aperçu. Voici les méthodes disponibles dans le package Runuran, tiré de la vignette:

Distributions continues:

  • Échantillonnage de rejet adaptatif
  • Rejet de densité transformée inverse
  • Interpolation polynomiale des CDF inverses
  • Méthode simple du ratio d'uniformes
  • Rejet de densité transformée

Distributions discrètes:

  • Inversion de réjection automatique discrète
  • Méthode Alias-Urn
  • Méthode table-guide pour l'inversion discrète

Distributions multivariées:

  • Algorithme Hit-and-Run avec la méthode Ratio-of-Uniforms
  • Méthode du ratio naif multivarié des uniformes

Exemple:

Pour un exemple rapide, supposons que vous vouliez générer une distribution normale limitée entre 0 et 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()fonction est une fonction d'encapsulation pratique. Je crois que dans les coulisses, il utilise la méthode d’interpolation polynomiale de CDF inverse, mais je ne suis pas sûr. Pour quelque chose de plus complexe, disons une distribution normale discrète limitée entre 0 et 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
la source
3

Il semble qu’un package R répondant à vos besoins a été publié hier! simstudy Par Keith Goldfeld

Simule des ensembles de données afin d'explorer des techniques de modélisation ou de mieux comprendre les processus de génération de données. L'utilisateur spécifie un ensemble de relations entre les covariables et génère des données en fonction de ces spécifications. Les ensembles de données finaux peuvent représenter des données provenant d'essais contrôlés randomisés, de conceptions de mesures répétées (longitudinales) et d'essais randomisés en grappes. Les manquements peuvent être générés à l'aide de divers mécanismes (MCAR, MAR, NMAR).

Tyelcie
la source
1
Ni dans la vignette, ni sur la page d'accueil du programme, on ne mentionne le respect exact des contraintes. Pourquoi pensez-vous que ce paquet satisfait à l'exigence de tirer des distributions conditionnelles?
gg le
2

C’est une réponse qui arrive si tard qu’elle n’a apparemment aucun sens, mais il existe toujours une solution MCMC à la question. À savoir, pour projeter la densité de joint de l'échantillon sur la variété définie par les contraintes, par exemple Le seul problème est alors de simuler des valeurs sur cette variété, c'est-à-dire de trouver un paramétrage de la dimension correcte. Un article publié en 2015 par Bornn, Shephard et Solgi étudie ce problème même (avec une réponse intéressante, voire ultime ).

i=1nf(xi)
i=1nxi=μ0i=1nxi2=σ02
Xi'an
la source
2

Cette réponse considère une autre approche du cas où vous voulez forcer les variables à se situer dans une plage spécifiée et dicter en plus la moyenne et / ou la variance.

Limitez notre attention à l'intervalle unitaire . Utilisons une moyenne pondérée pour la généralité, donc certaines pondérations avec ou définissez si vous souhaitez une pondération standard. Supposons que les quantités et représentent respectivement la variance souhaitée (pondérée) et la variance (pondérée). La limite supérieure de est nécessaire car c'est la variance maximale possible sur un intervalle d'unité. Nous sommes intéressés par dessiner des variables de avec ces restrictions de moment.[0,1]wk[0,1]k=1Nwk=1wk=1/Nμ(0,1)0<σ2<μ(1μ)σ2x1,...,xN[0,1]

Commençons par dessiner des variables de toute distribution, comme . Cette distribution affectera la forme de la distribution finale. Ensuite, nous les contraignons à l'intervalle unitaire aide d'une fonction logistique:y1,...,yNN(0,1)[0,1]

xk=11+e(ykvh)

Avant de faire cela, cependant, comme le montre l'équation ci-dessus, nous transformons les avec la traduction et l'échelle . Ceci est analogue à la première équation de la réponse de @ Macro. L'astuce consiste maintenant à choisir et pour que les variables transformées aient le ou les moments souhaités. C’est-à-dire qu’il faut l’un ou l’autre des éléments suivants: ykhvhvx1,...,xN

μ=k=1Nwk1+e(ykvh)σ2=k=1Nwk(1+e(ykvh))2(k=1Nwk1+e(ykvh))2

L'inversion analytique de ces équations pour et n'est pas réalisable, mais le faire numériquement est simple, d'autant plus que les dérivées par rapport à et sont faciles à calculer; il ne faut que quelques itérations de la méthode de Newton.vhvh

Comme premier exemple, disons que nous ne nous soucions que de limiter la moyenne pondérée et non la variance. Fixer , , , . Ensuite, pour les distributions sous-jacentes , et nous obtenons respectivement les histogrammes suivants, de sorte que la moyenne des variables est exactement égale à (même pour petit ):v = 1 w k = 1 / N N = 200 000 N ( 0 , 1 ) N ( 0 , 0,1 ) Unif ( 0 , 1 ) 0,8 Nμ=0.8v=1wk=1/NN=200000N(0,1)N(0,0.1)Unif(0,1) 0.8N

Exemple 1

Ensuite, limitons la moyenne et la variance. Prenez , , et considérez les trois écarts types souhaités . En utilisant la même distribution sous-jacente , voici les histogrammes pour chacun:w k = 1 / N N = 2000 σ = 0,1 , 0,05 , 0,01 N ( 0 , 1 )μ=0.2wk=1/NN=2000σ=0.1,0.05,0.01N(0,1)

Exemple 2

Notez que ceux-ci peuvent sembler un peu bêta-distribués, mais ils ne le sont pas.

Ian Hincks
la source
1

Dans ma réponse ici , j'ai énuméré trois paquets R pour faire ceci:

abalter
la source
Il doit y avoir un format pour un lien vers des références. Cela devrait-il être un commentaire?
Abalter