Como simular dados que satisfazem restrições específicas, como média e desvio padrão específicos?

56

Esta questão é motivada pela minha pergunta sobre meta-análise . Mas imagino que também seria útil no ensino de contextos em que você deseja criar um conjunto de dados que espelhe exatamente um conjunto de dados publicado existente.

Eu sei como gerar dados aleatórios a partir de uma determinada distribuição. Por exemplo, se eu li sobre os resultados de um estudo que teve:

uma média de 102,
um desvio padrão de 5,2 e
uma amostra de 72.

Eu poderia gerar dados semelhantes usando rnormem R. Por exemplo,

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

Obviamente, a média e o DP não seriam exatamente iguais a 102 e 5,2, respectivamente:

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

Em geral, estou interessado em como simular dados que satisfaçam um conjunto de restrições. No caso acima, as restrições são tamanho da amostra, média e desvio padrão. Em outros casos, pode haver restrições adicionais. Por exemplo,

um mínimo e um máximo nos dados ou na variável subjacente podem ser conhecidos.
sabe-se que a variável assume apenas valores inteiros ou apenas valores não negativos.
os dados podem incluir várias variáveis com inter-correlações conhecidas.

Questões

Em geral, como posso simular dados que satisfazem exatamente um conjunto de restrições?
Existem artigos escritos sobre isso? Existem programas no R que fazem isso?
Por uma questão de exemplo, como posso e devo simular uma variável para que ela tenha uma média e um sd específicos?

r dataset simulation random-generation Jeromy Anglim
fonte

11

Por que você deseja que eles sejam exatamente como os resultados publicados? Essas estimativas da média e desvio padrão da população não são dadas pela amostra de dados. Dada a incerteza nessas estimativas, quem pode dizer que a amostra que você mostra acima não é consistente com as observações deles?

Reintegrar Monica - G. Simpson

4

Como essa pergunta parece estar coletando respostas que não atendem ao objetivo (IMHO), gostaria de salientar que conceitualmente a resposta é direta: restrições de igualdade são tratadas como distribuições marginais e restrições de desigualdade são análogos multivariados de truncamento. O truncamento é relativamente fácil de manusear (geralmente com amostragem por rejeição); o problema mais difícil é encontrar uma maneira de provar essas distribuições marginais. Isso significa amostragem de marginais, dada a distribuição e a restrição, ou integração para encontrar a distribuição marginal e a amostragem.

whuber

4

BTW, a última pergunta é trivial para famílias de distribuição em escala de localização. Por exemplo, x<-rnorm(72);x<-5.2*(x-mean(x))/sd(x)+102faz o truque.

whuber

11

@whuber, como o cardeal alude em um comentário à minha resposta (que menciona esse "truque") e um comentário a outra resposta - esse método, em geral, não mantém as variáveis dentro da mesma família distributiva, já que você está dividindo pelo desvio padrão da amostra.

Macro

5

@ Macro Este é um bom ponto, mas talvez a melhor resposta seja: "é claro que eles não terão a mesma distribuição"! A distribuição que você deseja é a distribuição condicional às restrições. Em geral, isso não será da mesma família que a distribuição dos pais. Por exemplo, cada elemento de uma amostra de tamanho 4 com média 0 e DP 1 extraída de uma distribuição normal terá uma probabilidade quase uniforme em [-1,5, 1,5], porque as condições colocam limites superior e inferior nos valores possíveis.

whuber

26

Em geral, para tornar sua amostra média e variação exatamente igual a um valor pré-especificado, você pode alterar e dimensionar adequadamente a variável. Especificamente, se é uma amostra, então as novas variáveis $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}$

onde é a média da amostra $\overline{X} = \frac{1}{n} \sum_{i=1}^{n} X_i$ é a variância da amostra são tais que a média da amostra do's é exactamentee a sua variância da amostra é exactamente. Um exemplo de construção semelhante pode restringir o intervalo - $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)$

produzirá um conjunto de dados restrito ao intervalo . $B_1, ..., B_n$ $(a,b)$

Nota: Esses tipos de mudança / dimensionamento, em geral, alteram a família de distribuição dos dados, mesmo que os dados originais venham de uma família de escala de local.

Dentro do contexto da distribuição normal, a mvrnormfunção in R permite simular dados normais (ou normais multivariados) com uma média / covariância amostral pré-especificada por configuração empirical=TRUE. Especificamente, essa função simula dados da distribuição condicional de uma variável normalmente distribuída, considerando que a média da amostra e (co) variância é igual a um valor pré-especificado . Observe que as distribuições marginais resultantes não são normais, como apontado por @whuber em um comentário à questão principal.

Aqui está um exemplo simples e univariado, em que a média da amostra (de uma amostra de ) é restrita a 0 e o desvio padrão da amostra é 1. Podemos ver que o primeiro elemento é muito mais semelhante a uma distribuição uniforme do que uma normal distribuição: $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")

$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$ insira a descrição da imagem aqui

Macro
fonte

11

O não será normalmente distribuído, embora possa ser aproximadamente isso se o tamanho da amostra for grande. O primeiro comentário à resposta de @ Sean faz alusão a isso.

