Simulação envolvendo condicionamento na soma de variáveis ​​aleatórias

8

Eu estava lendo esta pergunta e pensei em simular a quantidade necessária. O problema é o seguinte: Se e são normais normais, o que é ? Então, eu quero simular . (para um valor escolhido de )B E ( A 2 | A + B ) E ( A 2 | A + B ) A + BABE(A2|A+B)E(A2|A+B)A+B

Eu tentei o seguinte código para conseguir isso:

n <- 1000000
x <- 1 # the sum of A and B

A <- rnorm(n)
B <- rnorm(n)

sum_AB = A+B

estimate <- 1/sum(sum_AB==x) * sum( (A[sum_AB==x])^2 )

O problema é que quase sempre não há valor em sum_ABque correspondências x(entre simulações). Se eu escolher algum elemento sum_AB, geralmente é a única instância de seu valor no vetor.

Em geral, como lidar com esse problema e executar uma simulação precisa para encontrar uma expectativa do formulário fornecido? ( e podem não ser necessariamente normalmente distribuídos ou da mesma distribuição.)BAB

Comp_Warrior
fonte
1
Sua edição recente altera substancialmente a questão, como indica nosso intercâmbio de comentários. Torna-se mais difícil responder na generalidade muito maior que você supõe agora. Por exemplo, existem técnicas especiais - e bastante envolvidas - apenas para responder quando o valor de é raro (em uma das caudas). A+B
whuber
@whuber Todos os valores não seriam relativamente raros quando lidamos com duas variáveis ​​aleatórias contínuas?
Comp_Warrior 25/08
1
Sim, mas faixas estreitas de valores - que geralmente são suficientes para tais simulações - nunca funcionariam nos caudas (nem em nenhuma outra região onde o PDF fica muito pequeno), enquanto que quando a densidade é relativamente grande, você pode executar facilmente um teste. cálculo de força bruta que é garantido para produzir um número decente de dados com perto o suficiente do valor desejado para permitir que algumas conclusões sejam tiradas da simulação. A+B
whuber
@ whuber entendo - você poderia dar alguma indicação em sua resposta das técnicas especiais mencionadas? Desculpas por não indicar o que eu estava interessado abaixo nos comentários.
Comp_Warrior
Comp_Warrior Estou anexando uma segunda solução que acredito ser a que @whuber está se referindo.
Dan

Respostas:

5

Meu comentário no tópico mencionado sugere uma abordagem eficiente: porque e são conjuntamente normais com covariância zero, eles são independentes, de onde a simulação precisa apenas gerar (que tem média e variação ) e construa . Neste exemplo, a distribuição de é examinada por meio do histograma de valores simulados.X=A+BY=ABY02A=(X+Y)/2A2|(A+B=3)105

x <- 3
y <- rnorm(1e5, 0, sqrt(2))
a <- (x+y)/2
hist(a^2)

A expectativa pode ser estimada como

mean(a^2)

A resposta deve estar próxima de .11/4=2.75

whuber
fonte
1
Obrigado - isso faz sentido. No entanto, estou certo ao entender que essa simplificação só funcionará se as duas variáveis ​​aleatórias em questão forem normais? E se eu tivesse um caso em que e fossem de outra distribuição (e possível separar um do outro)? AB
Comp_Warrior
1
Sua compreensão está correta. Esta é uma das razões pelas quais as variáveis ​​normais são tão populares, tanto na teoria quanto nos modelos de computador! No entanto, a idéia básica de procurar uma maneira de transformar as variáveis ​​em conjuntos de variáveis ​​independentes (ou facilmente relacionadas) será transferida para um cenário mais geral.
whuber
2

Uma maneira genérica de resolver esse problema é considerar a mudança de variáveis ​​de para . O Jacobiana deste transformar ser igual a um (1), a densidade de é Portanto, a densidade de condição de é com o termo de proporcionalidade sendo o inverso da densidade marginal de , . Como , uma transformação determinística, essa também é a densidade da junta de dada(A,B)(A,A+B=S)(A,S)

