Interpretação probabilística de splines de suavização de placas finas

8

TLDR: As splines de regressão em placas finas têm uma interpretação probabilística / bayesiana?

Dado pares de entrada-saída (xi,yi) , i=1,...,n ; Eu quero estimar uma função f() seguinte forma

f(x)u(x)=ϕ(xi)Tβ+i=1nαik(x,xi),
k(,)ϕ(xi)m<nαiβi
minαRn,βRm1nYΦβKαRn2+λαTKα,
Φϕ(xi)Ti,jKk(xi,xj)
α=λ1(I+λ1K)1(YΦβ)
β={ΦT(I+λ1K)1Φ}1ΦT(I+λ1K)1Y.
Supondo que k(,) seja uma função definida positiva do kernel, esta solução pode ser vista como o melhor preditor linear e imparcial do seguinte modelo bayesiano: h ( ) G P ( 0 , τ k ( , ) ) , β 1 ,
y | (β,h())  N(ϕ(x)β+h(x),σ2),
h()  GP(0,τk(,)),
β1,
que σ2/τ=λ e GP denota um processo gaussiano. Veja por exemplohttps://www.ncbi.nlm.nih.gov/pmc/articles/PMC2665800/

Minha pergunta é a seguinte. Suponha que eu deixe k(x,x):=|xx|2ln(|xx|) e ϕ(x)T=(1,x) , isto é, spline de placa fina regressão. Agora, k(,) não é uma função semidefinida positiva e a interpretação acima não funciona. O modelo acima e sua solução ainda têm uma interpretação probabilística, como no caso em que k(,) é semidefinido positivo?

MthQ
fonte
Você parece assumir o está em uma espaço dimensional com ou pelo menos é o inteiro é mesmo. d d = 2 dxdd=2d
Yves
Ok, então quais são as implicações?
MthQ 11/0318
2
Este foi apenas um comentário secundário, porque na questão pode-se pensar que são escalares. Mas, neste caso, o kernel de Duchon tem a forma com inteiro, e para o spline de suavização usual. Penso que a interpretação probabilística permanece quase inalterada, mas o GP não é estacionário: é uma Função Aleatória Intrínseca . Para o spline de suavização usual, esse acaba sendo um processo Wiener integrado. | x - x | 2 m - 1 m m = 2xi|xx|2m1mm=2
Yves
1
@ Yves isso parece interessante. Você pode expandir seu comentário para uma resposta, explicando um pouco mais o que é uma função aleatória intrínseca e adicionando o exemplo clássico do spline de suavização. Se você se preocupa em provar que o kernel do TPS gera um GP não estacionário, talvez uma simulação possa ser um compromisso útil, especialmente se você adicionar uma estimativa não paramétrica da variação da distribuição preditiva posterior.
DeltaIV
@DeltaIV. Obrigado. Vou tentar fazê-lo, ainda não é uma tarefa fácil. Estou certo de que isso quando as funções são polinômios adequados relacionados ao kernel, mas isso pode não ser mais verdade com arbitrário como no contexto mais clássico do GP. ϕ jϕjϕj
Yves

Respostas:

5

Deixe o modelo da pergunta ser escrito como que é um GP não observado com índice e é um termo de ruído normal com variância . Presume-se que o GP seja centrado, estacionário e não determinístico. Observe que o termo pode ser considerado como um GP (determinístico) com o kernel que h(x)xRdεiσ2φ(x)pφ(x)B

(1)Yi=ϕ(xi)β+h(xi)+εi
h(x)xRdεiσ2ϕ(x)βB B : = ρϕ(x)Bϕ(x)Bé uma matriz de covariância de `` valor infinito ''. De fato, tomando com , obtemos as equações de krigagem da pergunta. Isso geralmente é chamado de difuso anterior para . Um posterior adequado para resulta apenas quando a matriz possui classificação completa. Portanto, o modelo escreve e que é um GP . A mesma interpretação de Bayes pode ser usada com restrições quando não é mais um GP, mas é um ρ β β Φ Y i = ζ ( x i ) + ε i ζ ( x ) ζ ( x )B:=ρIρββΦ
(2)Yi=ζ(xi)+εi
ζ(x)ζ(x)Função aleatória intrínseca (IRF). A derivação pode ser encontrada no livro de G. Wahba. Apresentações legíveis do conceito de IRF estão, por exemplo, no livro de N. Cressie e no artigo de Mardia et al. Os IRFs são semelhantes aos conhecidos processos integrados no contexto de tempo discreto (como o ARIMA): um IRF é transformado em um GP clássico por um tipo de operação diferenciada.

Aqui estão dois exemplos de IRF para . Primeiramente, considere um processo de Wiener com sua condição inicial substituída por uma condição inicial difusa : é normal com uma variação infinita. Depois que um valor é conhecido, o IRF pode ser previsto como o Wiener GP. Em segundo lugar, considere um processo integrado de Wiener fornecido pela equação em que é um processo de Wiener. Para obter um GP, agora precisamos de dois parâmetros escalares: dois valores e paraζ ( x ) ζ ( 0 ) = 0 ζ ( 0 ) ζ ( x ) d 2 ζ ( x ) / d x 2 = d W ( x ) / d x W ( x ) ζ ( x ) ζ ( x ) x x ζ ( x )d=1ζ(x)ζ(0)=0ζ(0)ζ(x)

