language design for ODE solvers

[Bob Carpenter]

I updated the doc with some proposals if anyone has
any comments:

https://github.com/stan-dev/stan/wiki/ODE-Integrator-Support

Proposal 2 is by far the easiest to code, but I think Proposal 3 is
the most natural. I had originally thought we'd go with Proposal 1,
but I also think that'll be both too tricky to code and for users to
call and also too redundant.

- Bob

[Ben Goodrich]

In addition to the PKBugs applications, we should think about what is the best way to proceed when the density is defined in terms of an ODE, such as

http://web.udl.es/Biomath/Group/Treballs/2006CSDA.pdf

Ben

[Bob Carpenter]

I added the reference to the design doc, but I'm kind of lost
at formula (1) in the paper, which is

dy/dx = y * (mu - x) / (a + b * x + c * x^2)

I assume a, b, c and mu are parameters, because they describe
them as such. Is x the variate and y the density?

If so, I don't see how we can use an integrator for initial
value problems to solve for the density at a given x, mu,
a, b, and c. But then I'm so bad at diff eqs it hurts to
even think about this problem.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

[Ben Goodrich]

On Wednesday, January 29, 2014 8:42:41 PM UTC-5, Bob Carpenter wrote:

I added the reference to the design doc, but I'm kind of lost
at formula (1) in the paper, which is

   dy/dx = y * (mu - x) / (a + b * x + c * x^2)

I assume a, b, c and mu are parameters, because they describe
them as such.  Is x the variate and y the density?

Yes, for the Pearson system. But it has explicit solutions that we could just implement.

For the GS distribution that the paper advocates, the algorithm is described in section 3.1.

Ben

[Bill Harris]

Thanks, Bob. I admit to liking things that are close to MCSim; whether that's because I'm more familiar with it or because it's better, I'm not sure.

You are creating a simulator that can work for system dynamics modeling, too. As you ask for libraries of models, http://www.metasd.com/models/ is a collection of system dynamics models. Some of the Vensim models are just text, so you can examine those. Others may require the Vensim model reader from http://www.vensim.com/ or the iThink model player from http://www.iseesystems.com/, both of which are free. Follow the direction on Tom Fiddaman's page.

A system dynamics simulator pretty much requires the ability to specify a tabular nonlinear function with some sort of interpolation between defined points. I know Stan is not compatible with the GPL, but I use GSL Akima interpolation routines in MCSim. I'm confident Boost has something similar.

Finally, you might be interested in XMILE (http://xmile.systemdynamics.org/ and https://www.oasis-open.org/committees/tc_home.php?wg_abbrev=xmile), the XML Interchange Language for System Dynamics.

If you're not sure what system dynamics is, see http://www.systemdynamics.org/what-is-s/. Despite them using DE notation, too, most system dynamics models are written as integral equations, something like

state = integrate(netflow, initialconditions)

and better modeling systems allow for the use of units and units checking with suitable units algebra checking.

Bill

[Bob Carpenter]

Thanks for the thoughtful note! I added it to our design doc
along with Ben's notes on definitions of densities
through diffeqs. As much as I find this all a bit overwhelming,
it's probably better to have a heads-up now rather than later.

- Bob

On 1/30/14, 7:20 AM, Bill Harris wrote:
> http://www.metasd.com/models/ is a collection of system dynamics models. Some of the Vensim models are just text, so you can examine those. Others may require the Vensim model reader fromhttp://www.vensim.com/ or the iThink model player fromhttp://www.iseesystems.com/, both of which are free. Follow the direction on Tom Fiddaman's page.

>
> A system dynamics simulator pretty much requires the ability to specify a tabular nonlinear function with some sort of interpolation between defined points. I know Stan is not compatible with the GPL, but I use GSL Akima interpolation routines in MCSim. I'm confident Boost has something similar.
>

> Finally, you might be interested in XMILE (http://xmile.systemdynamics.org/ andhttps://www.oasis-open.org/committees/tc_home.php?wg_abbrev=xmile), the XML Interchange Language for System Dynamics.
>
> If you're not sure what system dynamics is, seehttp://www.systemdynamics.org/what-is-s/. Despite them using DE notation, too, most system dynamics models are written as integral equations, something like

[Michael Betancourt]

I have opinions!  :-)  Oh, and for future reference it's the "simple harmonic oscillator" or "sho".  "ho" just looks weird.  :-p

I immediately don't like (3) because I think it will lead to no end of errors from scope confusing and naming conflicts.  I am a big fan of more strict syntax that prevents these kinds of collisions.

So that leaves (1) and (2).  I like the simplicity of (2)'s syntax (I'm already fearing the ridiculous signatures that would arise from user models if a fully generic (1) were implemented as is) but can see the huge restriction of limiting our users to reals.

That said, is there a way to combine them?  I'm thinking of having data, parameter, and state classes which are just polymorphic containers for Stan variables.  So when you define

differential equations {
  harmonic_oscillator {
    ...
    parameters {
      real<lower=0> gamma;
    }
    …
}

the user would have access to the new type

harmonic_oscillator::parameters

which they could use as

harmonic_oscillator::parameters::gamma <- 5;

Then we get the simple signature of (2) but the user has the freedom of defining different types while also having a clear understanding of their scope.

- Bob

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Bob Carpenter]

That's similar to some ideas Andrew had, and my thought is
that it would be best to do by adding a new list-like type
to Stan.

I'm also not keen on (3), hence my blog post on Andrew's blog
about R's behavior.

- Bob

On 1/30/14, 12:53 PM, Michael Betancourt wrote:
> I have opinions! :-) Oh, and for future reference it's the "simple harmonic oscillator" or "sho". "ho" just looks
> weird. :-p
>
> I immediately don't like (3) because I think it will lead to no end of errors from scope confusing and naming conflicts.
> I am a big fan of more strict syntax that prevents these kinds of collisions.
>
> So that leaves (1) and (2). I like the simplicity of (2)'s syntax (I'm already fearing the ridiculous signatures that
> would arise from user models if a fully generic (1) were implemented as is) but can see the huge restriction of limiting
> our users to reals.
> That said, is there a way to combine them? I'm thinking of having data, parameter, and state classes which are just
> polymorphic containers for Stan variables. So when you define
>
> differential equations {
> harmonic_oscillator {
> ...
> parameters {
> real<lower=0> gamma;
> }
> …
> }
>
> the user would have access to the new type
>
> harmonic_oscillator::parameters
>
> which they could use as
>
> harmonic_oscillator::parameters::gamma <- 5;
>
> Then we get the simple signature of (2) but the user has the freedom of defining different types while also having a
> clear understanding of their scope.
>
>
>
> On Wed, Jan 29, 2014 at 8:20 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
>
> I updated the doc with some proposals if anyone has
> any comments:
>

> https://github.com/stan-dev/__stan/wiki/ODE-Integrator-__Support

> <https://github.com/stan-dev/stan/wiki/ODE-Integrator-Support>
>
> Proposal 2 is by far the easiest to code, but I think Proposal 3 is
> the most natural. I had originally thought we'd go with Proposal 1,
> but I also think that'll be both too tricky to code and for users to
> call and also too redundant.
>
> - Bob
>
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to

> stan-dev+unsubscribe@__googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

>
>
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Michael Betancourt]

It's exciting -- I can already see building up some cool models with a syntax like that.

On Thu, Jan 30, 2014 at 12:00 PM, Bob Carpenter <ca...@alias-i.com> wrote:

That's similar to some ideas Andrew had, and my thought is
that it would be best to do by adding a new list-like type
to Stan.

I'm also not keen on (3), hence my blog post on Andrew's blog
about R's behavior.

- Bob

On 1/30/14, 12:53 PM, Michael Betancourt wrote:

I have opinions!  :-)  Oh, and for future reference it's the "simple harmonic oscillator" or "sho".  "ho" just looks
weird.  :-p

I immediately don't like (3) because I think it will lead to no end of errors from scope confusing and naming conflicts.
  I am a big fan of more strict syntax that prevents these kinds of collisions.

So that leaves (1) and (2).  I like the simplicity of (2)'s syntax (I'm already fearing the ridiculous signatures that
would arise from user models if a fully generic (1) were implemented as is) but can see the huge restriction of limiting
our users to reals.
That said, is there a way to combine them?  I'm thinking of having data, parameter, and state classes which are just
polymorphic containers for Stan variables.  So when you define

differential equations {
   harmonic_oscillator {
     ...
     parameters {
       real<lower=0> gamma;
     }

     ...

}

the user would have access to the new type

harmonic_oscillator::parameters

which they could use as

harmonic_oscillator::parameters::gamma <- 5;

Then we get the simple signature of (2) but the user has the freedom of defining different types while also having a
clear understanding of their scope.

On Wed, Jan 29, 2014 at 8:20 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:

    I updated the doc with some proposals if anyone has
    any comments:

    https://github.com/stan-dev/__stan/wiki/ODE-Integrator-__Support

    <https://github.com/stan-dev/stan/wiki/ODE-Integrator-Support>

    Proposal 2 is by far the easiest to code, but I think Proposal 3 is
    the most natural.  I had originally thought we'd go with Proposal 1,
    but I also think that'll be both too tricky to code and for users to
    call and also too redundant.

    - Bob

    --
    You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
    To unsubscribe from this group and stop receiving emails from it, send an email to

    stan-dev+unsubscribe@__googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.
    For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Marcus Brubaker]

I also think (3) is a bad idea.  I think those kinds of weird scoping rules just lead to unexpected surprises.  Plus, so far Stan has been pretty strict in it's language (e.g., requiring explicit declaration of variables) so (1) or (2) definitely seems more "natural" in that sense.

Personally, i'd vote for (2) for a two reasons.  First, it's a familiar style of interface for ODE solvers so I think most people will be comfortable with it.  (1) is more expressive, but is going to be tricky to get right, e.g., what happens if a user defines constrained parameters, how are initial conditions passed in when there are multiple parameters, how is the result returned?  Second, I think in the process of implementing (2) it would be easy to at least lay the groundwork for more general user-defined functions in Stan which is a feature that's been on the wish-list for a long time and would help clean-up more complex models.

Cheers,

Marcus

On Thu, Jan 30, 2014 at 7:00 AM, Bob Carpenter <ca...@alias-i.com> wrote:

That's similar to some ideas Andrew had, and my thought is
that it would be best to do by adding a new list-like type
to Stan.

I'm also not keen on (3), hence my blog post on Andrew's blog
about R's behavior.

- Bob

On 1/30/14, 12:53 PM, Michael Betancourt wrote:

I have opinions!  :-)  Oh, and for future reference it's the "simple harmonic oscillator" or "sho".  "ho" just looks
weird.  :-p

I immediately don't like (3) because I think it will lead to no end of errors from scope confusing and naming conflicts.
  I am a big fan of more strict syntax that prevents these kinds of collisions.

So that leaves (1) and (2).  I like the simplicity of (2)'s syntax (I'm already fearing the ridiculous signatures that
would arise from user models if a fully generic (1) were implemented as is) but can see the huge restriction of limiting
our users to reals.
That said, is there a way to combine them?  I'm thinking of having data, parameter, and state classes which are just
polymorphic containers for Stan variables.  So when you define

differential equations {
   harmonic_oscillator {
     ...
     parameters {
       real<lower=0> gamma;
     }

     ...

}

the user would have access to the new type

harmonic_oscillator::parameters

which they could use as

harmonic_oscillator::parameters::gamma <- 5;

Then we get the simple signature of (2) but the user has the freedom of defining different types while also having a
clear understanding of their scope.

On Wed, Jan 29, 2014 at 8:20 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:

    I updated the doc with some proposals if anyone has
    any comments:

    https://github.com/stan-dev/__stan/wiki/ODE-Integrator-__Support

    <https://github.com/stan-dev/stan/wiki/ODE-Integrator-Support>

    Proposal 2 is by far the easiest to code, but I think Proposal 3 is
    the most natural.  I had originally thought we'd go with Proposal 1,
    but I also think that'll be both too tricky to code and for users to
    call and also too redundant.

    - Bob

    --
    You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
    To unsubscribe from this group and stop receiving emails from it, send an email to

    stan-dev+unsubscribe@__googlegroups.com <mailto:stan-dev%2Bunsu...@googlegroups.com>.
    For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Bob Carpenter]

On 1/30/14, 6:08 PM, Marcus Brubaker wrote:
> I also think (3) is a bad idea. I think those kinds of weird scoping rules just lead to unexpected surprises. Plus, so
> far Stan has been pretty strict in it's language (e.g., requiring explicit declaration of variables) so (1) or (2)
> definitely seems more "natural" in that sense.
>
> Personally, i'd vote for (2) for a two reasons. First, it's a familiar style of interface for ODE solvers so I think
> most people will be comfortable with it. (1) is more expressive, but is going to be tricky to get right, e.g., what
> happens if a user defines constrained parameters, how are initial conditions passed in when there are multiple
> parameters, how is the result returned? Second, I think in the process of implementing (2) it would be easy to at least
> lay the groundwork for more general user-defined functions in Stan which is a feature that's been on the wish-list for a
> long time and would help clean-up more complex models.

OK, I think everyone's in line with (2).

Of course, that system makes it seem easy. When I was running it
by Andrew, I was using the system from a consulting project we're doing
that has around 10 parameters and marshaing/unmarshaling the vectors
is much more of a pain.

The plan is definitely to implement general functions first.

- Bob

[Marcus Brubaker]

One way to ease this is to allow an arbitrary number of parameters to be passed in.  The state vector would still need to packed and unpacked but it would reduce the complexity otherwise.  For instance, you could do:

vector[] harmonic_oscillator(real t, vector[] state, real param1, vector[] param2, matrix[] param3) {

...

}

and then when calling

solutions <- ode_solve(harmonic_oscillator,

                       intial_state,

                       solution_times, 

                       param1, param2, param3);

Basically, allow the function to take an arbitrary number of parameters of any type and when ode_solve is called, anything after the initial state and solutions times gets passed in to the function.

Cheers,

Marcus

[Bob Carpenter]

On 1/30/14, 7:07 PM, Marcus Brubaker wrote:
...

> One way to ease this is to allow an arbitrary number of parameters to be passed in. The state vector would still need
> to packed and unpacked but it would reduce the complexity otherwise. For instance, you could do:
>
> vector[] harmonic_oscillator(real t, vector[] state, real param1, vector[] param2, matrix[] param3) {
> ...
> }
>
> and then when calling
>
> solutions <- ode_solve(harmonic_oscillator,
> intial_state,
> solution_times,
> param1, param2, param3);
>
> Basically, allow the function to take an arbitrary number of parameters of any type and when ode_solve is called,
> anything after the initial state and solutions times gets passed in to the function.

The complexity of this is that we need to distinguish data from
parameters for the ODE case. We need the gradient of solutions
w.r.t. the parameters, but not the data.

- Bob

[Marcus Brubaker]

Ah, of course.  I guess we can deduce things at compile time using template magic, but it does increase the complexity.

Cheers,

Marcus

 

- Bob

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

[Bob Carpenter]

On 1/30/14, 8:59 PM, Marcus Brubaker wrote:

> On Thu, Jan 30, 2014 at 1:11 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
>
>
>
> On 1/30/14, 7:07 PM, Marcus Brubaker wrote:
> ...
>
> One way to ease this is to allow an arbitrary number of parameters to be passed in. The state vector would
> still need
> to packed and unpacked but it would reduce the complexity otherwise. For instance, you could do:
>
> vector[] harmonic_oscillator(real t, vector[] state, real param1, vector[] param2, matrix[] param3) {
> ...
> }
>
> and then when calling
>
> solutions <- ode_solve(harmonic_oscillator,
> intial_state,
> solution_times,
> param1, param2, param3);
>
> Basically, allow the function to take an arbitrary number of parameters of any type and when ode_solve is called,
> anything after the initial state and solutions times gets passed in to the function.
>
>
> The complexity of this is that we need to distinguish data from
> parameters for the ODE case. We need the gradient of solutions
> w.r.t. the parameters, but not the data.
>
>
> Ah, of course. I guess we can deduce things at compile time using template magic, but it does increase the complexity.

I was more thinking of doing it at code-generation time, because
we'll have the type of the expression being used in the call. Specifically,
I was thinking of something like this as being massaged into a call
more like suggestion 2 so that it's easy to code in general from the ODE
and auto-diff side.

Templates would probably be hairy for this, but I can think about it.

- Bob

[Bill Harris]

Bob Carpenter <ca...@alias-i.com> writes:

> Thanks for the thoughtful note! I added it to our design doc
> along with Ben's notes on definitions of densities
> through diffeqs. As much as I find this all a bit overwhelming,
> it's probably better to have a heads-up now rather than later.

Bob,

Two more thoughts.

First, you keep writing of Jacobians. I think most system dynamicists
wouldn't know what to do about writing a Jacobian of their model, but
it's been so long since I've had to do anything practical with them that
I'll have to do some study.

I'm presuming that Jacobians are needed for HMC; I don't recall anything
about Jacobians in MCSim.

To the degree that the calculation of Jacobians can be done
automatically, that's great. To the degree they need to be coded, I
suspect a significant number of your potential users may need clear
help.

Second, I recall the biggest problem I have with MCSim and, for that
matter, to a degree with some commercial simulators: the ability to
encapsulate submodels.

For example, say I create a system dynamics (ODE) model of a supply
chain, and I have a planning function as part of that model. The supply
chain is fairly clear to code and diagram, but the planning function
complicates the picture for the reader. isee systems' iThink used to
have a DFD (Decision Flow Diamond) that would hide parts of models. In
this case, I could put the entire planning function inside a DFD and
expand it only when I needed to see it. Otherwise, I'd just look at the
main supply chain and realize there was a planning function that was
intended to take care of keeping it filled.

In a non-graphical system such as MCSim or Stan, it might be that you
can encapsulate dynamic bits (some of the ODEs) in a function that gets
called by the higher-level model.

I'll have to look tomorrow: I presume that state variables and their
derivatives can be non-scalar, right? MCSim has a decent array
functionality, as do most professional SD simulators. That's quite
handy when you have the same structure repeated but perhaps with
different paramater values.

Bill
--
Bill Harris
Facilitated Systems
http://makingsense.facilitatedsystems.com/

[Bob Carpenter]

On 1/31/14, 7:34 AM, Bill Harris wrote:
> Bob Carpenter <ca...@alias-i.com> writes:
>
>> Thanks for the thoughtful note! I added it to our design doc
>> along with Ben's notes on definitions of densities
>> through diffeqs. As much as I find this all a bit overwhelming,
>> it's probably better to have a heads-up now rather than later.
>
> Bob,
>
> Two more thoughts.
>
> First, you keep writing of Jacobians. I think most system dynamicists
> wouldn't know what to do about writing a Jacobian of their model, but
> it's been so long since I've had to do anything practical with them that
> I'll have to do some study.

What I mean here is that we need the Jacobian of the system.
And by system here, I just mean the function that takes the state
and provides the derivatives w.r.t. time of the state at a given
point. So if you have x[1] and x[2] as the state (inputs), then the
diff eq system defines y[1] = d/dt x[1] and y[2] = d/dt x[2] (outputs).
So the system is a function f(x) = y, and what we need is the Jacobian
of the function, namely dy / dx, which is just dy[n] / dx[m] for n in 1:2
and m in 1:2.

They're pretty standard in solving diff eqs. For instance, both Boost
and dSolve in R require user-supplied Jacobians to solve stiff systems.

> I'm presuming that Jacobians are needed for HMC; I don't recall anything
> about Jacobians in MCSim.

No, we only need gradients of the solutions w.r.t. the parameters
for HMC. So if that function f(x) = y depends on some parameters,
that is, f(x) = y depends on alpha and beta, then we need the
gradients of y w.r.t. alpha and beta.

> To the degree that the calculation of Jacobians can be done
> automatically, that's great.

We absolutely want to compute them automatically, so that users
will only have to specify the function f(x) = y in a way that
may involve parameters (as well as constants).

> To the degree they need to be coded, I
> suspect a significant number of your potential users may need clear
> help.
>
> Second, I recall the biggest problem I have with MCSim and, for that
> matter, to a degree with some commercial simulators: the ability to
> encapsulate submodels.
>
> For example, say I create a system dynamics (ODE) model of a supply
> chain, and I have a planning function as part of that model. The supply
> chain is fairly clear to code and diagram, but the planning function
> complicates the picture for the reader. isee systems' iThink used to
> have a DFD (Decision Flow Diamond) that would hide parts of models. In
> this case, I could put the entire planning function inside a DFD and
> expand it only when I needed to see it. Otherwise, I'd just look at the
> main supply chain and realize there was a planning function that was
> intended to take care of keeping it filled.

