[R-sig-dyn-mod] Help with function
Daniel Kaschek
daniel.kaschek at physik.uni-freiburg.de
Mon Mar 16 10:16:35 CET 2015
Dear Austin,
On So, 2015-03-15 at 12:38 -0600, Austin Mullings wrote:
> For example, below is my equation and a link to my dataset variable names.
> When I run the syntax for that model, it does not run on my data and give
> me a summary of the fit and other results.
I think, this is a misunderstanding. The deSolve-package is supposed to
generate a numerical solution of an ODE. It doesn't do any parameter
estimation, fitting or other things.
If you have data and want to estimate your parameters p[1] to p[10] from
that, you could to set up a function similar to the following
res <- function(p) {
prediction <- ode(y, times, func, p)
residuals <- as.vector(prediction[,-1] - data)
return(residuals)
}
and use it with nls(), i.e.
myfit <- nls(~res, start = list(p = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1)))
The data used in the function res() must have the same structure as the
output of ode(), i.e. one column per state, ordered in the same way as
the states and the time points must coincide between data and
prediction.
Good luck!
Best,
Daniel.
> I was wondering if you could
> help me learn to run it, and after this example, I would hopefully be able
> to run these models on my own.
>
> Link:
> http://s1074.photobucket.com/user/austinmullings89/media/Picture%20of%20Levy%20Model%20data.jpg.html
>
> Equation syntax:
>
> library(deSolve)
>
> myfn <- function(t, y, p) {
>
> dy <- numeric(6)
> dy[1] <- p[1] - p[2]*y[4] - p[3]*y[1]
> dy[2] <- y[5]*dy[1]*(1+y[6]) - p[5]*y[2] - (p[6]*y[3]*y[4]*y[1] - p[7])
> dy[3] <- p[8] - p[9]*y[2]
> dy[4] <- y[6]*dy[1] - p[10]*y[4] + p[11]*(p[6]*y[3]*y[4]*y[1] - p[7]) +
> p[4]
> dy[5] <- 1.0 - y[5]*(y[1] + y[4] + 1)
> dy[6] <- 1.0 - y[6]*(y[1] + 1)
>
> return(list(dy))
>
> }
>
> pars <- runif(11)
> yini <- runif(6)
> times <- seq(0, 10, .01)
>
>
> out <- ode(yini, times, myfn, pars)
> plot(out)
>
> On Thu, Mar 12, 2015 at 8:23 AM, Daniel Kaschek <
> daniel.kaschek at physik.uni-freiburg.de> wrote:
>
> > Hi Austin,
> >
> > On Mi, 2015-03-11 at 10:03 -0600, Austin Mullings wrote:
> > > I had one more question, which may be a novice question, but how do I
> > need
> > > to create my dataset in a way that this function will run off of it?
> > > Specifically, I am trying to run it on some mock data (until I learn more
> > > about simulation), and do not know how R and this package expects my
> > > function or dataset column/variable names to be for it to run. Here is a
> > > picture of my mock data, if you wouldn't mind letting me know whether I
> > > have the variable names ordered incorrectly or need to rename them
> > > accordingly.
> >
> > I think you cannot send pictures via the mailing list. Unfortunately, I
> > didn't get your question. What do you mean by "data to run your
> > simulation on"?
> >
> > Concerning the sorting of variables: In the example below, I
> > initialized
> >
> > yini <- runif(6)
> >
> > resulting in a vector yini of length 6. The i'th element of yini is y[i]
> > in myfn (at time point 0). This means that the i'th element is the
> > initial value of the i'th state, which in turn will be the i'th column
> > next to the time column of the ode() output.
> >
> > Cheers,
> > Daniel
> >
> > >
> > >
> > >
> > >
> > > On Tue, Feb 17, 2015 at 2:04 AM, Daniel Kaschek <
> > > daniel.kaschek at physik.uni-freiburg.de> wrote:
> > >
> > > > Hi Austin,
> > > >
> > > > one possibility to do this is the code below. In this example code, I
> > > > randomly initialize your parameters and the initial state. Instead of
> > > > using the indexes to access parameter values or states you can use the
> > > > names as well (see ?ode, first example). However, in my experience, the
> > > > function with() slows down your code considerably.
> > > >
> > > > Best wishes,
> > > > Daniel
> > > >
> > > >
> > > > library(deSolve)
> > > >
> > > > myfn <- function(t, y, p) {
> > > >
> > > > dy <- numeric(6)
> > > > dy[1] <- p[1] - p[2]*y[4] - p[3]*y[1]
> > > > dy[2] <- y[5]*dy[1]*(1+y[6]) - p[5]*y[2] - (p[6]*y[3]*y[4]*y[1] -
> > p[7])
> > > > dy[3] <- p[8] - p[9]*y[2]
> > > > dy[4] <- y[6]*dy[1] - p[10]*y[4] + p[11]*(p[6]*y[3]*y[4]*y[1] -
> > p[7]) +
> > > > p[4]
> > > > dy[5] <- 1.0 - y[5]*(y[1] + y[4] + 1)
> > > > dy[6] <- 1.0 - y[6]*(y[1] + 1)
> > > >
> > > > return(list(dy))
> > > >
> > > > }
> > > >
> > > > pars <- runif(11)
> > > > yini <- runif(6)
> > > > times <- seq(0, 10, .01)
> > > >
> > > >
> > > > out <- ode(yini, times, myfn, pars)
> > > > plot(out)
> > > >
> > > >
> > > > On Mo, 2015-02-16 at 13:15 -0700, Austin Mullings wrote:
> > > > > I am trying to get this equation to work with deSolve, and can't
> > seem to
> > > > > get it to work. Does someone know how one would write the function
> > for
> > > > this?
> > > > >
> > > > > dY1(t)/dt = a - bY3(t) Y4(t) - cY1(t),
> > > > > dY2(t)/dt = Y5(t) [dY1(t)/dt] [1 + Y6(t)] - eY2(t) - {f Y3(t) Y4(t)
> > > > Y1(t) -
> > > > > g},
> > > > > dY3(t)/dt = h - iY2(t),
> > > > > dY4(t)/dt = Y6(t) [dY1(t)/dt] - jY4(t) + k {f Y3(t) Y4(t) Y1(t) - g}
> > + d,
> > > > > dY5(t)/dt = 1.0 - Y5(t) [Y1(t) + Y4(t) + 1],
> > > > > dY6(t)/dt = 1.0 - Y6(t) [Y1(t) + 1].
> > > > >
> > > > > [[alternative HTML version deleted]]
> > > > >
> > > > > _______________________________________________
> > > > > R-sig-dynamic-models mailing list
> > > > > R-sig-dynamic-models at r-project.org
> > > > > https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> > > >
> > > > --
> > > > Daniel Kaschek
> > > > Institute of Physics
> > > > Freiburg University
> > > >
> > > > Room 210
> > > > Phone: +49 761 2038531
> > > >
> > > > _______________________________________________
> > > > R-sig-dynamic-models mailing list
> > > > R-sig-dynamic-models at r-project.org
> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> > > >
> > >
> > >
> > >
> >
> > --
> > Daniel Kaschek
> > Institute of Physics
> > Freiburg University
> >
> > Room 210
> > Phone: +49 761 2038531
> >
> > _______________________________________________
> > R-sig-dynamic-models mailing list
> > R-sig-dynamic-models at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> >
>
>
>
More information about the R-sig-dynamic-models
mailing list