naive implementation of linear sinesoidal regression

linear-sinesoidal-regression.R

# Kyle Andrews
# This code implements math derived by Jean Jacquelin
# License: GPL3+

rm(list = ls())
library(pracma)

x <- sort(runif(20, max = 10))
n <- length(x)
a0 <- 0
b0 <- 2
c0 <- -1
p0 <- 1.5
w0 <- pi/6
y <- a0 + p0*x + b0*sin(w0*x) + c0*cos(w0*x) + rnorm(n, sd = 1)

frm0 <- paste(round(a0, 3), paste0("(", round(b0, 3), ")*", "sin(", round(w0, 3), "*x)"), paste0("(", round(c0, 3), ")*", "cos(", round(w0, 3), "*x)"), sep = "+")

plot(x, y, main = paste("simulated data from", frm0))
grid()

S <- cumtrapz(x, y)
A <- cumtrapz(x, S)
E <- x^3
B <- x^2
C <- x
D <- rep(1, n)

imod <- lm(y ~ A + E + B + C + D - 1)


# initial estimates
rm(A, E, B, C, D)
co <- as.list(coef(imod))
w1 <- with(co, sqrt(-1*A))
a1 <- with(co, 2*B/w1^2)
p1 <- with(co, 6*E/w1^2)

pp1 <- x
aa1 <- rep(1, n)
be1 <- sin(w1*x)
nu1 <- cos(w1*x)

imod2 <- lm(y ~ aa1 + pp1 + be1 + nu1 - 1)
co2 <- as.list(coef(imod2))

b1 <- with(co2, be1)
c1 <- with(co2, nu1)

# thats all you need to get started

make_function <- function(a, p, b, c, w) function(x) a+p*x+b*sin(w*x)+c*cos(w*x)

x1 <- seq(min(x), max(x), length.out=100)
y1 <- make_function(a1, p1, b1, c1, w1)(x1)

lines(x1, y1, col = "red")

# change variables
r1 <- sqrt(b1^2 + c1^2)
ph1 <- atan(c1/b1)
if(b1 < 0) ph1 <- ph1 + pi

K <- round((w1*x+ph1)/pi)
r <- y - a1 - p1*x

th <- ifelse(r1^2 > r^2, (-1)^K*atan(r/sqrt(r1^2-r^2)) + pi*K,
             ifelse(r > 0, pi/2*(-1)^K + pi*K, -pi/2*(-1)^K + pi*K))

w <- x
ph <- rep(1, n)
mod3 <- lm(th ~ w + ph - 1)
co3 <- as.list(coef(mod3))

w2 <- with(co3, w)
ph2 <- with(co3, ph)

w3 <- w2
p3 <- p1
a3 <- a1
b3 <- r1*cos(ph2)
c3 <- r1*sin(ph2)

y3 <- make_function(a3, p3, b3, c3, w3)(x1)
lines(x1, y3, col = "blue")

legend(1, 12, legend=c("short", "full"),
       col=c("red", "blue"), lty=1, cex=0.8)