Is least squares the standard method to fit a 3 parameters Gaussian function to some x and y data?

Question

A participant in one experiment needs to decide whether a flash and a sound are simultaneous or not for many possible asynchronies between the flash and the sound (x in seconds). For each asynchrony, the flash and sound are presented 100 times. In the below graph the proportion of 'simultaneous' responses (y) is plotted as a function of the asynchrony

I want to fit a 3-parameter Gaussian distribution to these data. I used least squares (below). My question is whether least squares is the most standard/orthodox form to proceed for this kind of fitting. It makes sense to use maximum likelihood estimation here? How can be done?

  x<-seq(-.3,.3,.1)
  yes<-c(10,45,90,100,60,10,5) #simultaneous responses
  y<-yes/100

  leastSquares<-function(p) sum((y-p[1]*dnorm(x,p[2],p[3]))^2)
  p<-optim(c(.001,0.1,.210),leastSquares)$par

  xseq<-seq(-.3,.3,.01)
  yseq<-p[1]*dnorm(xseq,p[2],p[3])

  plot(x,y)
  lines(xseq,yseq)

Answer 1

Use logistic regression with a Gaussian link.

---

The count of simultaneous responses for a given value of $x$, written $y(x)$, is the outcome of $n=100$ independent Bernoulli trials whose chance of success is given by the Gaussian function. Letting $\theta$ stand for the three parameters (unknown, to be estimated), let's write the value of that Gaussian at $x$ as $\mu(x, \theta)$. Then the probability of observing $y(x)$ successes, with $0 \le y(x) \le n$, is

$${\Pr}_\theta(y(x)) = \binom{n}{y(x)} \mu(x,\theta)^{y(x)} \left(1-\mu(x,\theta)\right)^{n-y(x)}.$$

The likelihood of a set of independent observations arising from the same underlying Gaussian is the product of these expressions. The part of the product that varies with $\theta$ is obtained from the last two factors. They will tend to be extremely small (because we are dropping the binomial coefficients), so about the only reasonable way to handle them on a computer is through their logarithms. Thus the part of the log likelihood that varies with $\theta$ is

$$\Lambda(\theta) = \sum_{x}\left(y(x)\log(\mu(x,\theta)) + (n-y(x))\log(1-\mu(x,\theta))\right).$$

Logistic regression with a Gaussian link maximizes this log likelihood.

To make sure the Gaussian peak does not exceed $1$, we might choose to parameterize its amplitude using a function that ranges from $0$ to $1$. Here is one convenient parameterization:

$$\mu(x, \theta) = \mu(x, (m,s,a)) = \frac{\exp(-\frac{1}{2} \left(\frac{x-m}{s}\right)^2)}{1 + \exp(-a)}.$$

The parameter $m$ is the mode of the Gaussian, $s$ (a positive number) is its spread, and $a$ (some real number) determines the amplitude, increasing with increasing $a$.

Use a multivariate nonlinear optimization procedure suitable for smooth functions. (Avoid specifying any constraints at all by using $s^2$ instead of $s$ as a parameter, if necessary.) Given some vaguely reasonable estimates of the parameters, it should have no trouble finding the global optimum.

---

As an example, here is a detailed implementation of the fitting procedure in R using data from the question. It is modified from code for a four-parameter least-squares fit of a Gaussian shown in an answer at Linear regression best polynomial (or better approach to use)?.

The fit is good: the standardized residuals do not become extreme and given the small amount of data, they are reasonably close to zero.

The estimated values are $\hat{m} = -0.033$, $\hat{s} = 0.127$, and $\hat{a} = 16.8$ (whence the estimated amplitude is $1/(1+\exp(-\hat{a})) = 1 - 0.00000005$).

