Jak symulować dane spełniające określone ograniczenia, takie jak posiadanie określonej średniej i odchylenia standardowego?

To pytanie jest motywowane moim pytaniem dotyczącym metaanalizy . Ale wyobrażam sobie, że przydałoby się to również w nauczaniu kontekstów, w których chcesz utworzyć zestaw danych, który dokładnie odzwierciedla istniejący opublikowany zestaw danych.

Wiem, jak generować losowe dane z danej dystrybucji. Na przykład, jeśli przeczytam o wynikach badania, które:

• średnio 102,
• odchylenie standardowe 5,2, oraz
• wielkość próbki 72.

Mógłbym wygenerować podobne dane używając rnormR. Na przykład

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

Oczywiście średnia i SD nie byłyby dokładnie równe odpowiednio 102 i 5,2:

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

Ogólnie interesuje mnie, jak symulować dane, które spełniają zestaw ograniczeń. W powyższym przypadku stałymi są wielkość próby, średnia i odchylenie standardowe. W innych przypadkach mogą istnieć dodatkowe ograniczenia. Na przykład,

• minimalne i maksymalne dane lub zmienna bazowa mogą być znane.
• wiadomo, że zmienna przyjmuje tylko wartości całkowite lub tylko wartości nieujemne.
• dane mogą obejmować wiele zmiennych o znanych wzajemnych korelacjach.

pytania

• Ogólnie, jak mogę symulować dane, które dokładnie spełniają zestaw ograniczeń?
• Czy są na ten temat artykuły? Czy są jakieś programy w R, które to robią?
• Na przykład, w jaki sposób i powinienem symulować zmienną, aby miała określoną średnią i sd?

Ogólnie rzecz biorąc, aby średnia i wariancja w próbce były dokładnie równe z góry określonej wartości, można odpowiednio przesunąć i skalować zmienną. W szczególności, jeśli jest próbką, a następnie nowymi zmiennymi$X_1, X_2, ..., X_n$

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

gdzie $\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) \left( \frac{ X_i - \min (\{X_1, ..., X_n\}) }{\max (\{X_1, ..., X_n\}) - \min (\{X_1, ..., X_n\}) } \right)$$

utworzy zestaw danych który jest ograniczony do przedziału . $B_1, ..., B_n$$(a,b)$

Uwaga: Te typy przesunięcia / skalowania zasadniczo zmienią rodzinę dystrybucyjną danych, nawet jeśli oryginalne dane pochodzą z rodziny o skali lokalizacji.

W kontekście tego rozkładu normalnegomvrnorm funkcja w R pozwala symulować normalne (lub wielowymiarowych normalny) dane z góry określonej próbki oznaczać / kowariancji przez ustawienie empirical=TRUE. W szczególności funkcja ta symuluje dane z rozkładu warunkowego zmiennej normalnie rozłożonej, biorąc pod uwagę średnią próbkę i (ko) wariancję równą z góry określonej wartości . Zauważ, że wynikowe rozkłady krańcowe nie są normalne, jak zauważył @whuber w komentarzu do głównego pytania.

Oto prosty przykład z jedną zmienną, w którym średnia próbki (z próbki ) jest ograniczona do 0, a odchylenie standardowe próbki wynosi 1. Widzimy, że pierwszy element jest znacznie bardziej podobny do rozkładu równomiernego niż normalny dystrybucja:$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")

$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $

Jeśli chodzi o twoją prośbę o dokumenty, istnieją:

• Chatterjee, S. & Firat, A. (2007). Generowanie danych z identycznymi statystykami, ale z odmienną grafiką: kontynuacja zbioru danych Anscombe. The American Statistician, 61 , 3, s. 248–254.

Nie jest to dokładnie to, czego szukasz, ale może służyć jako młyn do młyna.

---

Jest inna strategia, o której nikt chyba nie wspomniał. Możliwe jest generowanie (pseudo) losowych danych ze zbioru o rozmiarze tak że cały zestaw spełnia ograniczenia, o ile pozostałe danych jest ustawione na odpowiednie wartości. Wymagane wartości powinny być możliwe do rozwiązania za pomocą układu równań , algebry i pewnego smaru łokciowego. $N-k$$N$$k$$k$$k$

Na przykład, aby wygenerować zestaw danych z rozkładu normalnego, który będzie miał podaną średnią próbkową, i wariancję, , musisz naprawić wartości dwóch punktów: i . Ponieważ średnia próbki to: musi być: Przykładowa wariancja to: ten sposób (po zamianie powyższego na , foliowanie / dystrybucja i zmiana kolejności ... ) otrzymujemy: $N$$\bar x$$s^2$$y$$z$