fA,S(a,s)=fA(a)fB(sa)
AS=s
fA|S(a|s)fA(a)fB(sa)
SfS(s)1B=SAS f A , B | S ( a , b | s ) f A ( a ) f B ( s - a ) I a + b = s(A,B)S
fA,B|S(a,b|s)fA(a)fB(sa)Ia+b=s
A geração de uma realização a partir desse destino pode ser feita diretamente se a forma for simples o suficiente ou por aceitar-rejeitar, Metropolis-Hastings, amostragem de fatia ou qualquer outro método de simulação padrão.
Xi'an
fonte
1

Você pode resolver esse problema usando amostras de inicialização. Por exemplo,

n <- 1000000

A <- rnorm(n)
B <- rnorm(n)
AB <- cbind(A,B)

boots <- 100
bootstrap_data <- matrix(NA,nrow=boots*n,ncol=2)


for(i in 1:boots){
    index <- sample(1:n,n,replace=TRUE)
    bootstrap_data[(i*n-n+1):(i*n),] <- cbind(A[index],B[index]) 
}

sum_AB <- bootstrap_data[,1] + bootstrap_data[,2]
x <- sum_AB[sample(1:n,1)]

idx <- which(sum_AB == x)

estimate <- mean(bootstrap_data[idx,1]^2)

Executando este código, por exemplo, obtenho o seguinte

> estimate
[1] 0.7336328
> x
[1] 0.9890429

Portanto, quando então .A+B=0.9890429E(A2|A+B=0.9890429)=0.7336328

Agora, para validar que essa deve ser a resposta, vamos executar o código do whuber em sua solução. Então, executando o código com x<-0.9890429resultados da seguinte maneira:

> x <- 0.9890429
> y <- rnorm(1e5, 0, sqrt(2))
> a <- (x+y)/2
> hist(a^2)
>
> mean(a^2)
[1] 0.745045

E, portanto, as duas soluções são muito próximas e coincidem uma com a outra. No entanto, minha abordagem para o problema deve permitir que você insira qualquer distribuição que desejar, em vez de confiar no fato de que os dados vieram de distribuições normais.


Uma segunda solução de força bruta, mais baseada no fato de que, quando a densidade é relativamente grande, você pode executar facilmente um cálculo de força bruta, é a seguinte

n <- 1000000

x <- 3  #The desired sum to condition on

A <- rnorm(n)
B <- rnorm(n)
sum_AB <- A+B

epsilon <- .01
idx <- which(sum_AB > x-epsilon & sum_AB < x+epsilon)
estimate <- mean(A[idx]^2)

estimate

Executando este código, obtemos o seguinte

> estimate
[1] 2.757067

Assim, a execução do código para resulta em que concorda com a verdadeira solução.A+B=3E(A2|A+B=3)=2.757067

Dan
fonte
1
Devo estar faltando alguma coisa: a questão pede ao usuário para especificar o valor de . Onde isso é feito no seu código? Como seria seu código no caso de precisar ser definido como , por exemplo? A+BA+B3
whuber
@whuber você está completamente correto. Só posso fazer isso por somas que sei que aparecerão.
Dan
0

parece-me que a questão se torna esta:

  1. como simular (X, Y) condicional em X + Y = ke depois
  2. use monte carlo para estimar UE (X, Y) para alguma função U (x, y)

vamos começar revisando a amostragem de importância :

EV(Z1)=V(z)f1(z)=V(z)f1(z)f2(z)f2(z)=EV(Z2)f1(Z2)f2(Z2)

onde os primeiros expectativas é em relação à variável aleatória com densidade e o segundo é um wrt com densidade .Z1f1(z)Z2f2(z)

Portanto, se você pode simular aleatoriamente 's de , faça uma estimativa usando ou alternativamente simule ' s de usandozif11niV(zi)zif21niV(zi)f1(zi)f2(zi)

Agora, voltemos ao nosso caso e são distribuídos como (X, Y) na condição X + Y = k, ou seja, e deixeU(x,y)=x2(X,Y)f(x,y)x+y=kf(x,y)A=x+y=kf(x,y)

então agora o procedimento é:

  1. gere cópias nid da densidade - e chame-as deg(x)Xi
  2. Seja observe que a distribuição deste (X, Y) é , onde é função indicadora g ( x ) I ( x + y = k ) I ( )Yi=kXig(x)I(x+y=k)I()
  3. a estimativa é
    1niU(xi,yi)f(xi,yi)Ag(xi)
pes
fonte
1
Sua solução não está correta desde . A=0
Xi'an