I'm having trouble adding a constraint to my nonlinear model. Suppose I have the following data that is roughly an integrated Gaussian:
x = 1:100
y = pnorm(x, mean = 50, sd = 15) + rnorm(length(x), mean = 0, sd = 0.03)
model <- nls(y ~ pnorm(x, mean = a, sd = b), start = list(a = 50, b = 15))
I can fit the data with nls, but I would like to add the constraint that my fit must fit the data exactly (i.e. have no residual) at y = 0.25 (or whatever point is closest to 0.25). I assume that I need to use glmc for this, but I can't figure out how to use it.
I know it's not necessarily kosher to make the fit adhere to the data like that, but I'm trying to replicate another person's work and this is what they did.
You could impose the restriction somewhat manually. That is, for any parameter b we can solve for a unique a (since the cdf of the normal distribution is strictly increasing) that the restriction would hold:
getA <- function(b, x, y)
optim(x, function(a) (pnorm(x, mean = a, sd = b) - y)^2, method = "BFGS")$par
Then, after finding (tx,ty), the observation of interest, with
idx <- which.min(abs(y - 0.25))
tx <- x[idx]
ty <- y[idx]
we can fit the model with a single parameter:
model <- nls(y ~ pnorm(x, mean = getA(b, tx, ty), sd = b), start = list(b = 15))
and get that the restriction is satisfied
resid(model)[idx]
# [1] -2.440452e-07
and the coefficient a is
getA(coef(model), tx, ty)
# [1] 51.00536