Como usar corretamente o código GPML Matlab para um problema real (não demo)?

Fiz o download do código GPML Matlab mais recente e o código GPML Matlab e li a documentação e executei a demonstração de regressão sem problemas. No entanto, estou tendo dificuldade em entender como aplicá-lo a um problema de regressão com o qual me deparo.

O problema de regressão é definido da seguinte maneira:

Seja um vetor de entrada e seja seu destino correspondente. O conjunto de entradas é organizado em uma matriz e seus destinos correspondentes são armazenados em uma matriz , com sendo o valor alvo médio em .$\mathbf{x}_i \in \mathbb{R}^{20}$$\mathbf{y}_i \in \mathbb{R}^{25}$$M$$\mathbf{X} = [\mathbf{x}_1, \dots, \mathbf{x}_M]^\top$$\mathbf{Y} = [\mathbf{y}_1 - \mathbf{\bar{y}}, \dots, \mathbf{y}_M-\mathbf{\bar{y}}]^\top$$\mathbf{\bar{y}}$$\mathbf{Y}$

Desejo treinar um modelo GPR usando a função exponencial ao quadrado:$\mathcal{G} = \lbrace \mathbf{X}, \mathbf{Y}, \theta \rbrace$

$k(\mathbf{x}_i, \mathbf{x}_j) = \alpha^2 \text{exp} \left( - \frac{1}{2\beta^2}(\mathbf{x}_i - \mathbf{x}_j)^2\right) + \gamma^2\delta_{ij}$ ,

onde é igual a se e caso contrário. Os hiperparâmetros são com sendo o nível de ruído assumido nos dados de treinamento e é a escala de comprimento.$\delta_{ij}$$1$$i = j$$0$$\theta = (\alpha, \beta, \gamma)$$\gamma$$\beta$

Para treinar o modelo, preciso minimizar a probabilidade marginal de log negativa em relação aos hiperparâmetros:

$-\text{log}\, p(\mathbf{Y} \mid \mathbf{X}, \theta) = \frac{1}{2} \text{tr}(\mathbf{Y}^\top\mathbf{K}^{-1}\mathbf{Y}) + \frac{1}{2}\text{log}\mid\mathbf{K}\mid + \,c,$

onde c é uma constante e a matriz é uma função dos hiperparâmetros (veja a equação k (xi, xj) = ...).$\mathbf{K}$

Com base na demonstração do site GPML, minha tentativa de implementar isso usando o código GPML Matlab está abaixo.

covfunc = @covSEiso;
likfunc = @likGauss;
sn = 0.1;
hyp.lik = log(sn);
hyp2.cov = [0;0];
hyp2.lik = log(0.1);
hyp2 = minimize(hyp2, @gp, -100, @infExact, [], covfunc, likfunc, X1, Y1(:, n));
exp(hyp2.lik)
nlml2 = gp(hyp2, @infExact, [], covfunc, likfunc, X1, Y1(:, n));
[m s2] = gp(hyp2, @infExact, [], covfunc, likfunc, X1, Y1(:, n), X2);
Y2r(:, n) = m;

X1 contém as entradas de treinamento

X2 contém as entradas de teste

Y1 contém as metas de treinamento

Y2r são as estimativas da aplicação do modelo

n é o índice usado para regredir cada elemento no vetor de saída

Dado o problema, esta é a maneira correta de treinar e aplicar o modelo GPR? Caso contrário, o que preciso alterar?

O GP faz um bom trabalho para os dados de treinamento do seu problema. No entanto, não é tão bom para os dados de teste. Você provavelmente já executou algo como o seguinte:

load('../XYdata_01_01_ab.mat');