Z_{i}

$Z_i$

cardeal

11

Bem, isso é uma coisa muito natural querer fazer ... e muitas vezes não causa também muita dificuldade.

cardeal

11

+1. No exemplo, o uniforme é a resposta exata , a propósito. (A aparente dropoff nas extremidades da trama é um artefacto de R como empates histogramas.)

whuber

11

@whuber, obrigado por motivar este exemplo. Dado o fato de que as distribuições marginais mudam quando você condiciona a média / variância da amostra, parece que a melhor "resposta" no espírito da pergunta do OP é apenas simular dados com média / variância populacional igual à relatada como amostra quantidades (como sugerido pelo próprio OP), não é? Dessa forma, você obtém quantidades de amostra "semelhantes" às desejadas e as distribuições marginais são o que você queria que elas fossem.

Macro

11

@ whuber, se sua amostra é normal, então tem uma distribuição , sim? A "nova" variável em questão será apenas uma combinação linear de .

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

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

t

$t$

T_{i}

$T_i$

Macro

22

Em relação à sua solicitação de trabalhos, há:

Chatterjee, S. & Firat, A. (2007). Gerando dados com estatísticas idênticas, mas com gráficos diferentes: Um acompanhamento do conjunto de dados Anscombe. The American Statistician, 61 , 3, pp. 248-254.

Isso não é exatamente o que você está procurando, mas pode servir como fonte para o moinho.

Há outra estratégia que ninguém parece ter mencionado. É possível gerar dados aleatórios (pseudo) a partir de um conjunto de tamanho modo que todo o conjunto atenda às restrições , desde que os dados restantes sejam fixados em valores apropriados. Os valores exigidos devem ser solucionáveis com um sistema de equações, álgebra e alguma graxa de cotovelo. $N-k$ $N$ $k$ $k$ $k$

Por exemplo, para gerar um conjunto de dados a partir de uma distribuição normal que terá uma dada média da amostra, , e a variância, , será necessário corrigir os valores de dois pontos: e . Como a média da amostra é: deve ser: A variação da amostra é: assim (após substituir o acima por , frustrar / distribuir e reorganizar ... ) Nós temos: $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$ Se tomarmos , , e como a negação da RHS, podemos resolver em usando a fórmula quadrática . Por exemplo, em , o seguinte código pode ser usado:

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

Há algumas coisas para entender sobre essa abordagem. Primeiro, não é garantido que funcione. Por exemplo, é possível que as suas iniciais dados são tais que não os valores e existentes que fará com que a variância dos resultantes definidas iguais . 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

Segundo, enquanto a padronização torna as distribuições marginais de todas as suas variáveis mais uniformes, essa abordagem afeta apenas os dois últimos valores, mas distorce suas distribuições marginais:

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

(inserir 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
}

(inserir gráfico)

Terceiro, a amostra resultante pode não parecer muito normal; pode parecer que tem 'outliers' (ou seja, pontos que provêm de um processo de geração de dados diferente do restante), já que esse é basicamente o caso. É menos provável que seja um problema com tamanhos de amostra maiores, pois as estatísticas da amostra dos dados gerados devem convergir para os valores necessários e, portanto, precisam de menos ajustes. Com amostras menores, você sempre pode combinar essa abordagem com um algoritmo de aceitação / rejeição que tenta novamente se a amostra gerada tiver estatísticas de forma (por exemplo, assimetria e curtose) que estão fora dos limites aceitáveis (cf., comentário do @ cardinal ) ou estender esta abordagem para gerar uma amostra com média fixa, variância, assimetria ecurtose (eu vou deixar a álgebra para você, no entanto). Como alternativa, você pode gerar um pequeno número de amostras e usar aquela com a menor (digamos) estatística 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

(adicionar plot)

- Reinstate Monica
fonte

10

A técnica geral é o 'Método de rejeição', no qual você apenas rejeita resultados que não atendem às suas restrições. A menos que você tenha algum tipo de orientação (como o MCMC), poderá gerar muitos casos (dependendo do cenário) que são rejeitados!

Onde você está procurando algo como um desvio médio e padrão e pode criar algum tipo de métrica de distância para dizer a que distância está de seu objetivo, pode usar a otimização para procurar as variáveis de entrada que fornecem a saída desejada valores.

Como um exemplo feio, onde procuraremos um vetor uniforme aleatório com comprimento 100 que tenha média = 0 e desvio padrão = 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
fonte

7

Restrições que ocorrem com probabilidade zero são difíceis de satisfazer. ;-) Para o exemplo específico em questão, uma mudança e uma dilatação apropriadas atingem facilmente os objetivos declarados , embora se queira analisar um pouco mais profundamente para ver como a distribuição dos dados é perturbada por essa operação.

cardeal

Obrigado. Certamente, seria fácil rejeitar observações menores que o mínimo e maiores que o máximo. E eu posso ver como você pode defini-lo como um problema de otimização. Seria ótimo ver alguns exemplos ou talvez ter algumas sugestões sobre o que ler a seguir.

21412 Jeremy

11

