Suavização de dados 2D

Os dados consistem em espectros ópticos (intensidade da luz em relação à frequência) gravados em momentos variados. Os pontos foram adquiridos em uma grade regular em x (tempo), y (frequência). Para analisar a evolução do tempo em frequências específicas (um aumento rápido, seguido de uma deterioração exponencial), eu gostaria de remover parte do ruído presente nos dados. Esse ruído, para uma frequência fixa, provavelmente pode ser modelado como aleatório com distribuição gaussiana. Em um horário fixo, no entanto, os dados mostram um tipo diferente de ruído, com picos espúrios grandes e oscilações rápidas (+ ruído gaussiano aleatório). Tanto quanto posso imaginar, o ruído ao longo dos dois eixos não deve ser correlacionado, pois tem origens físicas diferentes.

Qual seria um procedimento razoável para suavizar os dados? O objetivo não é distorcer os dados, mas remover artefatos ruidosos "óbvios". (e a suavização excessiva pode ser ajustada / quantificada?) Não sei se a suavização em uma direção independentemente da outra faz sentido ou se é melhor suavizar em 2D.

Eu li coisas sobre estimativa de densidade de kernel 2D, interpolação de polinômio / spline 2D, etc., mas não estou familiarizado com o jargão ou a teoria estatística subjacente.

Eu uso R, para o qual vejo muitos pacotes que parecem relacionados (MASS (kde2), campos (smooth.2d), etc.), mas não consigo encontrar muitos conselhos sobre qual técnica aplicar aqui.

Fico feliz em saber mais, se você tiver referências específicas para me indicar (ouvi dizer que MASS seria um bom livro, mas talvez seja muito técnico para um não estatístico).

Edit: Aqui está um espectrograma fictício representativo dos dados, com fatias ao longo das dimensões de tempo e comprimento de onda.

O objetivo prático aqui é avaliar a taxa de decaimento exponencial no tempo para cada comprimento de onda (ou caixas, se muito barulhento).

Você precisa especificar um modelo que separa o sinal do ruído.

Existe o componente de ruído no nível de medição que você assume gaussiano. Os outros componentes, dependentes nas medições:

• "Esse ruído, para uma frequência fixa, provavelmente pode ser modelado como aleatório com distribuição gaussiana". Precisa de esclarecimentos - o componente de ruído é comum a todos os pontos do tempo, dada a frequência? O desvio padrão é o mesmo para todas as frequências? Etc.

• "Em um horário fixo, no entanto, os dados mostram um tipo diferente de ruído, com grandes picos espúrios e oscilações rápidas" Como você separa isso do sinal, pois provavelmente você está interessado na variação da intensidade através da frequência. A variação interessante é de alguma forma diferente da variação desinteressante? Em caso afirmativo, como?

Oscilações espúrias ou ruído não gaussiano em geral não são um grande problema, se você tiver uma idéia realista de suas características. Pode ser modelado transformando os dados (e depois usando um modelo gaussiano) ou explicitamente usando uma distribuição de erro não gaussiana. O ruído de modelagem correlacionado sobre as medições é mais desafiador.

Dependendo de qual é o seu modelo de ruído e dados, você pode modelar os dados com uma ferramenta de uso geral, como os GAMs no pacote mgcv, ou pode precisar de uma ferramenta mais flexível, que leve facilmente a uma configuração bayesiana bastante personalizada . Existem ferramentas para esses modelos, mas se você não for estatístico, aprender a usá-los levará um tempo.

Eu acho que ou uma solução específica para análise espectral ou o pacote mgcv são suas melhores apostas.

Uma série temporal de espectros me sugere um experimento de cinética , e há uma quantidade bem estabelecida de literatura quimiométrica sobre isso.

O que você sabe sobre os espectros? Que tipo de espectro são eles? Você pode razoavelmente esperar que você tenha apenas duas espécies, eduto e produto?

Você pode razoavelmente assumir a bilinearidade, ou seja, os espectros medidos em um determinado momento são uma combinação linear das concentrações do componente vezes o espectro puro do componente : $\mathbf X$$\mathbf C$$\mathbf S$

$\mathbf X^{(nspc \times nwl)} = \mathbf C^{(nspc \times ncomp)} \mathbf S^{(ncomp \times nwl)}$

Você diz que deseja estimar uma deterioração exponencial (nas concentrações). Isso, juntamente com a bilinearidade, sugere-me a resolução de curvas multivariadas (MCR). Essa é uma técnica que permite que você use as informações que você possui (por exemplo, espectros de componentes puros de algumas substâncias ou suposições sobre o comportamento da concentração como decaimento exponencial) durante o ajuste do modelo.

• Um bom ponto de partida na literatura é:
  Anna de Juan, Marcel Maeder, Manuel Martínez Romà Tauler: Combinando modelagem rígida e flexível para resolver problemas cinéticos, Chemometrics and Intelligent Laboratory Systems 54, 2000. 123–141.
• Um bom ponto de partida na internet é http://www.mcrals.info/
• O UseR! Livro de quimiometria discute isso no final
• Implementação no R: package ALS

Até onde eu sei, é bastante comum suavizar as concentrações de acordo com algum modelo, por exemplo, cinético, mas é muito menos comum suavizar os espectros. No entanto, o algoritmo permite fazer isso. No verão, perguntei a Anna se elas impõem restrições de suavidade, mas ela me disse que não (e os bons espectroscopistas odeiam a suavização em vez de medir bons espectros ;-)). Muitas vezes, isso também não é necessário, porque a agregação de informações de todos os espectros já produzirá boas estimativas dos espectros de componentes puros.