d2ζ(x)/dx2=dW(x)/dx
W(x)ζ(x)ζ(x)xxou os valores e em algum escolhido . Podemos considerar que os dois parâmetros extras são em conjunto gaussianos com uma matriz de covariância infinita . Em ambos os exemplos, assim que um conjunto finito de observações estiver disponível, o IRF é quase como um GP. Além disso utilizou-se um operador diferencial: e , respectivamente. O espaço nulo é um espaço linear de funções tal que . Ele contém a função constante ζ(x)dζ(x)/dxx2×2L:=d/dxL:=d2/dx2Fϕ(x)Lϕ=0ϕ1(x)=1no primeiro caso e as funções e no segundo caso. Observe que no primeiro exemplo é GP para qualquer fixo no primeiro exemplo e da mesma forma é um GP no segundo caso.ϕ1(x)=1ϕ2(x)=xζ(x)ζ(x+δ)δζ(xδ)2ζ(x)+ζ(x+δ)

Para uma dimensão geral , considere um espaço linear de funções definidas em . Chamamos um incremento relativo a um conjunto finito de locais e pesos reais tal que Considere como o espaço nulo de nossos exemplos. Para o primeiro exemplo, podemos usar, por exemplo, com e arbitrários edFRdFsxiRdsνi

i=1sνiϕ(xi)=0 for all ϕF.
Fs=2x1x2[1,1] . Para o segundo exemplo, podemos considerar igualmente espaçados e . A definição de um IRF envolve um espaço de funções e uma função que é condicionalmente positiva wrt , o que significa que permanece assim que é um incremento wrt . De es=3xiν=[1,2,1]Fg(x,x)F
i=1sj=1sνiνjg(xi,xj)0
[νi,xi]i=1sFFg(x,x) podemos fazer um núcleo de covariância, portanto, um GP, como em Mardia et al. Podemos começar a partir de um operador diferencial linear e usar o espaço nulo como ; o IRF terá então conexão com a equação um ruído gaussiano.LFLζ=

O cálculo da previsão do IRF é quase o mesmo da pergunta, com substituído por , mas com o agora formando uma base de . A restrição extra deve ser adicionada no problema de otimização, que concederá esse . Ainda podemos adicionar mais funções básicas que não estão em se necessário; isso terá o efeito de adicionar um GP determinístico, digamos ao IRF k(x,x)g(x,x)ϕi(x)FΦα=0αKα0Fψ(x)γζ(x) em (2).

O spline de placa fina depende de um número inteiro tal que , o espaço contenha polinômios de baixo grau, com a dimensão dependendo de e . Pode-se mostrar que se é a seguinte função para depois define um wrt condicionalmente positivo . A construção refere-se a um operador diferencialmm>2dFp(m)mdE(r)r0

E(r):={(1)m+1+d/2r2mdlogrd even,r2mdd odd,
g(x,x):=E(xx)FL. Acontece que, para e a ranhura fina de chapa é nada além da ranhura cúbica natural usual, que se refere ao exemplo integrado de Wiener acima, com . Portanto, (2) nada mais é do que o modelo de spline de suavização usual. Quando e o espaço nulo tem dimensão e é gerado pelas funções , e .d=1m=2g(x,x)=|xx|3d=2m=2p(m)=31x1x2

Estatísticas de Cressie N para dados espaciais . Wiley 1993.

Mardia KV, Kent JT, Goodall CR e Little JA. Krigagem e splines com informações derivadas. Biometrika (1996), 83,1, pp. 207-221.

Modelos Wahba G Spline para dados observacionais . SIAM 1990.

Wang, Y Suavizando splines, métodos e aplicações . Chapman e Hall, 2011.

Yves
fonte
Muito obrigado pelo esforço realizado. Extremamente útil. Eu tenho uma pergunta adicional. Portanto, adicionar funções adicionais a (além das funções de ) não altera a interpretação de . O que notei, no entanto, é que a solução dada na minha pergunta acima sempre satisfaz , não apenas se . Como isso pode ser interpretado? F ζ ( ) α Φ α = 0 ϕ ( ) Fϕ()Fζ()αΦα=0ϕ()F
MthQ 19/0318
Sim. Nos dois casos, existem funções básicas na aproximação de , enquanto apenas observações são usadas. Portanto, temos algo como uma regressão com classificação deficiente com coeficientes e . Como a parte não é penalizada, tende a "absorver" mais da variação de do que a parte que traz restrições lineares. Note que nada proíbe o uso de algumas das "mudanças do kernel" funções como . Se usarmos todos eles, todosf ( x ) n β i α j β y α p n x k ( x , x i ) ϕ j ( x ) α jn+pf(x)nβiαjβyαpnxk(x,xi)ϕj(x)αjsão zero, o que parece sensato.
Yves