[R-sig-dyn-mod] Help with function

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.

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]]

--
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531

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]]

--
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531

--
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531

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 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]]

--
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531

--
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531

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]]

--
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531

--
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531