Suavizar datos 2D

Los datos consisten en espectros ópticos (intensidad de luz contra frecuencia) registrados en diferentes tiempos. Los puntos se adquirieron en una cuadrícula regular en x (tiempo), y (frecuencia). Para analizar la evolución del tiempo a frecuencias específicas (un aumento rápido, seguido de una disminución exponencial), me gustaría eliminar parte del ruido presente en los datos. Este ruido, para una frecuencia fija, probablemente se puede modelar como aleatorio con distribución gaussiana. Sin embargo, en un momento fijo, los datos muestran un tipo diferente de ruido, con picos espurios grandes y oscilaciones rápidas (+ ruido gaussiano aleatorio). Por lo que puedo imaginar, el ruido a lo largo de los dos ejes no debería estar correlacionado ya que tiene diferentes orígenes físicos.

¿Cuál sería un procedimiento razonable para suavizar los datos? El objetivo no es distorsionar los datos, sino eliminar los artefactos ruidosos "obvios". (¿y se puede ajustar / cuantificar el exceso de suavizado?) No sé si tiene sentido suavizar a lo largo de una dirección independientemente de la otra, o si es mejor suavizar en 2D.

He leído cosas sobre la estimación de la densidad del kernel en 2D, la interpolación polinómica / spline en 2D, etc., pero no estoy familiarizado con la jerga o la teoría estadística subyacente.

Uso R, para lo cual veo muchos paquetes que parecen relacionados (MASS (kde2), campos (smooth.2d), etc.) pero no puedo encontrar muchos consejos sobre qué técnica aplicar aquí.

Me alegra saber más, si tiene referencias específicas para señalarme (escuché que MASS sería un buen libro, pero quizás demasiado técnico para un no experto en estadística).

Editar: Aquí hay un espectrograma ficticio representativo de los datos, con cortes a lo largo del tiempo y las dimensiones de longitud de onda.

El objetivo práctico aquí es evaluar la tasa de descomposición exponencial en el tiempo para cada longitud de onda (o bins, si es demasiado ruidoso).

Necesita un para especificar un modelo que separe la señal del ruido.

Existe el componente de ruido en el nivel de medición que se supone gaussiano. Los otros componentes, que dependen de las mediciones:

• "Este ruido, para una frecuencia fija, probablemente se puede modelar como aleatorio con distribución gaussiana". Necesita aclaración: ¿es el componente de ruido común a todos los puntos de tiempo, dada la frecuencia? ¿La desviación estándar es la misma para todas las frecuencias? Etc.

• "Sin embargo, en un momento fijo, los datos muestran un tipo diferente de ruido, con picos espurios grandes y oscilaciones rápidas" ¿Cómo se separa eso de la señal, porque se supone que le interesa la variación de la intensidad a través de la frecuencia? ¿Es la variación interesante de alguna manera diferente de la variación poco interesante, y si es así, cómo?

Las oscilaciones espurias o el ruido no gaussiano en general no es un gran problema, si tiene una idea realista de sus características. Se puede modelar transformando los datos (y luego usando un modelo gaussiano) o explícitamente usando una distribución de error no gaussiana. El ruido de modelado que se correlaciona con las mediciones es más difícil.

Dependiendo de cómo sea su modelo de ruido y datos, es posible que pueda modelar los datos con una herramienta de uso general como los GAM en el paquete mgcv, o puede necesitar una herramienta más flexible, lo que fácilmente conduce a una configuración bayesiana bastante personalizada . Hay herramientas para tales modelos, pero si no eres un estadístico, aprender a usarlos tomará un tiempo.

Supongo que una solución específica para el análisis espectral o el paquete mgcv son sus mejores apuestas.

Una serie temporal de espectros me sugiere un experimento cinético , y existe una cantidad bien establecida de literatura quimiométrica sobre esto.

¿Qué sabes sobre los espectros? ¿Qué tipo de espectros son? ¿Puede esperar razonablemente que solo tenga dos especies, educt y producto?

