Как смоделировать данные, которые удовлетворяют определенным ограничениям, таким как наличие определенного среднего значения и стандартного отклонения?

56

Этот вопрос мотивирован моим вопросом о метаанализе . Но я полагаю, что это также было бы полезно при обучении контекстов, в которых вы хотите создать набор данных, который точно отражает существующий опубликованный набор данных.

Я знаю, как генерировать случайные данные из данного распределения. Так, например, если я читаю о результатах исследования, которое имело:

среднее значение 102,
стандартное отклонение 5,2, и
размер выборки 72.

Я мог бы сгенерировать аналогичные данные, используя rnormR. Например,

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

Конечно, среднее значение и SD не будут точно равны 102 и 5.2 соответственно:

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

В целом меня интересует, как симулировать данные, которые удовлетворяют набору ограничений. В приведенном выше случае составляющими являются размер выборки, среднее значение и стандартное отклонение. В других случаях могут быть дополнительные ограничения. Например,

минимум и максимум в данных или базовой переменной могут быть известны.
известно, что переменная принимает только целочисленные значения или только неотрицательные значения.
данные могут включать несколько переменных с известными взаимными корреляциями.

Вопросов

В общем, как я могу симулировать данные, которые точно удовлетворяют набору ограничений?
Есть ли статьи, написанные по этому поводу? Есть ли в R программы, которые это делают?
Для примера, как можно и нужно моделировать переменную, чтобы она имела определенное среднее значение и sd?

r dataset simulation random-generation Джером англим
источник

1

Почему вы хотите, чтобы они были в точности как опубликованные результаты? Разве это не оценки среднего значения и стандартного отклонения для популяции с учетом их выборки данных. Учитывая неопределенность в этих оценках, кто скажет, что приведенная выше выборка не соответствует их наблюдениям?

Восстановить Монику - Дж. Симпсон

4

Поскольку этот вопрос, по-видимому, собирает ответы, которые не попадают в точку (ИМХО), я хотел бы отметить, что концептуально ответ прост: ограничения на равенство рассматриваются как маргинальные распределения, а ограничения на неравенство являются многомерными аналогами усечения. Усечение относительно легко обрабатывать (часто с выборкой отклонения); более сложная проблема состоит в том, чтобы найти способ выборки этих предельных распределений. Это означает либо выборку маргинальных значений с учетом распределения и ограничения, либо интеграцию, чтобы найти предельное распределение и выборку из него.

whuber

4

Кстати, последний вопрос тривиален для семейств распределений масштаба расположения. Например, x<-rnorm(72);x<-5.2*(x-mean(x))/sd(x)+102делает трюк.

whuber

1

@whuber, как намекает кардинал в комментарии к моему ответу (в котором упоминается этот «трюк») и комментарии к другому ответу - этот метод, как правило, не будет хранить переменные в пределах одного и того же семейства распределений, поскольку вы делите по стандартному отклонению выборки.

Макрос

5

@Macro Это хороший момент, но, возможно, лучший ответ: «Конечно, они не будут иметь одинаковое распространение»! Распределение, которое вы хотите, является условием распределения по ограничениям. В общем, это не будет из той же семьи, что и родительский дистрибутив. Например, каждый элемент выборки размера 4 со средним значением 0 и SD 1, взятым из нормального распределения, будет иметь почти одинаковую вероятность на [-1,5, 1,5], потому что условия устанавливают верхнюю и нижнюю границы возможных значений.

whuber

26

В целом, чтобы среднее значение и дисперсия в вашей выборке были в точности равны предварительно заданному значению, вы можете соответствующим образом изменить и масштабировать переменную. В частности, если - образец, то новые переменные $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}$

где - среднее значение выборки, а - выборочная дисперсия такова, что среднее значение выборки для равно а их выборочная дисперсия равна точно . Аналогично построенный пример может ограничить диапазон - $\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)$

создаст набор данных который ограничен интервалом . $B_1, ..., B_n$ $(a,b)$

Примечание. Эти типы сдвига / масштабирования, как правило, изменяют семейство данных, даже если исходные данные принадлежат семейству масштабов местоположения.

В контексте нормального распределенияmvrnorm функция R позволяет моделировать нормальные (или многомерные нормальные) данные с заранее определенным образцом средней / ковариацией путем установкой empirical=TRUE. В частности, эта функция моделирует данные из условного распределения нормально распределенной переменной, учитывая, что среднее значение выборки и (со) дисперсия равна предварительно заданному значению . Обратите внимание, что получающиеся маргинальные распределения не являются нормальными, как указал @whuber в комментарии к основному вопросу.

Вот простой одномерный пример, где среднее значение выборки (из выборки ) ограничено равным 0, а стандартное отклонение выборки равно 1. Мы можем видеть, что первый элемент гораздо больше похож на равномерное распределение, чем нормальное распределение: $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")

$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$ введите описание изображения здесь

макрос
источник

1

не будет нормально распределен, хотя они могут быть примерно так , если размер выборки является большим. Первый комментарий к ответу @ Шона ссылается на это.

Z_{i}

$Z_i$

кардинал

1

