Gere números aleatórios normais dependentes de distribuição idêntica com soma pré-especificada

Como eu gero números aleatórios normais aleatoriamente distribuídos de forma idêntica, mas não independentes, de modo que sua soma caia dentro de um intervalo pré-especificado com probabilidade ?$n$$[a,b]$$p$

(Essa pergunta é motivada pela geração de uma caminhada aleatória que termina em um ponto pré-especificado: afinal, um processo aleatório não é tão aleatório (determinístico) . Como uma variável aleatória contínua tem uma probabilidade zero de atingir um número exato, fazemos a segunda melhor coisa e peça um intervalo inteiro para terminar.)

---

EDIT: A geração de amostras a partir da distribuição Gaussiana singular foi proposta como duplicada, que por sua vez é fechada como duplicata de Gerar números aleatórios normalmente distribuídos com matriz de covariância definida não positiva . Concordo que ambos são úteis. No entanto, o objetivo da pergunta atual (mais especificamente, da resposta) é primeiro descobrir que podemos usar uma distribuição normal multivariada para abordar a questão e, segundo, que tipo de matriz de covariância funciona. Como obter amostras de uma distribuição com essa covariância é uma terceira etapa, na qual os threads vinculados são úteis.

Geraremos normais multivariados com e forma que sua soma seja satisfatória nossa condição. Seja .$X\sim MN(\mu, \Sigma)$$\mu\in\mathbb{R}^n$$\Sigma\in\mathbb{R}^{n\times n}$$Z=X_1 + \dots + X_n$

Como um meio comum, escolhemos

$$\mu_1 = \dots = \mu_n = \frac{a+b}{2n}.$$

Para que com probabilidade , seu desvio padrão deva cumprir$Z\in[a,b]$$p$

$$\sigma_Z = \frac{b-a}{q_\alpha},$$

onde é o quantil normal padrão para o nível , aqui .$q_\alpha$$\alpha$$\alpha=1-\frac{1-p}{2}$

Agora precisamos especificar . Temos muita margem de manobra aqui. Vamos supor que queremos que a cada seja e a covariância seja para . A chave para criar um "bom" é esta resposta anterior por probabilityislogic . Isso resulta que a soma de nossos s tem variação$\Sigma$$X_i$$\sigma^2$$\text{cov}(X_i,X_j)=\tau$$i\neq j$$\Sigma$$X_i$

$$n\sigma^2 + n(n-1)\tau$$

então precisamos disso

$$n\sigma^2 + n(n-1)\tau = \frac{b-a}{q_\alpha}.$$

Também precisamos garantir que seja definitivo positivo, mas isso não é muito difícil. A maneira mais fácil de fazer isso é garantir que todas as entradas em sejam positivas, por exemplo, definindo$\Sigma$$\Sigma$

$$\sigma^2 := \frac{\sigma_Z^2}{2n}, \quad \tau := \frac{\sigma_Z^2}{2n(n-1)},$$

mas isso fornece valores muito pequenos e somas e trajetórias cumulativas muito chatas:

Menos chato é definir

$$\sigma^2 := 1, \quad \tau := \frac{1}{n-1}\big(\frac{\sigma_Z^2}{n}-\sigma^2\big),$$

que produz trajetórias muito mais interessantes:

Observe que definir isso de fato produz uma matriz de covariância válida, porque é então da forma , a saber$\Sigma$$\Sigma_{ij} = m(i-j)$

$$m(0) = \sigma^2, \quad m(j) = \tau\text{ for }j>0,$$

e nós temos isso

$$\sum_{j>0} |m(j)| = (n-1)|\tau| = \big|\frac{\sigma_Z^2}{n}-\sigma^2\big| = \big|\frac{\sigma_Z^2}{n}-1\big| < 1 = \sigma^2 = m(0),$$

que é uma condição suficiente para ser estritamente positivo definido pela Wikipedia (ponto 7 em "Outras propriedades" ).$\Sigma$

Código R abaixo, mas primeiro, vá para cima e vote na resposta da probabilityislogic .

n_steps <- 1000
target_min <- 1.99
target_max <- 2.01
target_prob <- 0.99

target_mean <- mean(c(target_min,target_max))
target_sd <- (target_max-target_mean)/qnorm(p=1-(1-target_prob)/2)

mm <- rep(target_mean/n_steps,n_steps)

# boring setting:
# sigma_sq <- target_sd^2/(2*n_steps)
# tau <- target_sd^2/(2*n_steps*(n_steps-1))

sigma_sq <- 1
tau <- (target_sd^2/n_steps-sigma_sq)/(n_steps-1)

CC <- matrix(tau,nrow=n_steps,ncol=n_steps)
diag(CC) <- sigma_sq

library(MASS)
foo <- mvrnorm(1,mu=mm,Sigma=CC)
sum(foo)

plot(cumsum(foo),type="l",xlab="",ylab="")
abline(h=target_mean,lty=2)