Simulando o posterior de um processo gaussiano

8

Pela primeira vez (desculpe imprecisões / erros) , observei os processos gaussianos e, mais especificamente, assisti a este vídeo de Nando de Freitas . As notas estão disponíveis online aqui .

Em algum momento, ele extrai 10 amostras aleatórias de um normal multivariado gerado pela construção de uma matriz de covariância baseada em um núcleo gaussiano (exponencial de distâncias quadradas no eixo x ). Essas amostras aleatórias formam os gráficos suaves anteriores que se tornam menos dispersos à medida que os dados se tornam disponíveis. Por fim, o objetivo é prever, modificando a matriz de covariância e obtendo a distribuição Gaussiana condicional nos pontos de interesse.

insira a descrição da imagem aqui

O código inteiro está disponível em um excelente resumo de Katherine Bailey aqui , que por sua vez credita um repositório de código de Nando de Freitas aqui . Eu publiquei o código Python aqui por conveniência.

Começa com (em vez de 10 acima) funções anteriores e introduz um "parâmetro de ajuste".310

Eu traduzi o código para Python e [R] , incluindo os gráficos:

Aqui está o primeiro pedaço de código em [R] e o gráfico resultante de três curvas aleatórias geradas por um kernel Gaussiano com base na proximidade dos valores no conjunto de testes:x

