Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
HarithaJagadeesh opened this issue Aug 17, 2023 · 1 comment

Comments

@HarithaJagadeesh
Copy link

HarithaJagadeesh commented Aug 17, 2023

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]+delta*x[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 wellharmonic_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
Copy link
Owner

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants