robjhyndman / forecast

Forecasting Functions for Time Series and Linear Models
http://pkg.robjhyndman.com/forecast
1.11k stars 341 forks source link

I want to fit a harmonic regression model to the Lotka-Volterra model #944

Closed HarithaJagadeesh closed 1 year ago

HarithaJagadeesh commented 1 year ago

I am trying o fit a harmonic regression to the training data of prey population, predict it in the test data and find the RMSE value so that I can compare it with the other fits ( like smoothing spline) and determine which one is better. Inorder to do that, I tried using the lm() command but the prediction was coming in a very small range and I could not figure out why. So I tried making it a ts() object and fit it but then I was getting errors with the fourier command saying K must be not be greater than period/2 and Number of periods does not match number of orders. Can anyone help me with this?

This is the data : `library(ggplot2) require(mvtnorm) library(KGode)

noise= 0.1 ## set the variance of noise SEED = 19537 set.seed(SEED)

Define ode function, we use lotka-volterra model in this example.

we have two ode states x[1], x[2] and four ode parameters alpha, beta,

gamma and delta.

LV_fun = function(t,x,par_ode){ alpha=par_ode[1] beta=par_ode[2] gamma=par_ode[3] delta=par_ode[4] as.matrix( c( alphax[1]-betax[2]x[1] , -gammax[2]+deltax[1]x[2] ) ) }

create an ode class object

kkk0 = ode$new(2,fun=LV_fun)

set the initial values for each state at time zero.

xinit = as.matrix(c(0.5,1)) tinterv = c(0,120)

solve the ode numerically using alpha=1, beta=1, gamma=4, delta=1.

kkk0$solve_ode(c(1,1,4,1),xinit,tinterv)

Add noise to the numerical solution of the ode model and

treat it as the noisy observation.

y_true= t(kkk0$y_ode) n_o = max( dim( y_true) ) t_no = kkk0$t #time points for each observation y_no = y_true + rmvnorm(n_o,c(0,0),noise*diag(2))

splitting to training and test set

set.seed(123) data = data.frame(t_no,y_no) train_prop =0.8 train_index <- createDataPartition(1:nrow(data), p = train_prop, list = FALSE)

Split the data into training and test sets

train_data <- data[train_index, ] test_data <- data[-train_index, ]

Separate the prey and predator populations in the training set

train_data_prey <- train_data[, c(1, 2)]
train_data_predator <- train_data[, c(1, 3)]

Separate the prey and predator populations in the test set

test_data_prey <- test_data[, c(1, 2)]
test_data_predator <- test_data[, c(1, 3)] This is the harmonic model I tried but the prediction were not coming out well harmonic_model <- lm(X2 ~ sin((2pit_no)/6) + cos((2pit_no)/6), data = train_data_predator) summary(harmonic_model)

Predict using the Harmonic Regression Model

test_data_predator$predicted <- predict.lm(harmonic_model, newdata = test_data_predator, type = "response")`

Can someone please help me? Thanks, Haritha

robjhyndman commented 1 year ago

This is not the place to ask for help. Please post your question at http://crossvalidated.com or at https://community.rstudio.com/