Dado un conjunto de datos, ¿cómo se realiza una regresión sinusoidal en papel? ¿Cuáles son las ecuaciones, algoritmos?

Question

La mayoría de las regresiones son fáciles. Triviales una vez que sabes cómo hacerlo. La mayoría de ellas implican sustituciones que transforman los datos en una regresión lineal. Pero aún no he logrado averiguar cómo hacer una regresión sinusoidal. Estoy buscando el concepto más allá de los resultados. No necesito respuestas de Excel, TI o CAS. Me gustaría ver ecuaciones, métodos, etc. ¿Cómo lo harías en papel si realmente estuvieras dispuesto a hacerlo en papel?

---

Para aclarar, me refiero a la regresión sinusoidal más general $$ y = A\sin(Bx+C) + D$$

No algún caso especial asumiendo los valores de una de estas constantes.

Answer 1

El algoritmo de Gauss-Newton aborda directamente este tipo de problemas. Dados m puntos de datos $(x_i,y_i)$ para regresión con una función de n parámetros $\vec \beta =(\beta_1,...,\beta_n)$ $$min_{\vec \beta }\ S(\vec \beta)\ donde\ S(\vec \beta)=\sum_{i=1}^m r_i(\vec \beta)^2=(y_i-f(\vec \beta,x_i))^2$$ Me salto la derivación del algoritmo que puedes encontrar en cualquier libro de texto (primero se usa la aproximación de Taylor y luego se utiliza el método de Newton). $$\Delta \vec \beta=\big(J^T\ J\big)^{-1}\ J^T\ \vec r$$ $$\vec \beta=\vec \beta + \alpha\ \Delta \beta$$ donde $\alpha$ es el coeficiente de amortiguación y $$J=\begin{pmatrix}\bigg(\frac{\partial f}{\partial \beta_1}\bigg)_{x=x_1}&...&\bigg(\frac{\partial f}{\partial \beta_n}\bigg)_{x=x_1}\\\ ...&...&...\\\ \bigg(\frac{\partial f}{\partial \beta_1}\bigg)_{x=x_m}&...&\bigg(\frac{\partial f}{\partial \beta_n}\bigg)_{x=x_m}\end{pmatrix}\quad \vec r=\begin{pmatrix}y_1-f(\vec \beta,x_1)\\\ ... \\\ y_m-f(\vec \beta,x_m) \end{pmatrix}$$ Para tu caso específico $$\frac{\partial f}{\partial A}=sin(Bx_i+C)$$ $$\frac{\partial f}{\partial B}=Ax_icos(Bx_i+C)$$ $$\frac{\partial f}{\partial C}=Acos(Bx_i+C)$$ $$\frac{\partial f}{\partial D}=1$$

En Matlab generé 60 puntos de muestra aleatorios uniformemente distribuidos. Utilicé estos puntos para calcular una curva sin conocida como $$y=0.5sin(1.2x+0.3)+0.6$$ Añadí términos de error N(0,0.2) a cada punto. Mi suposición inicial fue $A=0.1\ B=0.5\ C=0.9\ D=0.1$ y configuré el coeficiente de amortiguación en 0.01. El algoritmo determina la siguiente ecuación de aproximación después de 407 iteraciones $$\hat y=0.497sin(1.178x+0.352)+0.580$$ A continuación puedes ver el gráfico (negro - curva original $y(x)$, rojo - curva con términos de error $y(x)+\epsilon (x)$, azul - curva de aproximación $\hat y(x)$)

Answer 2

Veamos cómo funciona la regresión.

Normalmente, se elige una función objetivo para tu curva de regresión (por ejemplo, una línea, una parábola, una curva logarítmica, etc.), y desarrollas alguna noción del error.

Tienes parámetros desconocidos en tu curva (por ejemplo, coeficientes polinomiales), y calculas los valores de estos parámetros de manera que tu función de error se minimice.

Para la regresión lineal simple, esto suele ser solo un problema de mínimos cuadrados. Para polinomios, puedes usar una matriz de Vandermonde y resolver un sistema lineal equivalente sin problema. Pero los polinomios son fáciles: son solo combinaciones lineales de potencias sucesivas de tus mediciones de variable independiente. Entonces, ¿qué hacemos cuando queremos ajustar, por ejemplo, $\sin(kt + b)$ a nuestros datos?

