[R] Fwd: Using odesolve to produce non-negative solutions

[Martin Henry H. Stevens]

Hi Jeremy,
First, setting hmax to a small number could prevent a large step, if  
you think that is a problem. Second, however, I don't see how you can  
get a negative population size when using the log trick. I would  
think that that would prevent completely any negative values of N  
(i.e. e^-100000 > 0). Can you explain? or do you want to a void that  
trick? The only other solver I know of is rk4 and it is not recommended.
Hank
On Jun 11, 2007, at 11:46 AM, Jeremy Goldhaber-Fiebert wrote:

> 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 <spencer.graves at pdf.com> 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 at 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 at 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.



Dr. Hank Stevens, Assistant Professor
338 Pearson Hall
Botany Department
Miami University
Oxford, OH 45056

Office: (513) 529-4206
Lab: (513) 529-4262
FAX: (513) 529-4243
http://www.cas.muohio.edu/~stevenmh/
http://www.muohio.edu/ecology/
http://www.muohio.edu/botany/

"E Pluribus Unum"