[R-sig-dyn-mod] Help!

[Julian OVIEDO SANTANA]

Hello Prof. Petzoldt


I apologize for writing to you from a personal mail account, please receive this from my academic one. I am sending my problem to the mailing list as well.

I reduced the problem down to just one experiment and also removed the virtual (rnorm()) data I added before. In fact, this is the way the experiment should operate, without the last state variable.
Nonetheless, the sensitivities still change when modifying the sample time.
Could it be something wrong with the model?

Kind regards,

Julian.




# Data --------------------------------------------------------------------


## Import data from file and place them onto a data.frame
Data12 <- data.frame(cbind(A12_full$tmes2, A12_full$Xmes12, A12_full$NO3mes12, A12_full$QNmes12))


## Retrieve parameters from file
pars <- c( '15 parameters’... )


## Constant parameters
L   <- 0.2   ;  Nin <- 25*200  ;  K   <- 14/62000  ;  I   <- c(60, 160, 220, 120, 180, 180)



# Model -------------------------------------------------------------------


rates_batch <- function(t, state, pars) {
  with(as.list(c(state, pars)), {

    X <- state[[1]];
    N <- state[[2]];
    QN <- state[[3]];
    Istar <- state[[4]]

    ## Constants
    D <- 0

    ## Reaction kinetics
    gamma <- gamma_max*(kistar/(kistar + Istar))
    Chl   <- gamma*X*QN
    theta <- Chl/X
    KsI   <- KsIstar/theta
    Ibar  <- Isim*(Katt/(Katt + (a*Chl + b*X + c)*L))
    m_max <- mumaxa*(Istar/(KsI + Istar + (Istar^2/KiI)*exp(Chl*1000*(L/2))))
    mu    <- m_max*(1 - (qnmin/QN))
    rhon  <- rhonmax*(N/(N + kn*K))*(1 - QN/qnmax)

    ## Model RHS
    dx <- mu*X - D*X - R*X
    dn <- -rhon*X + D*(Nin*K - N)
    dq <- rhon - mu*QN
    dI <- delta*mu*(Ibar - Istar)

    return(list(c(dx, dn, dq, dI)))
  })
}



# Simulation --------------------------------------------------------------


## In here we run the model with and remove the last variable “I” from out12
times_2 <- seq(0, 20, by = 0.1)
Isim <- I[4]

model_out12 <- function(pars, times = times_2, Isim) {
  state12 <- c(X = X0_12, N = NO30_12*K, Q = QN0_12, I = Istar0_12)
  out12 <- data.frame(ode(y = state12, times = times, func = rates_batch, parms = pars, method = "lsoda", rtol=1e-8, atol=1e-9))
  out12[,3] <- out12[,3]/K ; out12 <- out12[,-5]
  return(out12)
}
out12 <- model_out12(pars)



# Cost function -----------------------------------------------------------


model_cost <- function(pars) {
  out12 <- model_out12(pars)
  cost <- modCost(model = out12, obs = Data12, weight = "std")
  return(cost)
}
model_cost(pars)$model



# Local sensitivity -------------------------------------------------------


## Sensitivity functions completely different when changing between model_cost and model_out12 in the argument ‘func’…shouldn’t be, right? … also when modifying the sample time (in ‘times_2’), this happens.

sFun <- sensFun(func = model_cost, parms = pars)
plot(sFun)






UMONS
____________________
Julián Oviedo
Ph.D Student
Université de Mons
Automatic Control Laboratory
Boulevard Dolez, 31
7000 Mons
Tel: +32(0)65 37 41 31
E-mail: julian.oviedosantana at umons.ac.be<mailto:julian.oviedosantana at umons.ac.be>
____________________

Automatic Control Laboratory  http://autom.fpms.ac.be<http://autom.fpms.ac.be/>
UMons  http://www.umons.ac.be

Make a change, start now!
https://goodlivingchoices.usana.com



El 20/01/2017, a las 10:09, Thomas Petzoldt <thomas.petzoldt at tu-dresden.de<mailto:thomas.petzoldt at tu-dresden.de>> escribió:

Dear Julian,

I think that your problem are the random numbers generated with rnorm in your example, so that the results differ, depending on the length of the generated data.

Unfortunately, your example is still not a **minimum** reproducible example, so I could not go into the details due to other responsibilities.

If the "rnorm" issue dos not solve your problem, please try the following:

1. reduce your example as much as possible, so that the observed effect is still visible. It is likely that you will find out what's going wrong yourself by doing this.

2. Increase motivation of the people to assist you by disclosing your affiliation.

3. If the problem remains, send it to the mailing list, not only to me.


As said, I think that the rnorm() call is a good candidate to look at.

Regards,

Thomas


	[[alternative HTML version deleted]]