$$\bar x = \frac{\sum_{i=1}^{N-2}x_i\; + \;y\!+\!z}{N}$$ $y$ $$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}$$ $y$ $$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$$ Jeśli weźmiemy , , oraz jako negacja RHS, możemy rozwiązać dla za pomocą wzoru kwadratowego . Na przykład można użyć następującego kodu: $a=-2$$b=2(N\bar{x} - \sum_{i=1}^{N-2}x_i)$$c$$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

Podejście to należy zrozumieć. Po pierwsze, nie ma gwarancji, że zadziała. Na przykład, możliwe jest, że początkowe Dane są takie, że nie ma wartości i tego, że istnieją będzie wariancji Otrzymaną równy . Rozważać: $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

Po drugie, podczas gdy standaryzacja sprawia, że ​​rozkłady krańcowe wszystkich twoich zmiennych są bardziej jednolite, to podejście wpływa tylko na dwie ostatnie wartości, ale sprawia, że ​​ich rozkłady krańcowe są wypaczone:

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

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
}

Po trzecie, uzyskana próbka może nie wyglądać bardzo normalnie; może to wyglądać tak, jakby zawierało „wartości odstające” (tj. punkty, które pochodzą z innego procesu generowania danych niż reszta), ponieważ tak jest w istocie. Jest to mniej prawdopodobne, że będzie to stanowić problem w przypadku większych próbek, ponieważ statystyki próbek z wygenerowanych danych powinny być zbieżne z wymaganymi wartościami, a zatem wymagać mniejszego dostosowania. Przy mniejszych próbkach zawsze możesz połączyć to podejście z algorytmem akceptowania / odrzucania, który próbuje ponownie, jeśli wygenerowana próbka ma statystyki kształtu (np. Skośność i kurtoza), które są poza dopuszczalnymi granicami (por. Komentarz kardynała ) lub rozszerzyć takie podejście do generowania próbki o ustalonej średniej, wariancji, skośności ikurtoza (ale zostawię tobie algebrę). Alternatywnie, możesz wygenerować niewielką liczbę próbek i użyć tej z najmniejszą (powiedzmy) statystyką Kołmogorowa-Smirnowa.

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

Ogólna technika to „Metoda odrzucenia”, w której po prostu odrzucasz wyniki, które nie spełniają twoich ograniczeń. Jeśli nie masz wskazówek (np. MCMC), możesz generować wiele przypadków (w zależności od scenariusza), które są odrzucane!

Jeśli szukasz czegoś w rodzaju średniej i odchylenia standardowego i możesz utworzyć metrykę odległości, aby określić, jak daleko jesteś od celu, możesz użyć optymalizacji do wyszukiwania zmiennych wejściowych, które dają pożądany wynik wartości.

Jako brzydki przykład, w którym szukamy losowego jednorodnego wektora o długości 100, który ma średnią = 0 i odchylenie standardowe = 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

Czy są jakieś programy w R, które to robią?

Pakiet Runuran R zawiera wiele metod generowania losowych zmiennych. Wykorzystuje biblioteki C z projektu UNU.RAN (Universal Non-Uniform RAndom Number generator). Moja własna wiedza na temat generowania zmiennych losowych jest ograniczona, ale winieta Runuran zapewnia ładny przegląd. Poniżej znajdują się dostępne metody w pakiecie Runuran, zaczerpnięte z winiety:

Ciągłe dystrybucje:

• Adaptacyjne odrzucanie próbek
• Odwrotne odwrócenie transformowanej gęstości
• Wielomianowa interpolacja odwrotnego CDF
• Prosta metoda stosunku mundurów
• Odrzucone przekształcenie gęstości

Dyskretne rozkłady:

• Dyskretna automatyczna inwersja odrzucania
• Metoda Alias-Urn
• Metoda tabeli dyskretnej inwersji

Rozkłady wielowymiarowe:

• Algorytm Hit-and-Run z metodą Ratio-of-Uniforms
• Metoda wielowymiarowego naiwnego stosunku mundurów

Przykład:

Na przykład, załóżmy, że chcesz wygenerować rozkład normalny ograniczony od 0 do 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)

Ta urnorm()funkcja jest wygodną funkcją owijania. Uważam, że za kulisami używa metody interpolacji wielomianowej odwrotnej CDF, ale nie jestem pewien. Dla czegoś bardziej złożonego, powiedzmy, dyskretny rozkład normalny ograniczony od 0 do 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)

