Wie können Daten simuliert werden, die bestimmte Bedingungen erfüllen, z. B. einen bestimmten Mittelwert und eine bestimmte Standardabweichung?

56

Diese Frage ist durch meine Frage zur Metaanalyse motiviert . Ich stelle mir jedoch vor, dass dies auch in Lehrkontexten nützlich ist, in denen Sie ein Dataset erstellen möchten, das genau einem vorhandenen veröffentlichten Dataset entspricht.

Ich weiß, wie man zufällige Daten aus einer bestimmten Distribution generiert. Wenn ich also zum Beispiel über die Ergebnisse einer Studie lese, die Folgendes hatte:

  • ein Mittelwert von 102,
  • eine Standardabweichung von 5,2 und
  • eine Stichprobengröße von 72.

Ich könnte ähnliche Daten mit rnormin R erzeugen . Zum Beispiel

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

Natürlich wären der Mittelwert und die SD nicht exakt gleich 102 bzw. 5,2:

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

Im Allgemeinen interessiert mich, wie man Daten simuliert, die eine Reihe von Einschränkungen erfüllen. Im obigen Fall sind die Konstanten Stichprobengröße, Mittelwert und Standardabweichung. In anderen Fällen kann es zu zusätzlichen Einschränkungen kommen. Zum Beispiel,

  • Möglicherweise sind ein Minimum und ein Maximum in den Daten oder der zugrunde liegenden Variablen bekannt.
  • Es ist möglicherweise bekannt, dass die Variable nur ganzzahlige oder nur nicht negative Werte annimmt.
  • Die Daten können mehrere Variablen mit bekannten Wechselbeziehungen enthalten.

Fragen

  • Wie kann ich im Allgemeinen Daten simulieren, die genau einer Reihe von Einschränkungen entsprechen?
  • Gibt es Artikel darüber? Gibt es Programme in R, die dies tun?
  • Wie könnte und sollte ich beispielsweise eine Variable so simulieren, dass sie einen bestimmten Mittelwert und eine bestimmte SD hat?
Jeromy Anglim
quelle
1
Warum sollen sie genau den veröffentlichten Ergebnissen entsprechen? Sind diese Schätzungen des Populationsmittelwerts und der Standardabweichung nicht gegeben, wenn ihre Stichprobe von Daten vorliegt? Wer ist angesichts der Unsicherheit bei diesen Schätzungen der Ansicht, dass die von Ihnen gezeigte Stichprobe nicht mit ihren Beobachtungen übereinstimmt?
Gavin Simpson
4
Da diese Frage anscheinend Antworten enthält, die die Marke verfehlen (IMHO), möchte ich darauf hinweisen, dass die Antwort konzeptionell einfach ist: Gleichheitsbeschränkungen werden wie Randverteilungen behandelt, und Ungleichheitsbeschränkungen sind multivariate Analoga der Kürzung. Das Abschneiden ist relativ einfach zu handhaben (häufig mit Rückweisungsabtastung); das schwierigere Problem besteht darin, einen Weg zu finden, um diese Randverteilungen abzutasten. Dies bedeutet, dass entweder Marginals unter Berücksichtigung der Verteilung und der Einschränkung abgetastet werden oder dass die Marginalverteilung ermittelt und daraus abgetastet wird.
whuber
4
Übrigens ist die letzte Frage für ortsbezogene Distributionsfamilien trivial. ZB x<-rnorm(72);x<-5.2*(x-mean(x))/sd(x)+102macht den Trick.
Whuber
1
@whuber, wie Kardinal in einem Kommentar zu meiner Antwort (in dem dieser "Trick" erwähnt wird) und einem Kommentar zu einer anderen Antwort anspielt - diese Methode hält die Variablen im Allgemeinen nicht in derselben Verteilungsfamilie, da Sie sich teilen um die Standardabweichung der Probe.
Makro
5
@Macro Das ist ein guter Punkt, aber vielleicht ist die beste Antwort: "Natürlich werden sie nicht die gleiche Verteilung haben"! Die gewünschte Verteilung ist die Verteilung , die von den Einschränkungen abhängig ist . Im Allgemeinen stammt dies nicht aus derselben Familie wie die übergeordnete Verteilung. Beispielsweise wird jedes Element einer Stichprobe der Größe 4 mit Mittelwert 0 und SD 1, das aus einer Normalverteilung gezogen wird, eine nahezu einheitliche Wahrscheinlichkeit für [-1,5, 1,5] haben, da die Bedingungen die möglichen Werte nach oben und unten begrenzen.
whuber