Bueno, depende. Digamos que tienes la creencia razonable de que tus datos se ajustan bien a una curva sinusoidal. Puedes calcular $\sin^{-1} y_i$ para todas tus mediciones (con una re-escala apropiada de $y_i$ al dominio de $\sin^{-1}$, y luego realizar una regresión lineal en $kt_i + b = \sin^{-1} y_i$.

Por supuesto, las funciones trigonométricas inversas suelen comportarse mal, por lo que esto podría no funcionar.

En su lugar, podrías pensar en cálculo, y recordar que $\sin$ es una función continua, e incluso, diferenciable. Una función diferenciable es bien aproximada por una función lineal en algún punto. Entonces, podrías calcular la derivada de tu función, y asumir que tus datos están localmente aproximados por muchas funciones lineales.

Otra forma es notar que $\sin$ tiene una expansión en serie de Taylor, y que suficientes términos deberían darte un polinomio que está bastante cerca de $\sin$ en algún dominio. Así que podrías realizar una expansión de la serie de Taylor de $n$ términos y hacer una regresión polinómica lineal regular sobre el resultado.

Si crees que tu función es una serie de senos, podrías escribir una expansión de serie de Fourier, y realizar un ajuste de mínimos cuadrados en los coeficientes de la serie de Fourier.

Si eso falla, podrías configurar una red neuronal para darte el resultado como una serie de senos. (En realidad, esto podría no funcionar sin pre-multiplicar cada término seno por, digamos, una Gaussiana, o una función de sombrero).

Puedes utilizar un algoritmo evolutivo para realizar una búsqueda estocástica en un espacio de parámetros, y usar la función de coste como tu criterio de supervivencia.

Finalmente, podrías emplear cualquier cantidad de algoritmos de mínimos cuadrados no lineales para estimar los parámetros de tu ajuste. Levenberg-Marquardt es un algoritmo comúnmente usado para estas cosas. Efectivamente, esto es lo mismo que hacer una linealización local de la función de regresión. Es una búsqueda basada en gradientes.

Esencialmente, todos los problemas de regresión son problemas de búsqueda: se buscan parámetros que den forma a la función objetivo de la manera más óptima. En consecuencia, cualquier algoritmo de búsqueda funcionará, pero no todos funcionarán bien. Desafortunadamente, no hay una respuesta bonita y cerrada elemental para calcular estos parámetros, como lo hay con la regresión lineal simple. Así que a menos que puedas transformar tus datos de alguna manera ingeniosa, es probable que no puedas hacerlo con papel y lápiz.

Answer 3

Um.

Si tienes un montón de puntos $(x_i,y_i)$ y quieres ajustar el mejor senoidal $A \sin (Bx +C),$ lo que haces es definir un error al cuadrado $$E = \sum_{i=1}^n \left( y_i - A \sin (B x_i + C) \right)^2.$$ Luego resuelves simultáneamente tres ecuaciones para el trío $(A,B,C),$ $$ \frac{\partial E}{\partial A} = 0, $$ $$ \frac{\partial E}{\partial B} = 0, $$ $$ \frac{\partial E}{\partial C} = 0. $$

Tendré que pensar si hay una solución en forma cerrada para las tres. Como se mencionó, en el caso lineal hay una solución en forma cerrada bien conocida. Si no la hay, bueno, dejaré esto aquí, ilustra un enfoque.

Answer 4

Por lo que sé, no hay un método exacto para la solución. Primero tienes que transformar $y$ de la siguiente manera $$y=A\sin(Bx+C)+D=A\sin(Bx)\cos(C)+B\cos(Bx)\sin(C)+D$$ $$y=\alpha \sin(Bx)+\beta \cos(Bx)+D\quad donde\ \alpha=A\cos(C)\ \beta=B\sin(C)$$ El primer método utiliza la reformulación de Prony para la frecuencia. En el método de Prony la función original es reemplazada por una función compleja. Introduce un polinomio complejo para encontrar las raíces en el dominio complejo y luego lo traduce de vuelta al dominio real. Después de elegir la frecuencia, el método determina otros parámetros mediante el "método de mínimos cuadrados". Tiene algunas desventajas. En primer lugar, debido a la reformulación de Prony, las muestras deben estar igualmente espaciadas. Además, durante la determinación de la frecuencia te encuentras con un sistema sobredeterminado y debes emplear algunas transformaciones que son bastante complejas (al menos, complejas para mí).

El segundo método reduce el problema al problema análogo para polinomios algebraicos y lo resuelve mediante el método de Forsythe. Usando este método puedes trabajar con muestras arbitrariamente distanciadas.

Por favor ten en cuenta que la información anterior proviene de mi investigación pasada; y el tema era sobre la facilidad de los métodos numéricos sobre los exactos; así que no soy un experto en estos métodos. Si deseas continuar tu búsqueda, te deseo mucha suerte. Por favor, ten en cuenta que la investigación es muy limitada en este tema.

Por otro lado, si deseas enfocarte en métodos numéricos puedo proporcionarte algoritmos (y métodos) para tu problema específico; y si me proporcionas puntos de datos un ejemplo sencillo (en Matlab).

Answer 5

Aquí está la versión de Matlab

damp = .01;

x = 0:.1:6;

y=0.5*sin(1.2*x+0.3)+0.6;

yn = y' + (rand(length(y),1) - .5)/10;

B=  [0.4; 1; 0.4; 0.5];

rows = length(yn);

for interation = 1:100000

    %construir J y r
    for row = 1:rows
         xk = x(row);

        J(row,1) = sin(B(2) * xk + B(3));

        J(row,2) = B(1) * xk * cos(B(2) * xk + B(3)); 

        J(row,3) = B(1) * cos(B(2) * xk + B(3)); 

        J(row,4) = 1.0;  

        r(row) =  -yn(row) +  B(1) * sin(B(2) * xk + B(3)) + B(4);
    end   

    Jt = transpose(J);    
    delta = (inv(Jt*J))*Jt*r';   
    B = B - damp*delta;