@ cardinal - concordou. Deve-se olhar para as distribuições (isto é, um histograma) dos números simulados de entrada e saída, pois às vezes eles podem parecer muito estranhos!

21412 Sean

9

Existem programas no R que fazem isso?

O pacote Runuran R contém muitos métodos para gerar variáveis aleatórias. Ele usa bibliotecas C do projeto UNU.RAN (gerador de número de uniforme universal não uniforme). Meu próprio conhecimento do campo da geração aleatória de variáveis é limitado, mas a vinheta de Runuran fornece uma boa visão geral. Abaixo estão os métodos disponíveis no pacote Runuran, retirados da vinheta:

Distribuições contínuas:

Amostragem de rejeição adaptativa
Rejeição de densidade transformada inversa
Interpolação polinomial de CDF inverso
Método Simples de Proporção de Uniformes
Rejeição de densidade transformada

Distribuições discretas:

Inversão discreta de rejeição automática
Método Alias-Urn
Método da tabela de guia para inversão discreta

Distribuições multivariadas:

Algoritmo Hit-and-Run com método de proporção de uniformes
Método multivariado de proporção de uniformes ingênuos

Exemplo:

Para um exemplo rápido, suponha que você deseje gerar uma distribuição Normal limitada entre 0 e 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)

A urnorm()função é uma função de invólucro conveniente. Eu acredito que nos bastidores ele usa o método de interpolação polinomial de CDF inverso, mas não tenho certeza. Para algo mais complexo, digamos, uma distribuição Normal discreta limitada entre 0 e 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
fonte

3

Parece que há um pacote R atendendo aos seus requisitos publicado ontem! simstudy Por Keith Goldfeld

Simula conjuntos de dados para explorar técnicas de modelagem ou entender melhor os processos de geração de dados. O usuário especifica um conjunto de relacionamentos entre covariáveis e gera dados com base nessas especificações. Os conjuntos de dados finais podem representar dados de ensaios clínicos randomizados, desenhos de medidas repetidas (longitudinais) e ensaios clínicos randomizados por cluster. A falta pode ser gerada usando vários mecanismos (MCAR, MAR, NMAR).

Tyelcie
fonte

11

Nem na vinheta nem na home page do programa são mencionadas exatamente as restrições. Por que você acha que este pacote atende ao requisito de extrair distribuições condicionais?

gg

2

Esta é uma resposta que chega tão tarde que, presumivelmente, não tem sentido, mas sempre há uma solução para o MCMC. Ou seja, para projetar a densidade da junta da amostra no coletor definido pelas restrições, por exemplo A única questão é simular valores sobre essa variedade, ou seja, encontrar uma parametrização da dimensão correta. Um artigo de 2015 de Bornn, Shephard e Solgi estuda esse mesmo problema (com uma resposta interessante, se não definitiva ).

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

2

Esta resposta considera outra abordagem para o caso em que você deseja forçar as variáveis a ficarem em um intervalo especificado e, adicionalmente, ditar a média e / ou variação.

Restrinja nossa atenção ao intervalo da unidade . Vamos usar uma média ponderada para generalidade, portanto, fixe alguns pesos com ou defina se desejar uma ponderação padrão. Suponha que as quantidades e representem a média (ponderada) desejada e a variação (ponderada), respectivamente. O limite superior em é necessário porque essa é a variação máxima possível em um intervalo de unidades. Estamos interessados em desenhar algumas variáveis de com essas restrições 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]$

Primeiro, desenhamos algumas variáveis de qualquer distribuição, como . Essa distribuição afetará o formato da distribuição final. Em seguida, as restringimos ao intervalo unitário usando uma função 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)}}$

Antes de fazermos isso, no entanto, como visto na equação acima, transformamos os 's com a tradução escala . Isso é análogo à primeira equação na resposta do @ Macro. O truque é agora de escolher e , para que as variáveis transformadas têm o momento desejado (s). Ou seja, exigimos um ou ambos dos seguintes : $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$

A inversão dessas equações para e analiticamente não é viável, mas fazê-lo numericamente é simples, especialmente porque derivadas com relação a e são fáceis de calcular; são necessárias apenas algumas iterações do método de Newton. $v$ $h$ $v$ $h$

Como primeiro exemplo, digamos que apenas nos preocupamos em restringir a média ponderada e não a variação. Correção , , , . Então, para as distribuições subjacentes , e , terminamos com os seguintes histogramas, respectivamente, e de modo que a média das variáveis seja exatamente (mesmo para pequeno ): $\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$

Em seguida, vamos restringir a média e a variância. Tome , , e considere os três desvios padrão desejados . Usando a mesma distribuição subjacente , aqui estão os histogramas para cada: $\mu=0.2$ $w_k=1/N$ $N=2000$ $\sigma=0.1,0.05,0.01$ $N(0,1)$

Observe que eles podem parecer um pouco distribuídos em beta, mas não são.

Ian Hincks
fonte

1

Na minha resposta aqui , listei três pacotes R para fazer isso:

abalter
fonte

Deve haver algum formato para um link para referências. Deve ser um comentário?

abalter 13/08

Como simular dados que satisfazem restrições específicas, como média e desvio padrão específicos?

Questões

Respostas: