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

Tue Mar 17 20:22:02 CET 2015

Yes nls can be used for parameter fitting, and there are also other 
packages, e.g. package FME and the following paper:

http://www.jstatsoft.org/v33/i03

or the following teaching example

http://desolve.r-forge.r-project.org/user2014/examples/FME/fit_twocomp.svg

and please note that "Fitting models is an art ..."

Thomas

Am 3/16/2015 um 10:16 AM schrieb Daniel Kaschek:
> 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