Skip to content

[R-sig-dyn-mod] deSolve ODE, interpretation/defintion of result parameter

2 messages · Armin Boehrer, Thomas Petzoldt

#
Hello,
?
I wanted to start using ODE in deResolve.
I thought an own simple example would be best to learn:
?
The ODE solving for an exponential? y = a * exp(k*x)
and its derivative dy/dx = a*k * exp(k*x)
should be easy:?
??? Just define??? dy <- ak_fit * y
? ? Run code, see blow,
? ? and it works!
?
What I wonder/My question:
I expected to get the result ak_fit as the product of my input parameters a and k,
but I see that I get just the decay-constant k.
If I had not known the result, I would have misinterprested the result.
?
Can you help? What is wrong with my naive interpretation?
Thanks and best wishes,
Armin B?hrer
?
?
?

#############################
### Solve_aexpkx_ODE_optim_Rcode.txt
### a simple exponential
#############################
rm(list = ls())
set.seed(12345)
my_a <- 5.4321
my_k <- -0.12345
my_x <- seq(0,10,1)
my_y <- my_a*exp(my_k*my_x) * rnorm(11,1,0.0000001)
my_y_true <- my_a*exp(my_k*my_x)
# print factor for derivative
my_a*my_k
# [1] -0.6705927
plot(my_x,my_y)
points(my_x,my_y_true,pch=2,col='green')
# Load required package
library(deSolve)
?
# Define the differential equation
dydx <- function(x, y, ak_fit) {
? dy <- ak_fit * y
? list(dy)
}
?
# Parameters
parameters <- c(ak_fit = -0.1)
# Define the sum of squares function
cost_function <- function(parameters) {
? out <- ode(y = state, times = my_x, func = dydx, parms = parameters)
? predicted_y <- out[, "y_fit"]
? cost_y <- sum((my_y - predicted_y)^2)
? total_cost <- cost_y
? assign("total_cost_SSE", total_cost, envir = .GlobalEnv)
? return(total_cost)
}
# Initial conditions , here usually a guess is used, lets take my_y(0)
start_value <- 5.4
state <- c(y_fit = start_value)
# Use optim to find the best value for parameters
optim_result <- optim(par = parameters, fn = cost_function)
# Warning message:
# In optim(par = parameters, fn = cost_function) :
#?? one-dimensional optimization by Nelder-Mead is unreliable:
# use "Brent" or optimize() directly
print(total_cost_SSE)
# [1] 0.002106345
final_parameters <- optim_result$par
print(final_parameters)
#???? ak_fit
# -0.1223438
outout <- ode(y = state, times = my_x, func = dydx, parms = final_parameters)
points(outout[,1],outout[,2],col='red'??? ,pch=3)
?
### ### optimize with optimize
# Define the sum of squares function for start value
cost_function_start <- function(start_value) {
? state <- c(y_fit = start_value)
? assign("state", state, envir = .GlobalEnv)
? optim_result <- optim(par = parameters, fn = cost_function)
? total_cost_start <- total_cost_SSE
? return(total_cost_start)
}
# Use optimize to find the best value for parameters
result_optimize <- optimize(cost_function_start, c(5.3,5.6), tol = 0.000001)
# Es gab 15 Warnungen (Anzeige mit warnings())
final_y0 <- result_optimize$minimum
final_y0
# [1] 5.432116
state <- c(y_fit = final_y0)
optim_result <- optim(par = parameters, fn = cost_function, method=c("Brent"), lower=c(-0.2), upper=c(-0.1))
print(total_cost_SSE)
# [1] 7.59912e-10
final_parameters <- optim_result$par
print(final_parameters)
# [1] -0.1234507
outout <- ode(y = state, times = my_x, func = dydx, parms = final_parameters)
points(outout[,1],outout[,2],col='blue'??? ,pch=16)
?
8 days later
#
Hello,

the original poster sent me a message, that we wants to cancel his post:

Am 27.08.2024 um 16:40 schrieb Armin Boehrer:
?[...]
Thomas Petzoldt
(moderator)



Am 19.08.2024 um 07:24 schrieb Armin Boehrer:
?[...]