¿Puede suponer razonablemente la bilinealidad, es decir, los espectros medidos en un momento dado son una combinación lineal de las concentraciones de componentes los espectros de componentes puros : $\mathbf X$$\mathbf C$$\mathbf S$

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

Dices que quieres estimar una disminución exponencial (en las concentraciones). Esto, junto con la bilinealidad, me sugiere una resolución de curva multivariada (MCR). Esta es una técnica que le permite utilizar la información que tiene (p. Ej., Espectros de componentes puros de algunas sustancias o suposiciones sobre el comportamiento de concentración, como la disminución exponencial) durante el ajuste del modelo.

• Un buen punto de partida en la literatura es:
  Anna de Juan, Marcel Maeder, Manuel Martínez Romà Tauler: Combinando modelado duro y blando para resolver problemas cinéticos, quimiometría y sistemas inteligentes de laboratorio 54, 2000. 123–141.
• Un buen punto de partida en Internet es http://www.mcrals.info/
• El UseR! El libro de quimiometría lo analiza al final
• Implementación en R: paquete ALS

Hasta donde yo sé, es bastante común suavizar las concentraciones de acuerdo con algunos, por ejemplo, el modelo cinético, pero es mucho menos común suavizar los espectros. Sin embargo, el algoritmo permite hacerlo. Le pregunté a Anna en verano si imponen restricciones de suavidad, pero ella me dijo que no (y los buenos espectroscopistas odian suavizar en lugar de medir buenos espectros ;-)). A menudo, tampoco es necesario, porque la agregación de la información de todos los espectros ya producirá buenas estimaciones de los espectros de componentes puros.

Suavicé "espectros de componentes" (de hecho, componentes principales) dos veces últimamente ( Dochow et al .: Dispositivo Raman en chip y fibras de detección con rejilla de fibra Bragg para análisis de soluciones y partículas, LabChip, 2013 y Dochow el al. : Chip microfluídico de cuarzo para la identificación de células tumorales por espectroscopía Raman en combinación con trampas ópticas, AnalBioanalChem, aceptado), pero en estos casos mi conocimiento espectroscópico me dijo que se me permite hacer esto. Regularmente aplico una interpolación de reducción de muestreo y suavizado a mis espectros Raman ( hyperSpec::spc.loess).

¿Cómo saber qué es demasiado suavizado? Creo que la única respuesta posible es "conocimiento experto sobre el tipo de espectroscopia y experimento".

editar: volví a leer la pregunta y dices que quieres estimar la desintegración en cada longitud de onda. Sin embargo, ¿es cierto o desea estimar la descomposición de diferentes especies con espectros superpuestos?

Los datos consisten en espectros ópticos (intensidad de luz contra frecuencia) registrados en> tiempos variables. Los puntos se adquirieron en una cuadrícula regular en x (tiempo), y (frecuencia).

Para mí, esto suena como un caso para el análisis de datos funcionales (FDA), aunque no tengo idea de la física detrás de su problema, y ​​podría estar completamente equivocado. Si puede considerar que el proceso detrás de sus datos es inherentemente fluido y continuo, es posible que desee utilizar una expansión de función base bivariada para capturar sus mediciones en la forma , siendo una suma de funciones básicas (p. ej. b-splines) y coeficientes. Un conjunto limitado de funciones básicas reduce directamente la rugosidad y, por lo tanto, cancela una buena parte del ruido blanco.$intensity = f(time, frequency)$$f$

He leído cosas sobre la estimación de densidad de kernel 2D, la interpolación polinómica / spline 2D, etc.

...

Uso R, para lo cual veo muchos paquetes que parecen relacionados (MASS (kde2), campos (smooth.2d), etc.) pero no puedo encontrar muchos consejos sobre qué técnica aplicar aquí.