> In a non-graphical system such as MCSim or Stan, it might be that you
> can encapsulate dynamic bits (some of the ODEs) in a function that gets
> called by the higher-level model.

That's the plan. Unless we get very directed funding or volunteers,
it's unlikely we're going to do anything specific for any of these
model classes.

> I'll have to look tomorrow: I presume that state variables and their
> derivatives can be non-scalar, right? MCSim has a decent array
> functionality, as do most professional SD simulators. That's quite
> handy when you have the same structure repeated but perhaps with
> different paramater values.

The plan for now is to have the system be arrays. But you could
pack anything into those arrays. This is what MCSim and Boost's
solvers do. For instance, MCSim doesn't have any matrix or 2D
array functionality. Part of what we're trying to think about is
how critical that'll be. It probably won't go in our first implementation
even if we do think it's critical. We just don't have much time to
work on this subsystem of Stan.

- Bob

[Bill Harris]

Bob,

Okay; I think I see. Simulators I have used essentially encode the
Jacobian in a set of equations that essentially describe each element
of the Jacobian in one of the equations. If an element isn't
described by some such equation, then the entry is de factor 0 for
that cell of the Jacobian, right?

SD modelers and perhaps others are taught to think structurally about
the problem, so that encoding is nice; a separate Jacobian matrix
definition would not be seen as a transparent model. (Most SD models
have associated "stock and flow" or "influence diagram" graphical
representations that show how stocks / state variables and flows /
derivatives connect together, something loosely like PBPK
compartmental model diagrams

Moreover, each element of a Jacobian / definition of the inflows and
outflows to or from a stock / derivative of (net inflow to) a stock
can be a complex, tabularly-defined nonlinear function, which might be
awkward to put into a matrix.

So I think we're saying much the same thing (and writing about
Jacobians just before going to sleep as I did last night is risky;
obviously that's what a Jacobian is). If you lead with Jacobians, I
suspect some will be scared off.

On a second topic, I don't see any need to do predefined complex
submodels up front or perhaps ever. There is an advantage to
providing an easy-to-use, expressive language and then making people
think through what they're doing (system dynamics languages, Stan,
MCSim, and R all fit that mold) There is a different, productivity
advantage to a language heavily supported by libraries (packages) such
as R. Vensim does provide such a thing for Vensim models; see
http://vensim.com/modeling-with-molecules-2-02/.

Oh, and have you seen http://en.wikipedia.org/wiki/SBML? MCSim works
with SBML, I think, but I've never used it.

Finally, http://en.wikipedia.org/wiki/DYNAMO_(programming_language)
speaks to the original SD language. Several books listed in the
bibliography describe that language, and section 3.7 of
http://eprints.cs.vt.edu/archive/00000363/01/TR-93-21.pdf shows an
example. Even though it's over 50 years old, it still has a strong
influence on SD modeling languages. FYI / FWIW.

Bill

> --
> You received this message because you are subscribed to a topic in the

> Google Groups "stan development mailing list" group.

> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/stan-dev/DWItD1i8_1Q/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> stan-dev+u...@googlegroups.com.

> For more options, visit https://groups.google.com/groups/opt_out.

--
Bill Harris
http://facilitatedsystems.com/weblog/

[Bob Carpenter]

I think just talking about diff-eqs is going to scare
most people off. I didn't even know what a Jacobian
was (or a Hessian) three years ago :-)

I haven't seen SBML. I'm not a huge fan of markup languages.
I know Frederic (designer of MCSim) is using some markup
languages for systems bio diff eqs because he was talking about them
the other day.

I was only being pedantic about the Jacobians because there
are actually a couple of different Jacobians in play. One is
w.r.t. the system inputs and one w.r.t. the parameters. They
serve different purposes.

We'll just lead with the language, as usual.

- Bob

[Bill Harris]

Bob Carpenter <ca...@alias-i.com> writes:

> I think just talking about diff-eqs is going to scare
> most people off.

That's a big reason why SD modelers almost always write integral
equations; they somehow seem more intuitive, and they seem more directly
related to what's going on in the world (I'm adding something here to
this stock [state variable], and I'm taking something away over there.).