Antworten:

26

Im Allgemeinen können Sie die Variable entsprechend verschieben und skalieren, damit der Mittelwert und die Varianz Ihrer Stichprobe genau einem vorgegebenen Wert entsprechen. Insbesondere dann , wenn ist eine Stichprobe, dann die neuen VariablenX1,X2,...,Xn

Zi=c1(XiX¯sX)+c2

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

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

erzeugt einen Datensatz , der auf das Intervall . B1,...,Bn(a,b)

Hinweis: Diese Arten der Verschiebung / Skalierung ändern im Allgemeinen die Verteilungsfamilie der Daten, auch wenn die Originaldaten aus einer Familie mit Ortsskalen stammen.

Im Rahmen der Normalverteilung können Sie mit der mvrnormFunktion in Rnormale (oder multivariate) Daten mit einem vorgegebenen Stichprobenmittelwert / einer vorgegebenen Kovarianz simulieren empirical=TRUE. Insbesondere simuliert diese Funktion Daten aus der bedingten Verteilung einer normalverteilten Variablen, vorausgesetzt, der Stichprobenmittelwert und die (Co) Varianz entsprechen einem vorgegebenen Wert . Beachten Sie, dass die resultierenden Randverteilungen nicht normal sind, wie @whuber in einem Kommentar zur Hauptfrage anmerkt.

Hier ist ein einfaches univariates Beispiel, bei dem der Stichprobenmittelwert (aus einer Stichprobe von ) auf 0 und die Stichprobenstandardabweichung auf 1 beschränkt ist. Wir können sehen, dass das erste Element einer gleichmäßigen Verteilung weitaus ähnlicher ist als einer Normalverteilung Verteilung: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")

                  Bildbeschreibung hier eingeben

Makro
quelle
1
Die werden nicht normal verteilt, obwohl dies ungefähr der Fall sein kann, wenn die Stichprobengröße groß ist. Der erste Kommentar zu @ Seans Antwort spielt darauf an. Zi
Kardinal
1
Nun, das ist eine ziemlich natürliche Sache, die man tun möchte ... und oft verursacht sie nicht allzu viel Ärger.
Kardinal
1
+1. Im Beispiel ist die Uniform übrigens die genaue Antwort. (Der offensichtliche Abfall an den Enden des Diagramms ist ein Artefakt dafür, wie R Histogramme zeichnet.)
whuber
1
@whuber, danke, dass du dieses Beispiel motiviert hast. Angesichts der Tatsache, dass sich die Grenzverteilungen ändern, sobald Sie den Stichprobenmittelwert / die Stichprobenvarianz bestimmen, scheint es die beste "Antwort" im Sinne der OP-Frage zu sein, nur Daten mit einem Populationsmittelwert / einer Populationsvarianz zu simulieren, die der als Stichprobe angegebenen entsprechen Mengen (wie vom OP selbst vorgeschlagen), nicht wahr? Auf diese Weise erhalten Sie Probenmengen, die den gewünschten "ähnlich" sind, und die Grenzverteilungen sind genau so, wie Sie sie haben wollten.
Makro
1
@whuber, Wenn deine Stichprobe normal ist, dann hat eine Verteilung, ja? Die fragliche "neue" Variable ist nur eine lineare Kombination von . Ti=(XiX¯)/stTi
Makro
22

