Re: [R] Fwd: Using odesolve to produce non-negative solutions
By the way, if someone could forward the original question to me (I'm subscribed to but not currently receiving R-help, as I found I was spending too much time reading it!) I might think of something more useful. (alternatively, when was it posted; I can find it on gmane, too). Woody R. Woodrow Setzer, Ph. D. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B205-01/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-1194 Martin Henry H. Stevens [EMAIL PROTECTED]To .eduSpencer Graves [EMAIL PROTECTED] 06/11/2007 01:02cc PM Jeremy Goldhaber-Fiebert [EMAIL PROTECTED], R-Help r-help@stat.math.ethz.ch, Woodrow Setzer/RTP/USEPA/[EMAIL PROTECTED] Subject Re: [R] Fwd: Using odesolve to produce non-negative solutions Hi Spencer, I have copied Woody Setzer. I have no idea whether lsoda can estimate parameters that could take imaginary values. Hank On Jun 11, 2007, at 12:52 PM, Spencer Graves wrote: in line Martin Henry H. Stevens wrote: 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. SG: Can lsoda estimate complex or imaginary parameters? Hmm. I have no idea. I would think that that would prevent completely any negative values of N (i.e. e^-10 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 snip __ 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. 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 __ 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.
Re: [R] Fwd: Using odesolve to produce non-negative solutions
Hi, all. lsoda can certainly not handle complex parameters. You can try (as Hank suggested) limiting hmax. You can also crank up relative and absolute precision by specifying smaller values of rtol and atol. I've seen similar problems in which the state variable becomes negative, with very small absolute value, when theoretically, the system has a non-negative solution. This is certainly due to imprecision in the numerical solution. Have you tried including an analytic jacobian? That could improve the numeric properties of the solution. Woody R. Woodrow Setzer, Ph. D. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B205-01/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-1194 Martin Henry H. Stevens [EMAIL PROTECTED]To .eduSpencer Graves [EMAIL PROTECTED] 06/11/2007 01:02cc PM Jeremy Goldhaber-Fiebert [EMAIL PROTECTED], R-Help r-help@stat.math.ethz.ch, Woodrow Setzer/RTP/USEPA/[EMAIL PROTECTED] Subject Re: [R] Fwd: Using odesolve to produce non-negative solutions Hi Spencer, I have copied Woody Setzer. I have no idea whether lsoda can estimate parameters that could take imaginary values. Hank On Jun 11, 2007, at 12:52 PM, Spencer Graves wrote: in line Martin Henry H. Stevens wrote: 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. SG: Can lsoda estimate complex or imaginary parameters? Hmm. I have no idea. I would think that that would prevent completely any negative values of N (i.e. e^-10 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 snip __ 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. 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 __ 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.
Re: [R] problems with optim, for-loops and machine precision
Without more detail - a reproducible example - it is hard to give you
concrete advice.
I wonder if the functions NLL23 and NLL21 depend on numerical solutions
of a system of ODEs, since you invoke the odesolve package? If so, try
switching to the Nelder-Mead optimizer, enforcing the parameter
constraints using transformation. Probably you are using the finite
difference derivatives calculated internally to optim for the gradient
information used in the L-BFGS-B optimizer. These can be unstable when
based on numerical solutions of odes, causing the optimizer to fail, or
sometimes to converge to a non-optimal point.
Some other points:
- you cannot change machine precision by changing values in .Machine.
To change the number of digits printed, use options(digits=8).
- use 'library()' instead of 'require()', unless you need to use the
return value from 'require()'
R. Woodrow Setzer, Ph. D.
National Center for Computational Toxicology
US Environmental Protection Agency
Mail Drop B205-01/US EPA/RTP, NC 27711
Ph: (919) 541-0128Fax: (919) 541-1194
Simon Ruegg
[EMAIL PROTECTED]
unizh.ch To
Sent by: r-help@stat.math.ethz.ch
[EMAIL PROTECTED]cc
tat.math.ethz.ch
Subject
[R] problems with optim,
01/10/2007 07:18 for-loops and machine precision
AM
Dear R experts,
I have been encountering problems with the optim routine using for
loops. I am determining the optimal parameters of several nested models
by
minimizing the negative Log-Likelihood (NLL) of a dataset.
The aim is to find the model which best describes the data. To this end,
I
am simulating artificial data sets based on the model with the least
number
of parameters (6) and the parameters determined with the field data. For
each artificial set I estimate the parameters of the model with 6
parameters
and the next more complex model with 7 parameters (two of these
parameters
are equal in the 6-parameter model) by minimizing the corresponding NLL
with
optim. In theory the 7-parameter model should fit the data either
equally
or better than the 6-parameter model. Therefore the difference of the
minimal NLLs should be 0 or larger.
For 500 data sets I use the following code:
require(odesolve)
res=matrix(nrow=500,ncol=18)
colnames(res)=c(a_23,beta_23,mu_23,d_23,I_23,M_23,NLL_23,
a_21,beta_21,mu_21,e_21,d_21,I_21,M_21,NLL_21,
NLL23_min_21,conv23,conv21)
for (s in 1:500)
{
assign(data,read.table(paste(populations/TE23simset_,s,.txt,sep=)),e
nv=MyEnv) #reading a data set
M23=optim(rep(0.1,6),NLL23,method=L-BFGS-B,lower=0,
upper=c(Inf,Inf,Inf,Inf,1,1),control=c(maxit=150))
if (M23$convergence==0)
{
M21=optim(rep(0.1,7),NLL21,method=L-BFGS-B,lower=0,
upper=c(Inf,Inf,Inf,Inf,Inf,1,1),control=c(maxit=150))
}
res[s,1]=M23$par[1]
res[s,2]=M23$par[2]
res[s,3]=M23$par[3]
res[s,4]=M23$par[4]
res[s,5]=M23$par[5]
res[s,6]=M23$par[6]
res[s,7]=M23$value
res[s,8]=M21$par[1]
res[s,9]=M21$par[2]
res[s,10]=M21$par[3]
res[s,11]=M21$par[4]
res[s,12]=M21$par[5]
res[s,13]=M21$par[6]
res[s,14]=M21$par[7]
res[s,15]=M21$value
res[s,16]=M23$value-M21$value
res[s,17]=M23$convergence
res[s,18]=M21$convergence
write.table(res,compare23_21TE061221.txt)
rm(M23,M21)
}
}
For some strange reason the results do not correspond to what I expect:
about 10% of the solutions have a difference of NLL smaller than 0. I
have
verified the optimisation of these results manually and found that a
minimal
NLL was ignored and a higher NLL was returned at convergence. To check
what was happening I inserted a printing line in the NLL function to
print
all parameters and the NLL as the procedure goes on. To my surprise
optim
then stopped at the minimal NLL which had been ignored before. I have
then
reduced the machine precision to .Machine$double.digits=8 thinking, that
the
printing was slowing down the procedure and by reducing the
Re: [R] How to solve differential equations with a delay (time lag)?
lsoda does not solve delay differential equations, which is what you
need to be able to do. A quick search on netlib (netlib.org) turned up
one match for delay differential equations: ode/ddverk.f, which might
help (though not in R).
R. Woodrow Setzer, Ph. D.
National Center for Computational Toxicology
US Environmental Protection Agency
Mail Drop B205-01/US EPA/RTP, NC 27711
Ph: (919) 541-0128Fax: (919) 541-1194
Ravi Varadhan
[EMAIL PROTECTED]
eduTo
Sent by: r-help@stat.math.ethz.ch
[EMAIL PROTECTED]cc
tat.math.ethz.ch
Subject
[R] How to solve differential
11/29/2006 04:44 equations with a delay (time
PM lag)?
Hi,
I would like to solve a system of coupled ordinary differential
equations,
where there is a delay (time lag) term. I would like to use the lsoda
function odesolve package. However, I am not sure how to specify the
delay term using the syntax allowed by odesolve.
Here is an example of the kind of problem that I am trying to solve:
library(odesolve)
yprime - function(t, y, p) { # this function
yd1 - p[k1] *(t = p[T]) - p[k2] * y[2]
yd2 - p[k3] * y[1](t - p[delay]) - p[k4] * y[2] # this is not
syntactically valid, but it is what I would like to do
list(c(yd1,yd2))
}
times - seq(0,30,by=0.1)
y0 - c(0,0)
parms - c(k1=0.7, k2=0.5, k3=0.2, k4=0.8, T=10, delay=5)
Is there a way to incorporate delay in odesolve?
Any hints would be much appreciated.
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
[[alternative HTML version deleted]]
__
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.
Re: [R] nls convergence problem
Joerg van den Hoff [EMAIL PROTECTED] wrote on 08/16/2006 08:22:03 AM: Earl F. Glynn wrote: [deleted] efg [deleted] (I think this is recognized by d. bates, but simply way down his 'to do' list :-(). joerg No doubt Doug Bates would gladly accept patches ... . __ 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. Woodrow Setzer, Ph. D. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B205-01/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-1194 __ 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.
Re: [R] Error when loading odesolve
Prof. Ripley: Thanks for diagnosing Wuming Gong's problem. I'm not sure I would have recognized the solution so quickly. I have uploaded to CRAN a new version of odesolve with Makevars fixed. R. Woodrow Setzer, Ph. D. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B205-01/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-1194 Prof Brian Ripley [EMAIL PROTECTED]To .ac.uk Wuming Gong Sent by: [EMAIL PROTECTED] [EMAIL PROTECTED]cc tat.math.ethz.ch r-help@stat.math.ethz.ch Subject Re: [R] Error when loading 08/04/2006 01:33 odesolve PM Looking at odesolve, src/Makevars is PKG_LIBS=$(BLAS_LIBS) Now, the documentation says that if you have $(BLAS_LIBS) you must also have $(FLIBS), so please change this to PKG_LIBS=$(BLAS_LIBS) $(FLIBS) and take this up with the package maintainer (which is what the posting guide asked you to do in the first place). On Sat, 5 Aug 2006, Wuming Gong wrote: Dear list, I installed odesolve package (0.5-15) in R 2.3.1 in a Solaris server (Generic_118558-11 sun4u sparc SUNW,Sun-Blade-1000). The installing progress completed without errors, though several warnings like Warning: Option -fPIC passed to ld, if ld is invoked, ignored otherwise were outputed. However, when loading the odesolve package by library(odesolve), following error messages pop out: library(odesolve) Error in dyn.load(x, as.logical(local), as.logical(now)) : unable to load shared library '/project/scratch/ligroup/R1/lib/R/library/odesolve/libs/odesolve.so': ld.so.1: R: fatal: relocation error: file /project/scratch/ligroup/R1/lib/R/library/odesolve/libs/odesolve.so: symbol __f90_ssfw: referenced symbol not found Error: package/namespace load failed for 'odesolve' Could any one tell me how to fix this problem? Thanks very much. Wuming __ 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. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ 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.
Re: [R] odesolve loading problem
Hi Ariel, The problem is that I specified the wrong dependency in the DEPENDS field of the DESCRIPTION file of the package: I specified R version of at least 2.2.1, but that should have been 2.3.1. You have two choices -- upgrade your R to 2.3.1, or install odesolve 0.5.13. I will send an updated version of the package to CRAN, with a note to the CRAN maintainers about the problem, but that won't help you if you need to use R version 2.2.1. There is an archive of older R packages on CRAN, linked to at the bottom of the Contributed Packages page on CRAN. Please accept my apology for the inconvenience -- I rushed through a change requested by a user, and did not take time to fully appreciate the consequences. Woody R. Woodrow Setzer, Ph. D. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B205-01/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-1194 Ariel Chernomoretz [EMAIL PROTECTED]To r-help@stat.math.ethz.ch Sent by:cc [EMAIL PROTECTED] tat.math.ethz.ch Subject [R] odesolve loading problem 07/28/2006 01:30 AM Hi, I get the following error message when loading the package odesolve ( R 2.2.1 - odesolve 0.5.14 - AMD64 - Linux Fedora Core 4 ) : library(odesolve) Error in library.dynam(lib,package,package.lib) : shared library 'TRUE' not found Error: package/namespace load failed for 'odesolve' Any help would be greatly appreciated Ariel./ __ 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.
Re: [R] odesolve with banded Jacobian [was no subject]
Dear Karline Soetaert, I've just returned from a week of travel, so have not had a great deal of time to look at your request. From a brief rereading of the original lsoda documentation, it looks as if all I need to do is set a flag to a different value (jt to 4), and leave it up to the user to construct the function that calculates the jacobian properly. If you'd contact me directly, ideally with a test model, I will see if the modification is really that simple; if so, I'll make the change and release an updated odesolve to CRAN. Woody PS: Thanks, Martin R. Woodrow Setzer, Jr. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B205-01/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-1194 Martin Maechler [EMAIL PROTECTED] ath.ethz.chTo Soetaert, Karline 11/14/2005 09:46 [EMAIL PROTECTED] AM cc R-help@stat.math.ethz.ch, Woodrow Setzer/RTP/USEPA/[EMAIL PROTECTED] Please respondSubject toRe: [R] odesolve with banded Martin Maechler Jacobian [was no subject] [EMAIL PROTECTED] ath.ethz.ch KSoet == Soetaert, Karline [EMAIL PROTECTED] on Mon, 14 Nov 2005 13:20:24 +0100 writes: KSoet Hi, I am trying to solve a model that consists of KSoet rather stiff ODEs in R. KSoet I use the package ODEsolve (lsoda) to solve these KSoet ODEs. KSoet To speed up the integration, the jacobian is also KSoet specified. KSoet Basically, the model is a one-dimensional KSoet advection-diffusion problem, and thus the jacobian is KSoet a tridiagonal matrix. KSoet The size of this jacobian is 100*100. KSoet In the original package LSODA it is possible to KSoet specify that the jacobian is banded, which makes its KSoet inversion very efficient. KSoet However, this feature seems to have been removed in KSoet the R version. KSoet Is there a way to overcome this limitation? Yes. But probably not a very easy one; maybe even a very cumbersome one... ;-) Note however that questions like these should typically be addressed at the package author - which you can always quickly find out via packageDescription(odesolve) Package: odesolve Version: 0.5-12 Date: 2004/10/25 Title: Solvers for Ordinary Differential Equations Author: R. Woodrow Setzer [EMAIL PROTECTED] Maintainer: R. Woodrow Setzer [EMAIL PROTECTED] Depends: R (= 1.4.0) Description: This package provides an interface for the ODE solver lsoda. ODEs are expressed as R functions or as compiled code. ... I've CC'ed this e-mail to Woodrow to help you for once .. KSoet [[alternative HTML version deleted]] KSoet __ KSoet . KSoet PLEASE do read the posting guide! KSoet http://www.R-project.org/posting-guide.html if you do read that guide, it will tell you - why you should always use a 'Subject' for your e-mails - why HTML-ified e-mails are not much liked and what you can do about it. Regards, Martin Maechler, ETH Zurich __ 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
Re: [R] changing memory limits to speed up lsoda
Hank, I don't understand why you think memory is the problem, here. I'd try writing my model in C or Fortran (there is an example in the odesolve package). That speeds things up a lot, and is what I do with slow systems. R. Woodrow Setzer, Jr. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B205-01/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-1194 Martin Henry H. Stevens [EMAIL PROTECTED]To .eduWoodrow Setzer/RTP/USEPA/[EMAIL PROTECTED] cc 10/26/2005 10:45 AM Subject changing memory limits to speed up lsoda Hi Woody, I sent this to r-help as well. Thanks in advance for any input. I am running R 2.2.0 on Mac OS 10.4.2, dual G5 processors with 8 Gig RAM. I am running a simulation with lsoda that requires ~378 s to complete one set of time intervals. I need to optimize the parameters, and so need to considerably speed up the simulation. I have tried to figure out how to change the appropriate memory allocation and have search R help and Introductory information and the archives, but connot figure out how to allocate more to the right place. I would greatly appreciate any pointers or tips or leads. Thank you, Hank Stevens Dr. Martin Henry H. 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 __ 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
Re: [R] odesolve/lsoda differences on Windows and Mac
I've been talking offline with Hank Stephens about this; note that in the example he quotes, he set hmin = 0.1, and the quoted error message says that the stepsize had reached hmin with no convergence. I believe that he intended to set hmax (because of the pulsed input). Then, Peter Dalgaard's explanation would make sense -- the PPC platform just needs a smaller stepsize than does the PC to achieve convergence. R. Woodrow Setzer, Jr. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B305-03/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-4284 Martin Henry H. Stevens [EMAIL PROTECTED]To .edu'R-Help' r-help@stat.math.ethz.ch 07/27/2005 12:36cc PM Thomas Petzoldt [EMAIL PROTECTED], Woodrow Setzer/RTP/USEPA/[EMAIL PROTECTED] Subject odesolve/lsoda differences on Windows and Mac Hi - I am getting different results when I run the numerical integrator function lsoda (odesolve package) on a Mac and a PC. I am trying to simulating a system of 10 ODE's with two exogenous pulsed inputs to the system, and have had reasonably good success with many model parameter sets. Under some parameter sets, however, the simulations fail on the Mac (see error message below). The same parameter sets, however, appear to run fine for our computational technician on his PC, generating apparently very reasonable data. Our tech is successfully running Dell Latitude D810, Windows XP Pro (Service Pack 2), 1Gb RAM. RGUI 2.1.1 I am running: R Version 2.1.1 (2005-06-20) on a Mac OS 10.3.9 Machine Model:Power Mac G5 CPU Type: PowerPC 970 (2.2) Number Of CPUs: 2 CPU Speed:2 GHz L2 Cache (per CPU): 512 KB Memory: 1.5 GB Bus Speed:1 GHz Boot ROM Version: 5.0.7f0 Serial Number:XB3472Q1NVS My Error Message system.time( + outAc2 - as.data.frame(lsoda(xstart,times, pondamph, parms, tcrit=170*730, hmin=.1)) + ) [1] 0.02 0.01 0.04 0.00 0.00 Warning messages: 1: lsoda-- at t (=r1) and step size h (=r2), the 2: corrector convergence failed repeatedly 3: or with abs(h) = hmin 4: Returning early from lsoda. Results are accurate, as far as they go Thanks for any input. Hank Stevens Dr. Martin Henry H. 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/botany/bot/henry.html http://www.muohio.edu/ecology/ http://www.muohio.edu/botany/ E Pluribus Unum __ 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
Re: [R] odesolve/lsoda differences on Windows and Mac
I will have added most of the solvers from the ode package odepack by Alan Hindmarsh, LLNL, to odesolve this year. These are all solvers in the lsode family. I would also like to add solver(s) for DAEs, like daskr (Brown, Hindmarsh, Petzold), but that may take a bit longer. If other folks are planning to contribute solvers, I would be happy to discuss including them in odesolve, so there would be a single main package for solving ODEs. R. Woodrow Setzer, Jr. National Center for Computational Toxicology US Environmental Protection Agency Mail Drop B305-03/US EPA/RTP, NC 27711 Ph: (919) 541-0128Fax: (919) 541-4284 Thomas Petzoldt [EMAIL PROTECTED] z.tu-dresden.deTo Peter Dalgaard 07/28/2005 05:05 [EMAIL PROTECTED] AM cc Martin Henry H. Stevens [EMAIL PROTECTED], 'R-Help' r-help@stat.math.ethz.ch, Woodrow Setzer/RTP/USEPA/[EMAIL PROTECTED], Thomas Petzoldt [EMAIL PROTECTED] Subject Re: [R] odesolve/lsoda differences on Windows and Mac On 27 Jul 2005, Peter Dalgaard wrote: One thought: Integrating across input pulses is a known source of turbulence in lsoda. You might have better luck integrating over intervals in which the input function is continuous. Tweaking the lsoda tolerances is another thing to try. Yes, that's also our experience. Where I am usually succesful when playing with the tolerances or the interpolation rule of external pulses, some of our students use the fixed step rk4 algorithm and some others wrote their own integrators in R. I have heared that several people had plans to provide alternative ODE integrators for R but I currently do not know about the state of these projects. It wold be nice if they might post this to the list in order to avoid double work. I haven't seen lsoda fail like that, but it's not too surprising that marginal cases show platform dependency (i.e. the integrator just fails on Mac and barely succeeds on PC). Aha, I see. It should be regarded carefully when publishing examples that result in marginal cases as the common user would expect that R is platform independent. Thomas P. __ 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
[R] Using varPower in gnls, an answer of sorts.
Back on January 16, a message on R-help from Ravi Varadhan described a problem with gnls using weights=varPower(). The problem was that the fit failed with error Error in eval(expr, envir, enclos) : Object . not found I can reliably get this error in version 2.0.1-patched 2004-12-09 on Windows XP and 2.0.1-Patched 2005-01-26 on Linux. The key feature of that example is that the data are being passed in the environment. Consider a modification of the example in the man page for gnls: First, something that should work: gnls(weight ~ Asym/(1 + exp((xmid - Time)/scal)),data=Soybean, + start=c(Asym=16,xmid=50,scal=7),weights=varPower()) Generalized nonlinear least squares fit Model: weight ~ Asym/(1 + exp((xmid - Time)/scal)) Data: Soybean Log-likelihood: -486.8973 Coefficients: Asym xmid scal 17.35681 51.87230 7.62052 Variance function: Structure: Power of variance covariate Formula: ~fitted(.) Parameter estimates: power 0.8815438 Degrees of freedom: 412 total; 409 residual Residual standard error: 0.3662752 ## Now, use with() to pass Soybean in the environment of gnls: with(Soybean,gnls(weight ~ Asym/(1 + exp((xmid - Time)/scal)), + start=c(Asym=16,xmid=50,scal=7),weights=varPower())) Error in eval(expr, envir, enclos) : Object . not found ## drop the weights argument from gnls, and the error message goes away. The problem is in a call to model.frame. When varPower() (and presumably other weight functions using '.' to represent a fitted model) is used, gnls() constructs a formula argument for model.matrix that looks like: ~.+weight+Time This works when the data are passed in a data frame, but not when they are in the environment. Look at the example in the man page for model.frame data.class(model.frame(~.+dist+speed, data=cars)) [1] data.frame data.class(with(cars,model.frame(~.+dist+speed))) Error in eval(expr, envir, enclos) : Object . not found data.class(with(cars,model.frame(~dist+speed))) [1] data.frame env - new.env(TRUE,NULL) assign(dist,cars$dist,envir=env) assign(speed,cars$speed,envir=env) data.class(model.frame(~dist+speed,data=env)) [1] data.frame data.class(model.frame(~.+dist+speed,data=env)) Error in eval(expr, envir, enclos) : Object . not found I'm sending this to r-help rather to the authors of nlme because I am not sure where the bug is: is model.frame misbehaving, or should gnls not include '.' in its call to model.frame? R. Woodrow Setzer, Jr.Phone: (919) 541-0128 Experimental Toxicology Division Fax: (919) 541-4284 Pharmacokinetics Branch NHEERL B143-01; US EPA; RTP, NC 27711 __ 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
Re: [R] odesolve: lsoda vs rk4
lsoda doesn't pass along the names attribute of the state vector, y. If
you want to use the names of the state vector in your code, you need to
reassign it inside your ode function. I should either fix this in the
code for lsoda, or at least document it! For example, I can run the
following modification of your code:
func2 - function(t, y, p)
{
names(y) - c(A,B,C,D) ### Here put the names attribute back
Ad - p[p2]*(p[p1]*y[A]*y[D])/(p[p2]+p[p3]) +
p[p6]*(p[p4]*y[B]*p[p10])/(p[p5]+p[p6]) -
p[p1]*y[A]*y[C]
Bd - p[p3]*(p[p1]*y[A]*y[D])/(p[p2]+p[p3]) +
p[p5]*(p[p4]*y[B]*p[p10])/(p[p5]+p[p6]) -
p[p4]*y[B]*p[p10]
Cd - (p[p1]+p[p7])*y[A]*y[D] -
p[p1]*y[A]*y[C]-p[p9]*y[C]
Dd -p[p9]*y[C] - p[p7]*y[A]*y[D]
list(c(Ad, Bd, Cd, Dd))
}
out2 - lsoda(y, times, func2, parms)
out2[c(1:10,1000),]
timeABC D
[1,] 0.00 2.50e-06 2.50e-06 1.70e-06 5.70e-07
[2,] 0.05 2.432561e-06 2.500379e-06 1.671689e-06 5.312515e-07
[3,] 0.10 2.368078e-06 2.499286e-06 1.639363e-06 4.980014e-07
[4,] 0.15 2.306621e-06 2.496841e-06 1.603709e-06 4.697523e-07
[5,] 0.20 2.248381e-06 2.493370e-06 1.566611e-06 4.451401e-07
[6,] 0.25 2.193360e-06 2.488923e-06 1.528312e-06 4.239702e-07
[7,] 0.30 2.141477e-06 2.483726e-06 1.489881e-06 4.053217e-07
[8,] 0.35 2.092710e-06 2.477833e-06 1.451555e-06 3.889875e-07
[9,] 0.40 2.046809e-06 2.471378e-06 1.413729e-06 3.744582e-07
[10,] 0.45 2.003632e-06 2.464438e-06 1.376627e-06 3.614428e-07
[11,] 49.95 2.675890e-06 5.563057e-08 1.595362e-09 -7.457989e-11
In the event the results of lsoda and rk4 differ, I would trust the
lsoda result over that of rk4 (see the warnings in the help file for
rk4).
On June 10, 2004, Chris Knight wrote:
[snip ...]
I'm trying to use odesolve for integrating various series of coupled 1st
order differential equations (derived from a system of enzymatic
catalysis and copied below, apologies for the excessively long set of
parameters).
The thing that confuses me is that, whilst I can run the function rk4:
out - rk4(y=y,times=times,func=func, parms=parms)
and the results look not unreasonable:
out-as.data.frame(out)
par(mfrow=c(4,1))
for (i in 2:(dim(out)[2]))plot(out[,1],out[,i], pch=., xlab=time,
ylab=names(out)[i])
If I try doing the same thing with lsoda:
out - lsoda(y=y,times=times,func=func, parms=parms, rtol=1e-1, atol=
1e-1)
I run into problems with a series of 'Excessive precision requested'
warnings with no output beyond the first time point.
Fiddling with rtol and atol doesn't seem to do very much.
What is likely to be causing this (I'm guessing the wide range of the
absolute values of the parameters can't be helping), is there anything I
can sensibly do about it and, failing that, can I reasonably take the
rk4 results as being meaningful?
Any help much appreciated,
Thanks in advance,
Chris
[code deleted...]
~~
Dr. Christopher Knight Tel:+44 1865 275111
Dept. Plant Sciences +44 1865 275790
South Parks Road
Oxford OX1 3RB Fax:+44 1865 275074
` · . , ,(((º
~~
R. Woodrow Setzer, Jr.Phone: (919) 541-0128
Experimental Toxicology Division Fax: (919) 541-4284
Pharmacokinetics Branch
NHEERL B143-01; US EPA; RTP, NC 27711
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] lsoda with arbitrary zero thresholds (with psuedo-solution)
Dear Hank,
Last question first: really, only you can say for sure if 4e-281 and
5e-11 are small enough; it depends on the units you measure your state
variables in. However, this strategy cannot get the state variables to
exactly 0. Obviously, you could get closer to 0.0 faster by setting the
derivatives even larger in absolute value. You may run into problems
with the solver when the derivatives are discontinuous functions of the
state variables.
There is a simple and elegant solution in theory, but not (yet)
available in odesolve. soda has a variant called lsodar that returns
whenever a function of the state variables satisfies a given set of
conditions (in your case, you could tell lsodar to return whenever any
state variable drops below 0.4). Once the call to lsodar returns, you'd
then reset all the state variables that were 0.4 to 0, and restart the
integrator at that point. I've been meaning to add lsodar to the
odesolve package for some time, but I never seem to have the week or so
of time I'd need to do it. You can simulate the action of lsodar by
breaking your simulation in short sections, and doing a search to
identify time intervals where a state variable drops below its critical
value (that is, suppose you note that at t1 y[1] 0.4, and at t3, y[1]
0.4]; then search the time interval between t1 and t3 for the value
where abs(y[1] - 0.4) eps, say t2 ). Stop the current integration at
t2, change the state variables, and restart. For any given problem,
you'd probably have to experiment with time reporting intervals (the
intervals between points in the 'times' vector) so as not to miss any
important events.
Woody
On Jun 9, 2004, at 1:43 PM, Martin Henry H. Stevens wrote:
I have a new and less distressing, but potentially more interesting,
problem.
I realized the major flaw my old solution and now have a solution
that kind of works but is rather inelegant and I think may be
problematic in difficult systems.
Borrowing from the lsoda example again I once again highlight the code
that I have changed to put in place arbitrary thresholds:
parms - c(k1=0.04, k2=1e4, k3=3e7)
my.atol - c(1e-6, 1e-10, 1e-6)
times - seq(0,1000)
lsexamp - function(t, y, p)
{
if(y[1] .4) yd1 - -y[1] ### These if, else statements are new
else yd1 - -p[k1] * y[1] + p[k2] * y[2]*y[3]
if(y[3] .4) yd3 - -y[3] ### These if,else statements are new
else yd3 - p[k3] * y[2]^2
list(c(yd1,-yd1-yd3,yd3),c(massbalance=sum(y)))
}
out - lsoda(c(.5,0,.5),times,lsexamp, parms, rtol=1e-4, atol=
my.atol,
hmax=.1)
matplot(out[,1],out[,2:5], type=l)
out[dim(out)[1],] # The intent of my could was to cause population 1
to
fall to zero as soon as it reached 0.4. However, the populations 1
and 2 reach approximations of 0 (4e-281 and 5e-11).
So, I have two questions:
Can I set thresholds in a more elegant and simpler way?
Are the approximate zero values close enough?
Thank you kindly, as ever.
Sincerely,
Hank
On Jun 9, 2004, at 12:45 PM, Martin Henry H. Stevens wrote:
using R 2.0.0
I am trying to do some population modeling with lsoda, where I set
arbitrary zero population sizes when values get close to zero, but am
having no luck.
As an example of what I have tried, I use code below from the help
page on lsoda in which I include my modification bordered by ###
parms - c(k1=0.04, k2=1e4, k3=3e7)
my.atol - c(1e-6, 1e-10, 1e-6)
times - seq(0,)
lsexamp - function(t, y, p)
{ ### The next line is where I try to insert the threshold
ifelse(y 0.4, 0, y)
## all else is unchanged
yd1 - -p[k1] * y[1] + p[k2] * y[2]*y[3]
yd3 - p[k3] * y[2]^2
list(c(yd1,-yd1-yd3,yd3),c(massbalance=sum(y)))
}
out - lsoda(c(.5,0,.5),times,lsexamp, parms, rtol=1e-4, atol=
my.atol) # Initial values differ from help page
matplot(out[,1],out[,2:5], type=l)
out[dim(out)[1],] # The intent of my could was to cause population 1
to fall to zero as soon as it reached 0.4
Any thoughts would be appreciated. Thanks!
Hank Stevens
Dr. Martin Henry H. 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/botany/bot/henry.html
http://www.muohio.edu/ecology/
http://www.muohio.edu/botany/
E Pluribus Unum
R. Woodrow Setzer, Jr.Phone: (919) 541-0128
Experimental Toxicology Division Fax: (919) 541-4284
Pharmacokinetics Branch
NHEERL B143-01; US EPA; RTP, NC 27711
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] ODE solvers in R (was:The Wrong Choice: Locked in by licenserestrictions)
I have plans to add dassl to odesolve. However, I won't have time in the immediate future. If someone wanted to take this on, I'd be happy to give advice and continue to manage the package ... . R. Woodrow Setzer, Jr.Phone: (919) 541-0128 Experimental Toxicology Division Fax: (919) 541-4284 Pharmacokinetics Branch NHEERL B143-05; US EPA; RTP, NC 27711 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] creating function bodies using body()
I'm having trouble figuring out how to create a function using body-
(). The help file for body() says that the argument should be a list of
R expressions. However if I try that I get an error:
tmpfun - function(a, b=2){}
body(tmpfun) - list(expression(z - a + b),expression(z^2))
Error in as.function.default(c(formals(f), value), envir) :
invalid formal argument list for function
Can someone give me a simple example for doing this? Thanks!
R. Woodrow Setzer, Jr.Phone: (919) 541-0128
Experimental Toxicology Division Fax: (919) 541-4284
Pharmacokinetics Branch
NHEERL B143-05; US EPA; RTP, NC 27711
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] question about model formula
I did not know about stepsize being an issue. I had thought that problems with convergence in this case were due to bad approximations of the finite difference gradient. I guessed that around the optimum, numerical errors would come to dominate the gradient calculations, causing convergence to fail. I've found that the issue about confusing optimizers that use gradients can sometimes be fixed by augmenting the original system of odes with what I believe engineers call the sensitivity equations. If your original equation is dg/dt = f(t, b), where b is a parameter to be estimated, then include d^2g/dtdb (be sure to remember the chain rule when doing this!) with the original equation. With some regularity assumptions, this integrates to dg/db, which can be used to give nls a gradient to work with. R. Woodrow Setzer, Jr.Phone: (919) 541-0128 Experimental Toxicology Division Fax: (919) 541-4284 Pharmacokinetics Branch NHEERL B143-05; US EPA; RTP, NC 27711 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Problems with Rcmd check on Win 2000 rw1062
When I run Rcmd check on a package on my Windows 2000 machine, I get a series of error messages like the following: * checking generic/method consistency ...c:\DOCUME~1\R5018~1.WOO\LOCALS~1\Temp/R utils138414013: cannot open c:DOCUME~1R5018~1.WOOLOCALS~1Temp/Rin138408157: no s uch file It looks as if a Windows style path to the temp directory is not being interpreted correctly, with backslashes not being properly escaped. However, I define three different environment variables to point to a temp directory (TMP, TEMP, TMPDIR), and all definitions use forward slashes, Unix-style, so Rcmd check has found some other way to construct the path to a temp directory. Does anyone know how I can fix this so that Rcmd check will work? --please do not edit the information below-- Version: platform = i386-pc-mingw32 arch = i386 os = mingw32 system = i386, mingw32 status = major = 1 minor = 6.2 year = 2003 month = 01 day = 10 language = R Windows 2000 Professional (build 2195) Service Pack 2.0 Search Path: .GlobalEnv, package:OPmodels, package:odesolve, package:ctest, Autoloads, package:base R. Woodrow Setzer, Jr.Phone: (919) 541-0128 Experimental Toxicology Division Fax: (919) 541-4284 Pharmacokinetics Branch NHEERL B143-05; US EPA; RTP, NC 27711 __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Problems with Rcmd check on Win 2000 rw1062
TMPDIR seems to be OK: echo %TMPDIR% C:/DOCUME~1/R5018~1.WOO/LOCALS~1/Temp I'll poke around in the Perl code to see what I can find. R. Woodrow Setzer, Jr.Phone: (919) 541-0128 Experimental Toxicology Division Fax: (919) 541-4284 Pharmacokinetics Branch NHEERL B143-05; US EPA; RTP, NC 27711 [EMAIL PROTECTED] c.uk To: Woodrow Setzer/RTP/USEPA/US@EPA cc: [EMAIL PROTECTED] 02/11/03 02:01 PMSubject: Re: [R] Problems with Rcmd check on Win 2000 rw1062 Rcmd check uses TMPDIR, so what exactly do you have that set to? It does look to me as if you do not have it set correctly, but you can debug the Perl to find out what it is doing. It certainly works for me iff TMPDIR is set correctly. On Tue, 11 Feb 2003 [EMAIL PROTECTED] wrote: When I run Rcmd check on a package on my Windows 2000 machine, I get a series of error messages like the following: * checking generic/method consistency ...c: \DOCUME~1\R5018~1.WOO\LOCALS~1\Temp/R utils138414013: cannot open c:DOCUME~1R5018~1.WOOLOCALS~1Temp/Rin138408157: no s uch file It looks as if a Windows style path to the temp directory is not being interpreted correctly, with backslashes not being properly escaped. However, I define three different environment variables to point to a temp directory (TMP, TEMP, TMPDIR), and all definitions use forward slashes, Unix-style, so Rcmd check has found some other way to construct the path to a temp directory. Does anyone know how I can fix this so that Rcmd check will work? -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Problems with Rcmd check on Win 2000 rw1062
I do not understand why, but changing the definition of TMPDIR to C:/tmp solved the problem. R. Woodrow Setzer, Jr.Phone: (919) 541-0128 Experimental Toxicology Division Fax: (919) 541-4284 Pharmacokinetics Branch NHEERL B143-05; US EPA; RTP, NC 27711 Woodrow Setzer/RTP/USEPA/US@EPA Sent by: [EMAIL PROTECTED] 02/11/2003 02:14 PM To: [EMAIL PROTECTED] cc: [EMAIL PROTECTED] Subject:Re: [R] Problems with Rcmd check on Win 2000 rw1062 TMPDIR seems to be OK: echo %TMPDIR% C:/DOCUME~1/R5018~1.WOO/LOCALS~1/Temp I'll poke around in the Perl code to see what I can find. R. Woodrow Setzer, Jr.Phone: (919) 541-0128 Experimental Toxicology Division Fax: (919) 541-4284 Pharmacokinetics Branch NHEERL B143-05; US EPA; RTP, NC 27711 [EMAIL PROTECTED] c.uk To: Woodrow Setzer/RTP/USEPA/US@EPA cc: [EMAIL PROTECTED] 02/11/03 02:01 PMSubject: Re: [R] Problems with Rcmd check on Win 2000 rw1062 Rcmd check uses TMPDIR, so what exactly do you have that set to? It does look to me as if you do not have it set correctly, but you can debug the Perl to find out what it is doing. It certainly works for me iff TMPDIR is set correctly. On Tue, 11 Feb 2003 [EMAIL PROTECTED] wrote: When I run Rcmd check on a package on my Windows 2000 machine, I get a series of error messages like the following: * checking generic/method consistency ...c: \DOCUME~1\R5018~1.WOO\LOCALS~1\Temp/R utils138414013: cannot open c:DOCUME~1R5018~1.WOOLOCALS~1Temp/Rin138408157: no s uch file It looks as if a Windows style path to the temp directory is not being interpreted correctly, with backslashes not being properly escaped. However, I define three different environment variables to point to a temp directory (TMP, TEMP, TMPDIR), and all definitions use forward slashes, Unix-style, so Rcmd check has found some other way to construct the path to a temp directory. Does anyone know how I can fix this so that Rcmd check will work? -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help [[alternate HTML version deleted]] __ [EMAIL PROTECTED] mailing list http://www.stat.math.ethz.ch/mailman/listinfo/r-help