Eu suavizei "espectros de componentes" (de fato, componentes principais) duas vezes recentemente ( Dochow et al .: Dispositivo Raman no chip e fibras de detecção com grade de Bragg de fibra para análise de soluções e partículas, LabChip, 2013 e Dochow el al. : Chip microfluídico de quartzo para identificação de células tumorais por espectroscopia Raman em combinação com armadilhas ópticas, AnalBioanalChem, aceito), mas nesses casos meu conhecimento espectroscópico me dizia que eu posso fazer isso. Aplico regularmente uma interpolação de downsampling e suavização em meus espectros Raman ( hyperSpec::spc.loess).

Como saber o que é muita suavização? Eu acho que a única resposta possível é "conhecimento especializado sobre o tipo de espectroscopia e experimento".

editar: reli a pergunta e você diz que deseja estimar a deterioração em cada comprimento de onda. No entanto, isso é verdade ou você deseja estimar a deterioração de diferentes espécies com espectros sobrepostos?

Os dados consistem em espectros ópticos (intensidade da luz em relação à frequência) gravados em> tempos variados. Os pontos foram adquiridos em uma grade regular em x (tempo), y (frequência).

$intensity = f(time, frequency)$$f$sendo uma soma de funções básicas (por exemplo, splines-b) e coeficientes. Um conjunto limitado de funções básicas reduz diretamente a rugosidade e, portanto, cancela boa parte do ruído branco.

Eu li coisas sobre estimativa de densidade de kernel 2D, interpolação de polinômio / spline 2D, etc.

...

Eu uso R, para o qual vejo muitos pacotes que parecem relacionados (MASS (kde2), campos (smooth.2d), etc.), mas não consigo encontrar muitos conselhos sobre qual técnica aplicar aqui.

Você mencionou a interpolação de spline, mas não mencionou o pacote fda, que implementa de maneira bastante fácil e fácil a expansão da função básica que mencionei acima. O conjunto de medições simultâneas de tempo, frequência e intensidade (ordenadas como uma matriz tridimensional) pode ser capturado como um objeto de dados funcionais bivariados, ver. por exemplo, a função 'Data2fd'. Além disso, vários procedimentos de suavização estão disponíveis na embalagem, todos projetados para cancelar o ruído branco ou a "rugosidade" nas medições de processos inerentemente suaves.

O artigo da Wikipedia descreve o problema do ruído branco no FDA da seguinte maneira:

Os dados podem ser tão precisos que o erro pode ser ignorado, pode estar sujeito a um erro de medição substancial ou até ter um relacionamento indireto complexo com a curva que eles definem. ... registros diários de precipitação em uma estação meteorológica são tão variáveis ​​que exigem análises cuidadosas e sofisticadas para extrair algo como uma curva de precipitação média.

O FDA fornece as ferramentas para esses casos. Isso não se traduz no seu caso?

... mas não estou familiarizado com o jargão ou a teoria estatística subjacente ...

... mas não consigo encontrar muitos conselhos sobre qual técnica aplicar aqui ...

Com relação à fda: eu não era nenhum dos dois, mas o livro de Ramsay e Silverman sobre a FDA (2005) torna o básico muito bem acessível e Ramsay Hooker e Graves (2009) traduzem diretamente os insights do livro em código R. Ambos os volumes devem estar disponíveis como e-books em uma biblioteca universitária de estatística, biociências, climatologia ou psicologia. O Google também exibirá mais alguns links que não posso postar aqui.

Lamento não poder fornecer uma solução mais direta para o seu problema. No entanto, o FDA me ajudou muito depois que eu descobri para que serve.

Sendo um físico simples, não um especialista em estatística, eu adotaria uma abordagem simples. As duas dimensões são de naturezas diferentes. Faria sentido suavizar ao longo do tempo com um algoritmo e suavizar ao longo do comprimento de onda com outro.

Os algoritmos reais que eu usaria: para comprimento de onda, Savitzky-Golay com uma ordem superior, 6 talvez 8.

Com o tempo, se esse exemplo for típico, esse salto repentino e um declínio mais ou menos exponencial o tornam complicado. Eu tive dados experimentais e imagens barulhentas, exatamente assim. Se métodos simples e diretos não ajudarem o suficiente, tente um Gaussiano mais suave, mas suprima seu efeito próximo ao salto, como detectado por um detector de borda. Suavize e amplie a saída do detector de bordas, normalize-a para ir de 0,0 a 1,0 e use-a para selecionar entre a imagem original e a imagem suavizada de Gauss, pixel por pixel.

@ batiste: Fico feliz que você adicionou o enredo como sugeri. Isso ajuda muito:

Portanto, se eu entendi direito, seu objetivo prático é avaliar a taxa de decaimento exponencial para cada comprimento de onda; então vamos fazer exatamente isso! Defina uma função que você deseja minimizar para cada comprimento de onda separadamente e minimize-a.

Vejamos um único comprimento de onda, como no gráfico inferior direito.

$\tau$

$\tau$$\tau$

Posteriormente, se você acredita que o comprimento de onda adjacente deve ter constantes de decaimento semelhantes, você pode incorporar isso a um critério de otimização mais elaborado.

Se alguma coisa, eu sugiro que você leia um livro de otimização que deve ser lido: a otimização convexa de Boyd .

Espero que isto ajude!