In Bezug auf Ihre Anfrage für Papiere gibt es:

Dies ist nicht ganz das, wonach Sie suchen, könnte aber als Mahlgut für die Mühle dienen.


Es gibt eine andere Strategie, die anscheinend niemand erwähnt hat. Es ist möglich, (Pseudo) Zufallsdaten aus einer Menge der Größe zu erzeugen, so dass die gesamte Menge Bedingungen erfüllt , solange die verbleibenden Daten auf geeignete Werte festgelegt sind. Die erforderlichen Werte sollten mit einem System aus Gleichungen, Algebra und etwas Ellbogenfett lösbar sein . NkNkkk

Um beispielsweise einen Satz von Daten aus einer Normalverteilung mit einem bestimmten Stichprobenmittelwert ( ) und einer Varianz ( zu generieren , müssen Sie die Werte von zwei Punkten festlegen: und . Da der Stichprobenmittelwert ist: muss sein: Die Stichprobenvarianz ist: also (nach Ersetzen von das Obige , Folieren / Verteilen & Umordnen ... ) wir bekommen: 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
Wenn wir , ist und als Negation der RHS können wir mit der quadratischen Formel nach auflösen . In könnte beispielsweise der folgende Code verwendet werden: 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

Es gibt einige Dinge zu verstehen, über diesen Ansatz. Erstens ist es nicht garantiert zu arbeiten. Zum Beispiel ist es möglich, dass Ihre anfänglichen Daten so sind, dass keine Werte und existieren, die die Varianz der resultierenden Menge gleich . Erwägen: N2yzs2

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

Zweitens: Während durch die Standardisierung die Randverteilungen aller Ihrer Variablen einheitlicher werden, wirkt sich dieser Ansatz nur auf die letzten beiden Werte aus, führt jedoch zu einer Verschiebung der Randverteilungen:

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

(Handlung einfügen)

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
}

(Handlung einfügen)

Drittens sieht die resultierende Stichprobe möglicherweise nicht ganz normal aus. Es könnte so aussehen, als hätte es Ausreißer (dh Punkte, die aus einem anderen Datenerzeugungsprozess stammen als der Rest), da dies im Wesentlichen der Fall ist. Dies ist bei größeren Stichproben weniger wahrscheinlich, da die Stichprobenstatistik aus den generierten Daten auf die erforderlichen Werte konvergieren sollte und daher weniger Anpassungen erforderlich sind. Bei kleineren Stichproben können Sie diesen Ansatz immer mit einem Annahme- / Ablehnungsalgorithmus kombinieren, der es erneut versucht, wenn die generierte Stichprobe Formstatistiken (z. B. Schiefe und Kurtosis) aufweist, die außerhalb akzeptabler Grenzen liegen (vgl. @ Kardinals Kommentar ) oder erweitert werden dieser Ansatz zur Erzeugung einer Stichprobe mit einem festen Mittelwert, einer festen Varianz, einer festen Schiefe undKurtosis (die Algebra überlasse ich Ihnen). Alternativ können Sie eine kleine Anzahl von Stichproben generieren und die mit der kleinsten Kolmogorov-Smirnov-Statistik verwenden.

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

(Grundstück hinzufügen)

gung - Wiedereinsetzung von Monica
quelle
10

Die allgemeine Technik ist die "Ablehnungsmethode", bei der Sie nur Ergebnisse ablehnen, die nicht Ihren Vorgaben entsprechen. Wenn Sie keine Anleitung (wie MCMC) haben, können Sie (abhängig von Ihrem Szenario) viele Fälle generieren, die abgelehnt werden!

Wenn Sie nach einem Mittelwert und einer Standardabweichung suchen und eine Entfernungsmetrik erstellen können, die angibt, wie weit Sie von Ihrem Ziel entfernt sind, können Sie mithilfe der Optimierung nach den Eingabevariablen suchen, die Ihnen die gewünschte Ausgabe liefern Werte.

