I've had pretty good luck solving problems like this when the post
included a simple, reproducible example, as suggested in the posting
guide (www.R-project.org/posting-guide.html). Without that, I'm close
to clueless. I've encountered this kind of error myself before, and I
think the way I solved it was just simplifying the example in steps
until the problem went away. By doing this, I could isolate the source
of the problem.
Have you tried "trace=1", as described in the nlminb help page? With
this, nlminb prints "The value of the objective function and the
parameters [at] every trace'th iteration." Also, have you worked all
the examples in the nlminb help page? If that didn't show me how to
solve my problem, I think I'd next try having nlminb estimate the
gradient and hessian, rather than supplying them myself.
hope this helps,
spencer graves
Greg Tarpinian wrote:
> All,
>
> I have been able to successfully use the optim( ) function with
> "L-BFGS-B" to find reasonable parameters for a one-compartment
> open pharmacokinetic model. My loss function in this case was
> squared error, and I made no assumptions about the distribution
> of the plasma values. The model appeared to fit pretty well.
>
> Out of curiosity, I decided to try to use nlminb( ) applied to
> a likelihood function that assumes the plasma values are normally
> distributed on the semilog scale (ie, modeling log(conc) over
> time). nlminb( ) keeps telling me that it has converged, but
> the estimated parameters are always identical to the initial
> values.... I am certain that I have committed "ein dummheit"
> somewhere in the following code, but not sure what... Any help
> would be greatly appreciated.
>
> Kind regards,
>
> Greg
>
>
>
> model2 <- function(parms, dose, time, log.conc)
> {
> exp1 <- exp(-parms[1]*time)
> exp2 <- exp(-parms[2]*time)
> right.hand <- log(exp1 - exp2)
> numerator <- dose*parms[1]*parms[2]
> denominator <- parms[3]*(parms[2] - parms[1])
> left.hand <- log(numerator/(denominator))
> pred <- left.hand + right.hand
>
> # defining the distribution of the values
> const <- 1/(sqrt(2*pi)*parms[4])
> exponent <- (-1/(2*(parms[4]^2)))*(log.conc - pred)^2
> likelihood <- const*exp(exponent)
>
> #defining the merit function
> -sum(log(likelihood))
> }
>
> deriv2
> <- deriv( expr = ~ -log(1/(sqrt(2*pi)*S)*exp((-1/(2*(S^2)))*
> (log.conc-(log(dose*Ke*Ka/(Cl*(Ka-Ke)))
> +log(exp(-Ke*time)-exp(-Ka*time))))^2)),
> namevec = c("Ke","Ka","Cl","S"),
> function.arg = function(Ke, Ka, Cl, S, dose, time, log.conc) NULL )
>
> gradient2.1compart <- function(parms, dose, time, log.conc)
> {
> Ke <- parms[1]; Ka <- parms[2]; Cl <- parms[3]; S <- parms[4]
> colSums(attr(deriv2.1compart(Ke, Ka, Cl, S, dose, time, log.conc),
> "gradient"))
> }
>
> attach(foo.frame)
> inits <- c(Ke = .5,
> Ka = .5,
> Cl = 1,
> S = 1)
>
> #Trying out the code
> nlminb(start = inits,
> objective = model2,
> gradient = gradient2,
> control = list(eval.max = 5000, iter.max = 5000),
> lower = rep(0,4),
> dose = DOSE,
> time = TIME,
> log.conc = log(RESPONSE))
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