Wygląda na to, że pakiet R spełniający twoje wymagania został opublikowany wczoraj! simstudy Keith Goldfeld

Symuluje zestawy danych w celu poznania technik modelowania lub lepszego zrozumienia procesów generowania danych. Użytkownik określa zestaw relacji między zmiennymi towarzyszącymi i generuje dane na podstawie tych specyfikacji. Ostateczne zestawy danych mogą reprezentować dane z randomizowanych prób kontrolnych, powtarzalnych pomiarów (podłużnych) i losowych prób skupień. Brak można wygenerować za pomocą różnych mechanizmów (MCAR, MAR, NMAR).

To odpowiedź przychodzi tak późno, że prawdopodobnie nie ma ona znaczenia, ale zawsze istnieje rozwiązanie MCMC. Mianowicie, aby rzutować gęstość połączenia próbki na kolektor zdefiniowany przez ograniczenia, na przykład Jedynym problemem jest wówczas symulacja wartości w tym kolektorze, tj. znalezienie parametryzacji właściwego wymiaru. Artykuł z 2015 roku autorstwa Bornna, Shepharda i Solgi analizuje ten bardzo problem (z ciekawą, jeśli nie ostateczną odpowiedzią ).

$$\prod_{i=1}^n f(x_i)$$ $$\sum_{i=1}^n x_i=\mu_0\qquad\sum_{i=1}^n x_i^2=\sigma_0^2$$

Ta odpowiedź rozważa inne podejście do przypadku, w którym chcesz zmusić wariaty do położenia się w określonym zakresie i dodatkowo dyktować średnią i / lub wariancję.

Ogranicz naszą uwagę do interwału jednostkowego . Użyjmy średniej ważonej dla ogólności, więc napraw niektóre wagi pomocą lub ustaw jeśli chcesz standardowej wagi. Załóżmy, że ilości i reprezentują odpowiednio pożądaną (ważoną) średnią i (ważoną) wariancję. Górna granica jest konieczna, ponieważ jest to maksymalna możliwa wariancja w jednostkowym przedziale. Jesteśmy zainteresowani rysowaniem niektórych wariantów z z tymi ograniczeniami momentu.$[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]$

Najpierw narysujemy kilka odmian z dowolnego rozkładu, np. . Rozkład ten wpłynie na kształt rozkładu końcowego. Następnie ograniczamy je do przedziału jednostek za pomocą funkcji logistycznej:$y_1,...,y_N$$N(0,1)$$[0,1]$

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

Jednak zanim to zrobimy, jak widać w powyższym równaniu, przekształcamy za pomocą translacji i skali . Jest to analogiczne do pierwszego równania w odpowiedzi @ Macro. Sztuką jest wybrać i tak, że przekształcone zmienne mieć pożądaną chwilę (y). Oznacza to, że do przechowywania potrzebujemy jednego lub obu następujących elementów: $y_k$$h$$v$$h$$v$$x_1,...,x_N$

$$\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$$

Analityczne odwrócenie tych równań dla i nie jest możliwe, ale wykonanie liczbowe jest proste, zwłaszcza że pochodne w odniesieniu do i są łatwe do obliczenia; zajmuje tylko kilka iteracji metody Newtona.$v$$h$$v$$h$

Jako pierwszy przykład załóżmy, że zależy nam jedynie na ograniczeniu średniej ważonej, a nie na wariancji. Fix , , , . Następnie dla leżących u podstaw rozkładów , i otrzymujemy odpowiednio następujące histogramy i takie, że średnia zmiennych wynosi dokładnie (nawet dla małych ):v = 1 w k = 1 / N N = 200000 N ( 0 , 1 ) N ( 0 , 0,1 ) Unif ( 0 , 1 ) 0,8 $\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$

Następnie ograniczmy zarówno średnią, jak i wariancję. Weźmy , , i rozważmy trzy pożądane odchylenia standardowe . Stosując ten sam podstawowy rozkład , tutaj są histogramy dla każdego:w k = 1 / N N = 2000 σ = 0,1 , 0,05 , 0,01 N ( 0 , 1 $\mu=0.2$$w_k=1/N$$N=2000$$\sigma=0.1,0.05,0.01$$N(0,1)$

Pamiętaj, że mogą one wyglądać nieco rozproszone w wersji beta, ale nie są.

W mojej odpowiedzi tutaj wymieniłem trzy pakiety R do zrobienia tego:

• SimCorrMix
• SimMultiCorrData
• simrel