Als ein hässliches Beispiel, in dem wir nach einem zufälligen gleichförmigen Vektor mit einer Länge von 100 suchen, der einen Mittelwert von 0 und eine Standardabweichung von 1 hat.

# 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
quelle
7
Einschränkungen, die mit der Wahrscheinlichkeit Null auftreten, sind schwer zu erfüllen. ;-) Für das konkrete Beispiel werden die festgelegten Ziele durch eine geeignete Verschiebung und Erweiterung leicht erreicht , obwohl man möglicherweise eine eingehendere Analyse vornehmen möchte, um zu untersuchen, wie die Verteilung der Daten durch eine solche Operation gestört wird.
Kardinal
Vielen Dank. Sicherlich wäre es einfach, Beobachtungen abzulehnen, die kleiner als die min und größer als die max sind. Und ich kann sehen, wie Sie es als Optimierungsproblem definieren können. Es wäre toll, einige Beispiele zu sehen oder Vorschläge zu haben, was als nächstes zu lesen ist.
Jeromy Anglim
1
@ Kardinal - einverstanden. Man sollte sich die Verteilungen (dh ein Histogramm) sowohl der eingegebenen simulierten Zahlen als auch der ausgegebenen Zahlen ansehen, da diese manchmal tatsächlich sehr seltsam aussehen können!
Sean
9

Gibt es Programme in R, die dies tun?

Das Runuran R-Paket enthält viele Methoden zum Generieren von Zufallsvariablen. Es verwendet C-Bibliotheken aus dem UNU.RAN-Projekt (Universal Non-Uniform RAndom Number Generator). Mein eigenes Wissen über das Gebiet der Zufallsgenerierung ist begrenzt, aber die Runuran- Vignette bietet einen guten Überblick. Im Folgenden sind die verfügbaren Methoden des Runuran-Pakets aufgeführt, die der Vignette entnommen wurden:

Kontinuierliche Verteilungen:

  • Adaptive Rückweisungsabtastung
  • Inverse Ablehnung der transformierten Dichte
  • Polynominterpolation von Inverse CDF
  • Einfaches Verhältnis der Uniformen
  • Ablehnung der transformierten Dichte

Diskrete Verteilungen:

  • Diskrete automatische Zurückweisungsumkehr
  • Alias-Urn-Methode
  • Guide-Table-Methode für diskrete Inversion

Multivariate Verteilungen:

  • Hit-and-Run-Algorithmus mit der Ratio-of-Uniforms-Methode
  • Multivariate naive Verhältnis-von-Uniformen-Methode

Beispiel:

Angenommen, Sie möchten eine Normalverteilung zwischen 0 und 100 generieren:

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)

Die urnorm()Funktion ist eine praktische Wrapper-Funktion. Ich glaube, dass es hinter den Kulissen die Polynomial Interpolation of Inverse CDF-Methode verwendet, bin mir aber nicht sicher. Für etwas komplexeres, sagen wir eine diskrete Normalverteilung, die zwischen 0 und 100 liegt:

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
quelle
3

Es scheint, dass es ein R-Paket gibt, das Ihren Anforderungen entspricht, das erst gestern veröffentlicht wurde! simstudy Von Keith Goldfeld

Simuliert Datensätze, um Modellierungstechniken zu untersuchen oder Datenerzeugungsprozesse besser zu verstehen. Der Benutzer gibt eine Reihe von Beziehungen zwischen Kovariaten an und generiert Daten basierend auf diesen Spezifikationen. Die endgültigen Datensätze können Daten aus randomisierten Kontrollversuchen, (Längsschnitt-) Versuchsplänen mit wiederholten Messungen und randomisierten Clusterversuchen darstellen. Fehlzeiten können mit verschiedenen Mechanismen (MCAR, MAR, NMAR) erzeugt werden.

Tyelcie
quelle
1
Weder in der Vignette noch auf der Homepage des Programms wird die genaue Einhaltung der Auflagen erwähnt. Warum erfüllt dieses Paket Ihrer Meinung nach die Anforderung, aus bedingten Distributionen zu ziehen?
gg
2