for N = 1 : 25
    % normalize
    m = mean(Y1(N,:));
    s = std(Y1(N,:));
    Y1(N,:) = 1/s * (Y1(N,:) - m);
    Y2(N,:) = 1/s * (Y2(N,:) - m);

    covfunc = @covSEiso;
    ell = 2;
    sf = 1;
    hyp.cov = [ log(ell); log(sf)];

    likfunc = @likGauss;
    sn = 1;
    hyp.lik = log(sn);

    hyp = minimize(hyp, @gp, -100, @infExact, [], covfunc, likfunc, X1', Y1(N,:)');
    [m s2] = gp(hyp, @infExact, [], covfunc, likfunc, X1', Y1(N,:)', X1');    
    figure;    
    subplot(2,1,1); hold on;    
    title(['N = ' num2str(N)]);    
    f = [m+2*sqrt(s2); flipdim(m-2*sqrt(s2),1)];
    x = [1:length(m)];
    fill([x'; flipdim(x',1)], f, [7 7 7]/8);
    plot(Y1(N,:)', 'b');
    plot(m, 'r');
    mse_train = mse(Y1(N,:)' - m);

    [m s2] = gp(hyp, @infExact, [], covfunc, likfunc, X1', Y1(N,:)', X2');
    subplot(2,1,2); hold on;
    f = [m+2*sqrt(s2); flipdim(m-2*sqrt(s2),1)];
    x = [1:length(m)];
    fill([x'; flipdim(x',1)], f, [7 7 7]/8);    
    plot(Y2(N,:)', 'b');
    plot(m, 'r');
    mse_test = mse(Y2(N,:)' - m);

    disp(sprintf('N = %d -- train = %5.2f   test = %5.2f', N, mse_train, mse_test));
end

Ajustando os hiperparâmetros manualmente e não usando a função de minimização, é possível equilibrar o trem e testar o erro, mas ajustar o método observando o erro de teste não é o que você deve fazer. Eu acho que o que está acontecendo é uma super adaptação aos seus três sujeitos que geraram os dados do treinamento. Nenhum método imediato fará um bom trabalho aqui, e como poderia? Você fornece os dados de treinamento, para que o método tente obter o melhor possível sobre os dados de treinamento sem superajustar. E, de fato, não superajustar no sentido clássico. Ele não se adapta aos dados, mas se adapta aos três sujeitos do treinamento. Por exemplo, a validação cruzada com o conjunto de treinamento nos diria que não há ajuste excessivo. Ainda assim, seu conjunto de testes será mal explicado.

O que você pode fazer é:

1. Obtenha dados de mais assuntos para treinamento. Dessa forma, sua quarta pessoa terá menos probabilidade de parecer uma "discrepância", como atualmente. Além disso, você tem apenas uma sequência de cada pessoa, certo? Talvez ajude a gravar a sequência várias vezes.

2. De alguma forma, incorpore conhecimentos prévios sobre sua tarefa que impediriam que um método se ajustasse a assuntos específicos. Em um GP, isso pode ser feito através da função de covariância, mas provavelmente não é tão fácil de fazer ...

3. Se não me engano, as sequências são de fato séries temporais. Talvez faça sentido explorar as relações temporais, por exemplo, usando redes neurais recorrentes.

Definitivamente, há mais, mas essas são as coisas que consigo pensar agora.

Eu acho que o problema pode ser um erro de especificação do modelo. Se seus alvos são ângulos ajustados a + -180 graus, o "processo de ruído" para seus dados pode ser suficientemente não-guassiano, para que a evidência baysiana não seja uma boa maneira de otimizar os hiperparâmetros. Por exemplo, considere o que acontece quando o "ruído" causa a quebra do sinal. Nesse caso, pode ser aconselhável executar a seleção do modelo, minimizando o erro de validação cruzada (há uma implementação de domínio público do método simplex Nelder-Mead aquise você não tiver a caixa de ferramentas de otimização). A estimativa de desempenho com validação cruzada não é tão sensível à especificação incorreta do modelo, pois é uma estimativa direta do desempenho do teste, enquanto a probabilidade marginal do modelo é a evidência no suporte ao modelo, uma vez que as premissas do modelo estão corretas. Veja a discussão que começa na página 123 do livro de Rasmussen e Williams.

Outra abordagem seria re-codificar as saídas para que um modelo de ruído gaussiano seja mais apropriado. Uma coisa que você pode fazer é alguma forma de redução de dimensionalidade não supervisionada, pois existem relações não lineares entre seus alvos (como há apenas uma maneira limitada pela qual um corpo pode se mover), haverá uma variedade de dimensões mais baixas que seu os alvos permanecem vivos e seria melhor regredir as coordenadas dessa variedade em vez dos próprios ângulos (também pode haver menos alvos dessa maneira).

Além disso, algum tipo de análise de Procrustes pode ser uma boa idéia para normalizar as diferenças entre os sujeitos antes de treinar o modelo.

Você pode encontrar parte do trabalho realizado por Neil Lawrence em recuperação de interesse humano. Lembro-me de ver uma demonstração disso em uma conferência há alguns anos e fiquei muito impressionado.