# How to fit an implicit function in R?

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")

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:

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.