Ну, это довольно естественно , что хотят сделать ... и часто не вызывает слишком много проблем.

кардинал

1

+1. В примере, кстати , униформа является точным ответом. (Кажущийся спад в конце графика - это артефакт того, как R рисует гистограммы.)

whuber

1

@whuber, спасибо за мотивирование этого примера. Принимая во внимание тот факт, что предельные распределения меняются, когда вы определяете среднее / дисперсионное значение выборки, кажется, что лучший «ответ» в духе вопроса ФП - просто моделировать данные со средним / дисперсионным значением популяции, равным представленному в выборке. количества (как предполагает сам ОП), не так ли? Таким образом, вы получаете количество образцов, «похожее» на желаемое, а предельные распределения - это то, что вы хотели.

Макрос

1

@whuber, если у вас нормальная выборка, то у есть -распределение, да? Рассматриваемая «новая» переменная будет просто линейной комбинацией .

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

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

t

$t$

T_{i}

$T_i$

Макрос

22

Что касается вашего запроса на документы, есть:

Чаттерджи, С. и Фират, А. (2007). Генерация данных с идентичной статистикой, но с разной графикой: продолжение набора данных Anscombe. Американский статистик, 61 , 3, с. 248-254.

Это не совсем то, что вы ищете, но может служить зерном для мельницы.

Есть еще одна стратегия, о которой никто, кажется, не упомянул. Можно генерировать (псевдо) случайных данных из набора размера , чтобы весь набор удовлетворял ограничениям, при условии, что оставшиеся данных фиксированы в соответствующих значениях. Требуемые значения должны быть решаемы с помощью системы уравнений, алгебры и некоторого количества коленчатого смазочного материала. $N-k$ $N$ $k$ $k$ $k$

Например, чтобы сгенерировать набор из данных из нормального распределения, который будет иметь заданное среднее значение выборки, , и дисперсию , вам нужно будет зафиксировать значения двух точек: и . Поскольку среднее значение выборки: должно быть: Пример дисперсии: таким образом (после замены вышеуказанного на , срыв / распределение и перестановка ... ) мы получили: $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$ Если мы берем , , и как отрицание RHS, мы можем решить для используя квадратную формулу . Например, в , следующий код может быть использован:

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

Есть несколько вещей, чтобы понять об этом подходе. Во-первых, это не гарантирует работу. Например, возможно, что ваши исходные данные таковы, что не существует значений и , из-за которых дисперсия полученного набора будет равна . Рассмотреть возможность: $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

Во-вторых, в то время как стандартизация делает предельные распределения всех ваших вариантов более равномерными, этот подход влияет только на последние два значения, но делает их предельные распределения перекошенными:

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
}

(вставить сюжет)

В-третьих, полученный образец может выглядеть не совсем нормально; может показаться, что у него есть «выбросы» (т. е. точки, полученные из процесса генерации данных, отличного от остальных), поскольку это, по сути, так. Это менее вероятно, будет проблемой с большими размерами выборки, так как статистика выборки из сгенерированных данных должна сходиться к требуемым значениям и, следовательно, нуждается в меньшей корректировке. С небольшими выборками вы всегда можете объединить этот подход с алгоритмом принятия / отклонения, который повторяет попытку, если сгенерированный образец имеет статистику формы (например, асимметрию и эксцесс), которые находятся за пределами допустимых границ (см., Комментарий @ cardinal ), или расширяете этот подход для генерации выборки с фиксированным средним, дисперсией, асимметрией иkurtosis (я оставлю алгебру до вас, хотя). Кроме того, вы можете сгенерировать небольшое количество выборок и использовать одну с наименьшей (скажем) статистикой Колмогорова-Смирнова.

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

(добавить сюжет)

Gung - Восстановить Монику
источник

10

Общий метод - «Метод отклонения», при котором вы просто отклоняете результаты, которые не соответствуют вашим ограничениям. Если у вас нет какого-либо руководства (например, MCMC), то вы можете генерировать множество случаев (в зависимости от вашего сценария), которые отклоняются!

Там, где вы ищете что-то вроде среднего и стандартного отклонения, и вы можете создать метрику расстояния для определения расстояния до цели, вы можете использовать оптимизацию для поиска входных переменных, которые дают желаемый результат. ценности.

В качестве ужасного примера, где мы будем искать случайный равномерный вектор длиной 100, который имеет среднее значение = 0 и стандартное отклонение = 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

Шон
источник

7

Ограничения, которые возникают с нулевой вероятностью, трудно удовлетворить. ;-) Для конкретного примера под рукой соответствующий сдвиг и расширение легко достигают поставленных целей, хотя может потребоваться более глубокий анализ, чтобы посмотреть, как такая операция нарушает распределение данных.

кардинал

Благодарю. Конечно, было бы легко отклонить наблюдения меньше минимума и больше максимума. И я вижу, как вы можете определить это как проблему оптимизации. Было бы здорово увидеть некоторые примеры или, может быть, есть предложения о том, что читать дальше.

Jeromy Anglim

1

@ cardinal - согласен. Нужно посмотреть на распределение (то есть гистограмму) как симулированных входных чисел, так и выходных, поскольку иногда они могут выглядеть очень странно!