Usted mencionó la interpolación de spline, pero no mencionó el paquete fda que implementa muy bien y fácilmente accesible la expansión de la función base que mencioné anteriormente. El conjunto de mediciones simultáneas de tiempo, frecuencia e intensidad (ordenadas como una matriz tridimensional) podría capturarse como un objeto de datos funcionales bivariado, ver. por ejemplo, la función 'Data2fd'. Además, hay varios procedimientos de suavizado disponibles en el paquete, todos diseñados para cancelar el ruido blanco o la "rugosidad" en las mediciones de procesos inherentemente suaves.

El artículo de Wikipedia expresa el problema del ruido blanco en la FDA de la siguiente manera:

Los datos pueden ser tan precisos que el error puede ignorarse, puede estar sujeto a un error de medición sustancial o incluso tener una relación indirecta compleja con la curva que definen. ... los registros diarios de precipitación en una estación meteorológica son tan variables que requieren análisis cuidadosos y sofisticados para extraer algo así como una curva de precipitación media.

La FDA proporciona las herramientas para estos casos. ¿No se traduce esto en tu caso?

... pero no estoy familiarizado con la jerga o la teoría estadística subyacente ...

... pero no puedo encontrar muchos consejos sobre qué técnica aplicar aquí ...

Con respecto a la FDA: yo tampoco, pero el libro de Ramsay y Silverman sobre la FDA (2005) hace que los conceptos básicos sean muy accesibles y Ramsay Hooker y Graves (2009) traducen directamente las ideas del libro al código R. Ambos volúmenes deben estar disponibles como libros electrónicos en una biblioteca universitaria para estadísticas, biociencias, climatología o psicología. Google también mostrará algunos enlaces más que no puedo publicar aquí.

Lo siento, no puedo proporcionar una solución más directa para su problema. Sin embargo, la FDA me ayudó mucho una vez que descubrí para qué sirve.

Siendo un físico simple, no un experto en estadística, tomaría un enfoque simple. Las dos dimensiones son de diferentes naturalezas. Tendría sentido suavizar el tiempo con un algoritmo y suavizar la longitud de onda con otro.

Los algoritmos reales que usaría: para longitud de onda, Savitzky-Golay con un orden superior, 6 quizás 8.

A lo largo del tiempo, si ese ejemplo es típico, ese salto repentino y una disminución más o menos exponencial lo hacen complicado. He tenido datos experimentales e imágenes ruidosas, así como así. Si los métodos simples y directos no ayudan lo suficiente, pruebe con un suavizador gaussiano pero suprima su efecto cerca del salto, como lo detecta un detector de bordes. Suavice y amplíe la salida del detector de bordes, normalícelo para que pase de 0.0 a 1.0, y úselo para seleccionar entre la imagen original y la suavizada por Gauss, píxel por píxel.

@baptiste: Me alegro de que hayas agregado la trama como sugerí. Ayuda mucho:

Entonces, si entiendo correctamente, su objetivo práctico es evaluar la tasa de disminución exponencial para cada longitud de onda; entonces hagamos eso! Defina una función que desee minimizar para cada longitud de onda por separado y minimícela.

Veamos una sola longitud de onda dada, como en su gráfica inferior derecha.

Primero, por simplicidad, desechemos todos los valores antes de 0.2 segundos, porque contienen una discontinuidad masiva (nuestro enfoque se puede aumentar para lidiar con eso más adelante). Luego, defina el siguiente criterio de optimización, que tiene como objetivo encontrar la constante de desintegración :$\tau$

Puede resolver este problema de optimización analíticamente diferenciando wrt , igualando a cero y resolviendo para ; o puedes usar un solucionador.$\tau$$\tau$

Más adelante, si cree que la longitud de onda adyacente debería tener constantes de disminución similares, puede incorporar esto en un criterio de optimización más elaborado.

En todo caso, le sugiero que lea un libro de optimización que debe leer: la optimización convexa de Boyd .

¡Espero que esto ayude!