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

7

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.

Stephan Kolassa
fonte
2
veja também a literatura sobre pontes brownianas ?
precisa
@BenBolker: parece uma ideia muito melhor que a minha. Você estaria interessado em escrever uma resposta?
Stephan Kolassa
talvez eu consiga, mas qualquer pessoa que esteja lendo isso deve se sentir à vontade para entrar e escrever uma resposta. Eu não me importo.
Ben Bolker
11
(Desculpe, é claro, quero dizer "BB-BB".) #
Stephan Kolassa

Respostas:

7

Geraremos normais multivariados com e forma que sua soma seja satisfatória nossa condição. Seja .XMN(μ,Σ)μRnΣRn×nZ=X1++Xn

Como um meio comum, escolhemos

μ1==μn=a+b2n.

Para que com probabilidade , seu desvio padrão deva cumprirZ[a,b]p

σZ=baqα,

onde é o quantil normal padrão para o nível , aqui .qααα=11p2

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ΣXiσ2cov(Xi,Xj)=τijΣXi

nσ2+n(n1)τ

então precisamos disso

nσ2+n(n1)τ=baqα.

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ΣΣ

σ2:=σZ22n,τ:=σZ22n(n1),

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

entediante

Menos chato é definir

σ2:=1,τ:=1n1(σZ2nσ2),

que produz trajetórias muito mais interessantes:

interessante

Observe que definir isso de fato produz uma matriz de covariância válida, porque é então da forma , a saberΣΣij=m(ij)

m(0)=σ2,m(j)=τ for j>0,

e nós temos isso

j>0|m(j)|=(n1)|τ|=|σZ2nσ2|=|σZ2n1|<1=σ2=m(0),

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

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)
Stephan Kolassa
fonte