Dies ist eine Antwort, die so spät kommt, dass sie vermutlich bedeutungslos ist, aber es gibt immer eine MCMC-Lösung für die Frage. Nämlich, um die Verbindungsdichte der Probe auf den durch die Bedingungen definierten Verteiler zu projizieren , zum Beispiel Das einzige Problem ist dann, Werte über diesen Verteiler zu simulieren, dh eine Parametrisierung der richtigen Dimension zu finden. Eine Veröffentlichung von Bornn, Shephard und Solgi aus dem Jahr 2015 untersucht genau dieses Problem (mit einer interessanten, wenn nicht endgültigen Antwort ).

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

Diese Antwort berücksichtigt einen anderen Ansatz für den Fall, dass Sie die Variationen zwingen möchten, in einem bestimmten Bereich zu liegen, und zusätzlich den Mittelwert und / oder die Varianz vorgeben möchten .

Beschränken Sie unsere Aufmerksamkeit auf das Einheitsintervall . Lassen Sie uns einen gewichteten Mittelwert für die Allgemeinheit verwenden, also setzen Sie einige Gewichte mit , oder setzen Sie wenn Sie eine Standardgewichtung wünschen. Angenommen, die Größen und repräsentieren den gewünschten (gewichteten) Mittelwert bzw. die gewünschte (gewichtete) Varianz. Die Obergrenze für ist erforderlich, da dies die maximal mögliche Varianz für ein Einheitsintervall ist. Wir sind daran interessiert, einige Variablen aus mit diesen Momenteinschränkungen zu zeichnen .[0,1]wk[0,1]k=1Nwk=1wk=1/Nμ(0,1)0<σ2<μ(1μ)σ2x1,...,xN[0,1]

Zuerst zeichnen wir einige Variablen aus einer beliebigen Verteilung wie . Diese Verteilung beeinflusst die Form der endgültigen Verteilung. Dann beschränken wir sie mit einer logistischen Funktion auf das Einheitsintervall :y1,...,yNN(0,1)[0,1]

xk=11+e(ykvh)

Bevor wir dies tun, transformieren wir jedoch, wie in der obigen Gleichung gezeigt, die mit der Translation und der Skala . Dies ist analog zu der ersten Gleichung in der Antwort von @ Macro. Der Trick besteht nun darin, und so zu wählen , dass die transformierten Variablen die gewünschten Momente haben. Das heißt, wir benötigen eine oder beide der folgenden : ykhvhvx1,...,xN

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

Das analytische Invertieren dieser Gleichungen für und ist nicht möglich, aber numerisch ist unkompliziert, zumal Ableitungen in Bezug auf und leicht zu berechnen sind. es dauert nur ein paar Iterationen von Newtons Methode.vhvh

Als erstes Beispiel wollen wir nur den gewichteten Mittelwert und nicht die Varianz einschränken. Fix , , , . Dann erhalten wir für die zugrunde liegenden Verteilungen , und die folgenden Histogramme, und zwar so, dass der Mittelwert der Variablen genau beträgt (auch für kleine ):v = 1 w k = 1 / N N = 200000 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

Beispiel 1

Als nächstes beschränken wir sowohl den Mittelwert als auch die Varianz. Nehmen Sie , , und betrachten Sie die drei gewünschten Standardabweichungen . Unter Verwendung der gleichen zugrunde liegenden Verteilung sind hier die Histogramme für jedes: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)

Beispiel 2

Beachten Sie, dass diese möglicherweise ein bisschen Beta-verteilt aussehen, aber nicht.

Ian Hincks
quelle
1

In meiner Antwort hier habe ich drei R-Pakete dafür aufgelistet:

abalter
quelle
Es muss ein Format für einen Link zu Referenzen geben. Sollte es stattdessen ein Kommentar sein?
abalter