Skip to content

Instantly share code, notes, and snippets.

@wdkrnls
Created January 4, 2018 06:19
Show Gist options
  • Save wdkrnls/42a690e70703be205db4311491b8c4b4 to your computer and use it in GitHub Desktop.
Save wdkrnls/42a690e70703be205db4311491b8c4b4 to your computer and use it in GitHub Desktop.
naive implementation of linear sinesoidal regression
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment