1

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.

Gabe
  • 1,722
  • 1
  • 11
  • 7

1 Answers1

0

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
Julius Vainora
  • 47,421
  • 9
  • 90
  • 102
  • So for each parameter `b` the constraint means that there is exactly one pnorm line that goes through the proper point. Thus, the constraint effectively reduces the model down to a 1-parameter model. Good observation. – Gabe Mar 28 '15 at 05:21
  • @Gabe, yes, since the cdf or normal distribution is strictly increasing, we have exactly one solution. – Julius Vainora Mar 28 '15 at 11:22