! [insira a descrição da imagem aqui

O segundo pedaço de código R é mais cabeludo e começa simulando quatro pontos de dados de treinamento, o que eventualmente ajudará a diminuir a propagação entre as possíveis curvas (anteriores) em torno das áreas onde esses pontos de dados de treinamento estão. A simulação do valor para esses pontos de dados é como uma função sin ( ) . Podemos ver o "aperto das curvas em torno dos pontos":ypecado()

insira a descrição da imagem aqui

O terceiro pedaço do código R trata da plotagem da curva dos valores médios estimados (o equivalente da curva de regressão), correspondendo aos valores de μ (veja o cálculo abaixo) e seus intervalos de confiança:50. μ

insira a descrição da imagem aqui

PERGUNTA: Quero pedir uma explicação das operações que ocorrem quando se passa do GP anterior para o posterior.

Especificamente, eu gostaria de entender esta parte do código R (no segundo bloco) para obter os meios e o sd:

# Apply the kernel function to our training points (5 points):

K_train = kernel(Xtrain, Xtrain, param)                          #[5 x 5] matrix

Ch_train = chol(K_train + 0.00005 * diag(length(Xtrain)))        #[5 x 5] matrix

# Compute the mean at our test points:

K_trte = kernel(Xtrain, Xtest, param)                            #[5 x 50] matrix
core = solve(Ch_train) %*% K_trte                                #[5 x 50] matrix
temp = solve(Ch_train) %*% ytrain                                #[5 x 1] matrix
mu = t(core) %*% temp                                            #[50 x 1] matrix

umaumaK_trainΣumaumaCh_traineuumaumaumaeK_trteΣumaeμ^50.

(Eq.1)μ^=[euumauma-1[5×5]Σumae[5×50.]]Teuumauma-1[5×5]ytr[5×1]dimensões=[50.×1]
# Compute the standard deviation:

tempor = colSums(core^2)                                          #[50 x 1] matrix

# Notice that all.equal(diag(t(core) %*% core), colSums(core^2)) TRUE

s2 = diag(K_test) - tempor                                        #[50 x 1] matrix
stdv = sqrt(s2)                                                   #[50 x 1] matrix

(Eq.2)var^=diag(Σee)-diag[[euumauma-1[5×5]Σumae[5×50.]]T[euumauma-1[5×5]Σumae[5×50.]]]=d[1111]-d[[euumauma-1[5×5]Σumae[5×50.]]T[euumauma-1[5×5]Σumae[5×50.]]]dimensões=[50.×1]

Como é que isso funciona?

μ^

Ch_post_gener = chol(K_test + 1e-6 * diag(n) - (t(core) %*% core))
m_prime = matrix(rnorm(n * 3), ncol = 3)
sam = Ch_post_gener %*% m_prime
f_post = as.vector(mu) + sam 

(Eq.3)fpostar=μ^+[euee[50.×50.]-[[euumauma-1[5×5]Σumae[5×50.]]T[euumauma-1[5×5]Σumae[5×50.]]]][N(0 0,1)[50.×3]]dimensões=[50.×3]
Antoni Parellada
fonte
No último gráfico, os intervalos de confiança não devem "beliscar" nos pontos conhecidos?
GeoMatt22
@ GeoMatt22 Eles meio que fazem, você não acha?
Antoni Parellada

Respostas:

2

eumaumae

[umae]N([μumaμe],[ΣumaumaΣumaeΣumaeTΣee])

E(x1|x2)=μ1+Σ12Σ22-1(x2-μ2)[50.×50.]Σumauma[50.×5]Σumae, uma transposição será necessária para tornar as matrizes congruentes em:

E(e|uma)=μe+ΣumaeTΣumauma-1(y-μuma)
μuma=μe=0 0

E(e|uma)=ΣumaeTΣumauma-1ytr

Entre na decomposição de Cholesky (que novamente codificarei em laranja como no OP):

E(e|uma)=ΣumaeTΣumauma-1ytr<--α-->=ΣumaeT(euumaumaeuumaumaT)-1ytr=ΣumaeTeuumauma-Teuumauma-1ytr(*)=ΣumaeTeuumauma-Teuumauma-1ytr<-m->

m=euumauma-1ytreuumaumam=ytrm


insira a descrição da imagem aqui


BTUMAT=(UMAB)T

μ^=[euumauma-1[5×5]Σumae[5×50.]]Teuumauma-1[5×5]ytr[5×1]=(ΣumaeTeuumauma-T)(euumauma-1ytr)dimensões=[50.×1]

dado que

[euumauma-1[5×5]Σumae[5×50.]]T=Σumae[50.×5]Teuumauma-1T[5×5]

Um raciocínio semelhante seria aplicado à variação, começando com a fórmula da variação condicional em um gaussiano multivariado:

vumar(x1|x2)=Σ11-Σ12Σ22-1Σ21

que no nosso caso seria:

varμ^e=Σee-ΣumaeTΣumauma-1Σumae=Σee-ΣumaeT[euumaumaeuumaumaT]-1Σumae=Σee-ΣumaeT[euumauma-1]Teuumauma-1Σumae=Σee-[euumauma-1Σumae]Teuumauma-1Σumae

e chegando na Eq. (2):

varμ^e=d[Kee-[euumauma-1[5×5]Σumae[5×50.]]T[euumauma-1[5×5]Σumae[5×50.]]]dimensões=[50.×1]

Podemos ver que a Eq. (3) no OP é uma maneira de gerar curvas aleatórias posteriores condicionais aos dados (conjunto de treinamento) e utilizar um formulário de Cholesky para gerar três sorteios aleatórios normais multivariados :

fpostar=μ^+[varμ^e][rnorm(0 0,1)]=μ^+[euee[50.×50.]-[[euumauma-1[5×5]Σumae[5×50.]]T[euumauma-1[5×5]Σumae[5×50.]]]][rand.norm's[50.×3]]dimensões=[50.×3]
Antoni Parellada
fonte
1
Isso é de um livro ou papel? Você tem uma maneira robusta de calcular média e variação condicionais quando a matriz de covariância está EXTREMAMENTE mal condicionada (mas sem excluir ou mesclar pontos de dados quase dependentes (próximos)) em dupla precisão? A precisão múltipla em software funciona, mas possui uma desaceleração de magnitude de 2,5 a 3 vezes em comparação com a precisão dupla de hardware; portanto, mesmo um algoritmo de precisão dupla "lento" será bom. Eu não acho que Cholesky corta isso. Também não acho que o QR o faça quando a matriz de covariância está muito mal condicionada. Usando backsolves padrão, parece precisar de precisão de ocutuple.
Mark L. Stone