y <- c(10,45,90,100,60,10,5) # Counts of successes at each `x`.
N <- length(y)
n <- rep(100, N)             # Numbers of trials at each `x`.
x <- seq(-3,3,1)/10
#
# Define a Gaussian function (of three parameters {m,s,a}).
#
gaussian <- function(x, theta)  { 
  m <- theta[1]; s <- theta[2]; a <- theta[3];
  exp(-0.5*((x-m)/s)^2) / (1 + exp(-a))
}
#
# Compute an unnormalized log likelihood (negated).
# `y` are the observed counts,
# `n` are the trials,
# `x` are the independent values,
# `theta` is the parameter.
#
likelihood.log <- function(theta, y, n, x) {
  p <- gaussian(x, theta)
  -sum(y*log(p) + (n-y)*log(1-p))
}
#
# Estimate some starting values.
#
m.0 <- x[which.max(y)]; s.0 <- (max(x)-min(x))/4; a.0 <- 0
theta <- c(m.0, s.0, a.0)
#
# Do the fit.
#
fit <- nlm(likelihood.log, theta, y, n, x)
#
# Plot the results.
#
par(mfrow=c(1,1))
plot(c(min(x),max(x)), c(0,1), main="Data", type="n", xlab="x", ylab="y")
curve(gaussian(x, fit$estimate), add=TRUE, col="Red", lwd=2) #$
points(x, y/n, pch=19)
#
# Compute residuals.
#
mu.hat <- gaussian(x, fit$estimate)
residuals <- (y - n*mu.hat) / sqrt(n * mu.hat * (1-mu.hat))
boxplot(residuals, horizontal=TRUE, xlab="Standardized residuals")
#
# (Compute the variance-covariance matrix in terms of the Hessian
# of the log likelihood at the estimates, etc., etc.)

(Although R has a way of performing these calculations automatically (by means of a custom link function for its glm Generalized Linear Model function), the interface is so cumbersome and so uninstructive I have elected not to use it for this illustration.)

Answer 2

One answer to your question is that least squares is a standard way to fit Gaussians to "some x and y data".

However, your data are special in that they are fractions (probabilities) that must lie in [0,1]. What you are doing is fitting a curve which is for probability densities, so it is unaware of the upper limit. The fitted maximum exceeds 1 and is for a value of x just less than 0. So, from one point of view you fitted an impossible curve. Densities greater than 1 are certainly possible but probabilities greater than 1 certainly are not.

You could scale your probabilities to densities, assuming that each is for a bin width 0.1. However, that would then ignore the fact that your largest value is the largest possible value, so you can't win everything either way.

What you should be doing is a largely open question, as you would need to tell us more about your problem and your data to get better advice. But various comments spring to mind:

1. Even if you have a really good reason to be focusing on Gaussians, I would think of other models too, asymmetric as well as symmetric.

2. The data are presented as if already on a scale with time interval 0.1 s. If your data are finer, you should use exact times. If this is the actual time resolution, then be advised that the data are rather coarse for discriminating between possible models. Naturally, there may be excellent scientific or technical reasons for your time interval.

3. One model popular in ecology for proportional abundances has been called the Gaussian logit, which in your terms is a logit for y fitted as a quadratic, with predictors x and x$^2$. Taking each value to represent 100 repetitions, I fitted this model in Stata:

clear 
set obs 7 
gen x = (_n - 4)/10 
mat y = (10,45,90,100,60,10,5) 
gen y = y[1,_n]
gen xsq = x^2
glm y x xsq , link(logit) family(binomial 100)
gen py = y/100 
twoway function invlogit(_b[_cons]+ _b[x]*x + _b[xsq]*x^2), ra(x) || scatter py x , scheme(s1color) legend(order(1 "predicted" 2 "observed")) xla(-.3(.1).3) ytitle(proportion) 

A reference for this model is Jongman, R.H.G., ter Braak, C.J.F. and van Tongerren, O.F.R. 1995. Data analysis in community and landscape ecology. Cambridge: Cambridge University Press.

The model fit is not especially encouraging. Note that in turn a limitation of, or more neutrally a fact about, this model is that it can never predict a maximum value of exactly 1 (or a minimum value of exactly 0). This model does concur with your least-square result in fitting a maximum for a negative x.

However, you (apparently) don't have enough data for much more flexible models to be tried comfortably, unless there are parsimonious physically- or psychometrically-based models tailored to the purpose that can be tried.