Шон

9

Есть ли в R программы, которые это делают?

Пакет Runuran R содержит много методов для генерации случайных переменных. Он использует библиотеки C из проекта UNU.RAN (Универсальный генератор неоднородных случайных чисел). Мои собственные знания в области генерации случайных вариаций ограничены, но рунурская виньетка дает хороший обзор. Ниже приведены доступные методы в пакете Runuran, взятые из виньетки:

Непрерывные распределения:

Адаптивная выборка отклонения
Обратное преобразование измененной плотности
Полиномиальная интерполяция обратных CDF
Простой метод отношения форм
Преобразование трансформированной плотности

Дискретные распределения:

Дискретная автоматическая инверсия отклонения
Метод псевдонимов-урн
Метод направляющих таблиц для дискретной инверсии

Многомерные распределения:

Алгоритм удара и бега с помощью метода Ratio-of-Uniforms
Метод многомерного наивного соотношения форм

Пример:

Для быстрого примера, предположим, что вы хотите сгенерировать нормальное распределение, ограниченное от 0 до 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)

Эта urnorm()функция является удобной функцией-оберткой. Я считаю, что за кулисами он использует метод полиномиальной интерполяции обратного CDF, но я не уверен. Для чего-то более сложного, скажем, дискретного нормального распределения, ограниченного от 0 до 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
источник

3

Кажется, есть пакет R, отвечающий вашим требованиям, опубликованный только вчера! симстуди Кит Голдфельд

Имитирует наборы данных, чтобы изучить методы моделирования или лучше понять процессы генерации данных. Пользователь задает набор отношений между ковариатами и генерирует данные на основе этих спецификаций. Окончательные наборы данных могут представлять данные из рандомизированных контрольных испытаний, схем повторных измерений (продольных) и кластерных рандомизированных испытаний. Отсутствие может быть получено с использованием различных механизмов (MCAR, MAR, NMAR).

Tyelcie
источник

1

Ни в виньетке, ни на домашней странице программы точное упоминание ограничений не упоминается. Как вы думаете, почему этот пакет соответствует требованию рисования из условных распределений?

гг

2

Это ответ, приходящий так поздно, что он, по-видимому, не имеет смысла, но всегда есть решение MCMC для этого вопроса. А именно, чтобы спроецировать общую плотность образца на многообразии, определяемом ограничениями, например, Единственная проблема заключается в моделировании значений по этому многообразию, т. е. в поиске параметризации правильного измерения. Статья Bornn, Shephard и Solgi, опубликованная в 2015 году, исследует эту проблему (с интересным, если не окончательным ответом ).

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

Сиань
источник

2

Этот ответ рассматривает другой подход к случаю, когда вы хотите, чтобы переменные лежали в указанном диапазоне и дополнительно определяли среднее значение и / или дисперсию.

Ограничим наше внимание единичным интервалом . Давайте используем взвешенное среднее для общности, поэтому зафиксируем некоторые веса с помощью или установите если вы хотите стандартное взвешивание. Предположим, что величины и представляют желаемую (взвешенную) среднюю и (взвешенную) дисперсию соответственно. Верхняя граница необходима, потому что это максимально возможная дисперсия на единичном интервале. Нам интересно нарисовать некоторые переменные из с этими моментными ограничениями. $[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]$

Сначала мы нарисуем некоторые переменные из любого распределения, например . Это распределение будет влиять на форму окончательного распределения. Затем мы ограничим их единичным интервалом с помощью логистической функции: $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)}}$

Однако, прежде чем мы это сделаем, как видно из приведенного выше уравнения, мы преобразуем с переводом и масштабом . Это аналогично первому уравнению в ответе @ Macro. Хитрость заключается в том, чтобы выбрать и чтобы преобразованные переменные имели желаемый момент (ы). То есть нам требуется одно или оба из следующих условий: $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$

Аналитически инвертировать эти уравнения для и не представляется возможным, но сделать это численно просто, особенно потому, что производные по и легко вычислить; это займет всего несколько итераций метода Ньютона. $v$ $h$ $v$ $h$

В качестве первого примера, скажем, мы заботимся только об ограничении взвешенного среднего значения, а не о дисперсии. Fix , , , . Тогда для базовых распределений , и мы получаем следующие гистограммы, соответственно, и такие, что среднее значение переменных составляет точно (даже для маленьких ): $\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$

Далее, давайте ограничим среднее и дисперсию. Возьмите , , и рассмотрите три желаемых стандартных отклонения . Используя одно и то же базовое распределение , вот гистограммы для каждого: $\mu=0.2$ $w_k=1/N$ $N=2000$ $\sigma=0.1,0.05,0.01$ $N(0,1)$

Обратите внимание, что они могут выглядеть немного бета-распределенными, но это не так.

Ян Хинкс
источник

1

В моем ответе здесь я перечислил три пакета R для этого:

abalter
источник

Там должен быть какой-то формат для ссылки на ссылки. Это должен быть комментарий?

abalter

Как смоделировать данные, которые удовлетворяют определенным ограничениям, таким как наличие определенного среднего значения и стандартного отклонения?

Вопросов

Ответы: