Hi Spencer,
Thank you for your response. I also did not see anything on the lsoda help page
which is the reason that I wrote to the list.
>From your response, I am not sure if I asked my question clearly.
I am modeling a group of people (in a variety of health states) moving through
time (and getting infected with an infectious disease). This means that the
count of the number of people in each state should be positive at all times.
What appears to happen is that lsoda asks for a derivative at a given point in
time t and then adjusts the state of the population. However, perhaps due to
numerical instability, it occasionally lower the population count below 0 for
one of the health states (perhaps because it's step size is too big or
something).
I have tried both the logarithm trick and also changing the relative and
absolute tolerance inputs but I still get the problem for certain combinations
of parameters and initial conditions.
It occurs both under MS Windows XP Service Pack 2 and on a Linux cluster so I
am pretty sure it is not platform specific.
My real question to the group is if there is not a work around in lsoda are
there other ode solvers in R that will allow the constraint of solutions to the
ODEs remain non-negative?
Best regards,
Jeremy
>>> Spencer Graves <[EMAIL PROTECTED]> 6/8/2007 9:51 AM >>>
On the 'lsoda' help page, I did not see any option to force some
or all parameters to be nonnegative.
Have you considered replacing the parameters that must be
nonnegative with their logarithms? This effective moves the 0 lower
limit to (-Inf) and seems to have worked well for me in the past.
Often, it can even make the log likelihood or sum of squares surface
more elliptical, which means that the standard normal approximation for
the sampling distribution of parameter estimates will likely be more
accurate.
Hope this helps.
Spencer Graves
p.s. Your example seems not to be self contained. If I could have
easily copied it from your email and run it myself, I might have been
able to offer more useful suggestions.
Jeremy Goldhaber-Fiebert wrote:
> Hello,
>
> I am using odesolve to simulate a group of people moving through time and
> transmitting infections to one another.
>
> In Matlab, there is a NonNegative option which tells the Matlab solver to
> keep the vector elements of the ODE solution non-negative at all times. What
> is the right way to do this in R?
>
> Thanks,
> Jeremy
>
> P.S., Below is a simplified version of the code I use to try to do this, but
> I am not sure that it is theoretically right
>
> dynmodel <- function(t,y,p)
> {
> ## Initialize parameter values
>
> birth <- p$mybirth(t)
> death <- p$mydeath(t)
> recover <- p$myrecover
> beta <- p$mybeta
> vaxeff <- p$myvaxeff
> vaccinated <- p$myvax(t)
>
> vax <- vaxeff*vaccinated/100
>
> ## If the state currently has negative quantities (shouldn't have), then
> reset to reasonable values for computing meaningful derivatives
>
> for (i in 1:length(y)) {
> if (y[i]<0) {
> y[i] <- 0
> }
> }
>
> S <- y[1]
> I <- y[2]
> R <- y[3]
> N <- y[4]
>
> shat <- (birth*(1-vax)) - (death*S) - (beta*S*I/N)
> ihat <- (beta*S*I/N) - (death*I) - (recover*I)
> rhat <- (birth*(vax)) + (recover*I) - (death*R)
>
> ## Do we overshoot into negative space, if so shrink derivative to bring
> state to 0
> ## then rescale the components that take the derivative negative
>
> if (shat+S<0) {
> shat_old <- shat
> shat <- -1*S
> scaled_transmission <- (shat/shat_old)*(beta*S*I/N)
> ihat <- scaled_transmission - (death*I) - (recover*I)
>
> }
> if (ihat+I<0) {
> ihat_old <- ihat
> ihat <- -1*I
> scaled_recovery <- (ihat/ihat_old)*(recover*I)
> rhat <- scaled_recovery +(birth*(vax)) - (death*R)
>
> }
> if (rhat+R<0) {
> rhat <- -1*R
> }
>
> nhat <- shat + ihat + rhat
>
> if (nhat+N<0) {
> nhat <- -1*N
> }
>
> ## return derivatives
>
> list(c(shat,ihat,rhat,nhat),c(0))
>
> }
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
