4
$\begingroup$

I need to fit experimental data $x$ (dependent variable, with random error in log scale) vs $t \ge 0$ (independent variable, with ~ no error), where however $x(t)$ is not analytically known, because:

$$t = \frac b k \cdot \lbrack (a \cdot (b-1) -1) \cdot \ln {(1-\frac {1-x} b) + a \cdot (1-x)}\rbrack$$

which, as far as I can tell, is not invertible.

Note that $a > 0, 0 < b < 1, k > 0$.
From the above function and constraints it also follows that $1-b < x \le 1$.

My goal is to find the values of $a, b, k$ s.t. $\sum_i (\ln x_{exp,i}- \ln x_{pred,i})^2$ is minimal.

See below my attempt, R code.
I don't particularly like the result. Parameter $a$ is really badly fitted, I don't know why, and despite the fact that I used logs for the residuals, the points at large values of $t$ are really badly fitted, too.

I don't understand if I am doing something wrong but using the correct method, or if I am using a wrong method in the first place.
In particular, I am not sure if it is legitimate to use a numerical zero finder (uniroot) to invert $t(x)$.

Some posts I saw on the web say for such cases one should do a 'constrained optimization with nonlinear constraints', where $g(x,t) = 0$ is the constraint.
See for instance https://stat.ethz.ch/pipermail/r-help/2010-October/255116.html.
But I never tried that, I would not know where to start.

Any suggestions / ideas / pointers to relevant posts or literature?

The R code:

# Define the true parameters

a = 0.1 # a > 0, preferably a < 1
b = 0.7 # 0 < b < 1
k = 0.1 # k > 0

# Define a function that finds x from t and the parameters

xt <- function(t, x_start = 1, parms) {
  with(as.list(parms), {
  uniroot(function(x) t - b / k * ((a * (b - 1) - 1) * 
     log(1 - (1 - x) / b) + a * (1 - x)), c(1 - b, 1), 
      tol = 1.e-6)$root
  })
}

# Plot the true function

times <- seq(0,60,1)
plot(times, sapply(times, xt, parms = c(a = a, b = b, k = k)), 
     type = "l", log = "y", ylab = "x", xlab = "t")

# Simulate and plot the experimental x(t) data

#times_exp <- c(0, 5, 15, 45)
times_exp <- seq(0, 60, 5)
x_true <- sapply(times_exp, xt, parms = c(a = a, b = b, k = k))
set.seed(918234)
x_exp <- exp(log(x_true) + rnorm(length(times_exp), sd = 0.2))
points(times_exp, x_exp, pch = 16, col = "red")

# Fit the experimental data

require(FME)

cost <- function(ps) {
  a <- ps[1]
  b <- ps[2]
  k <- ps[3]
  x_fit <- sapply(times_exp, xt, parms = c(a = a, b = b, k = k))
  log(x_fit) - log(x_exp)
}

b_lower <- 1 - min(x_exp)
k_start <- -1 * coef(lm(log(x_exp) ~ times_exp))[[2]]

mod_out <- modFit(f = cost, p = c(1, 0.1 + b_lower, k_start),  
                  lower = c(0, b_lower, 0), 
                  upper = c(Inf, 1, Inf))

# Plot the fitted function

fitted_parms <- setNames(mod_out$par, c("a","b","k"))
lines(times, sapply(times, xt, parms = fitted_parms), 
      lty = 2, col = "blue")

enter image description here


EDIT after removing the lower bound for parameter $b$

In fact, I found that it's wrong to put a lower bound on $b$.
When I change the fitting part of the code to:

p_start <- c(1, 0.5, k_start)
mod_out <- modFit(f = cost, p = p_start, lower = c(0, 0, 0), 
                  upper = c(Inf, 1, Inf))

I get a much better picture:

enter image description here

There is still some numerical instability in some cases, though, due to uniroot finding NA at the extremes. The starting conditions are not always good. I seem to gather that this is a recurring issue with nonlinear fitting.

I'm still interested in possible alternatives or improvements, if anyone has got suggestions.

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.