Re: [R] Fwd: Using odesolve to produce non-negative solutions 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.