7
$\begingroup$

I am conducting a non-linear least squares regression fit using the python scipy.optimize.curve_fit function, and am trying to better understand the weights that go into this method.

I have a distribution of raw data points that I wish to fit to a Gaussian cumulative distribution function. I created a function for this that takes three parameters: an average, a standard deviation, and a scale factor for if my distribution doesn't quite approach 1.

My confidence in each raw data point is based on a separate instrumental count that doesn't necessarily have anything to do with the value of the data point, so I'm trying to include this in my fit using weights. Specifically, small and large values of x have less certainty, so I want them to matter less in the regression. When I conduct a fit by passing these raw counts as weights, the fit is not particularly good, whereas if I pass them as (1/counts) the fit improves. I have plotted the raw data, fits, and normalized weights for these two options.

What I am trying to understand is how to interpret the weights. I would have thought that higher values for the weights imply more importance in the regression. However, it seems that it is actually correct to use the weights where the "bad" raw data points have a higher weight. Why is this and how should I be interpreting the weights? Also, is there a better resource for understanding weights in non-linear regression? Most of the sources I have found have not explained weighting in a way I can comprehend.

Edit: I added a second plot that now shows actual (non-normalized) counts, along with the corrected fit (weighted according to counts) according the proper fitting technique shown below.

enter image description here

enter image description here

$\endgroup$
  • $\begingroup$ (1) Something is not quite right about the graphic. It has two vertical axes, but the right hand one cannot possibly be for counts or reciprocal counts (it would show them as negative, which they cannot be). We must conclude it shows the values of the fits. However, with only the three parameters mentioned in the text the fits necessarily approach an asymptote of $0$ at the left, whereas these both approach $-2$. What's going on? (2) What do the counts mean? Don't you have some prior understanding of how they relate to the precision of the measurements? (3) What is "x"? $\endgroup$ – whuber Jul 25 '13 at 0:47
  • 1
    $\begingroup$ Sorry, I should have explained better. The right axis shows normalized variables so that counts and inverse counts can fit onto the same axis. The actual counts can vary from 0 to a some larger value (here between about 2 and 85) and are what I am actually trying to use as my weights at each data point. The "x" axis is a size where independent measurements are made for each raw data point. These counts relate to the precision and accuracy of my measurement, so as counts get smaller I expect both more noise and less accuracy in the raw data values shown on the left axis. $\endgroup$ – Thursdays Coming Jul 25 '13 at 1:05
  • 1
    $\begingroup$ Thanks. Are you sure you haven't mixed up the two fits? The correct fit will be weighted by the counts, which themselves strongly downweight the measurements at the right side (for large $x$), indicating the count-weighted fit must be the blue line, not the green line. Either way, the fit is bad and contradicts some of your assumptions about the data, so you have deeper issues to deal with than choosing appropriate weights. $\endgroup$ – whuber Jul 25 '13 at 14:05
  • $\begingroup$ I think you are correct that I must have done something wrong with the fits. I confirmed that the green line in the first plot was indeed using the actual weights and they weren't mixed up. But I think that my problem must be the function of the distribution I was trying to fit against? When I did a standard scipy gaussian CDF, the fit weighted by the standard counts rather than inverse counts seems to fit better. I would still like to try to fit the data to a function that does not approach 1, but I think I must be going about it the wrong way. $\endgroup$ – Thursdays Coming Jul 25 '13 at 17:29
  • $\begingroup$ Clearly in your new bottom image you are not allowing a vertical rescaling: only two parameters (location and spread) have been varied. $\endgroup$ – whuber Jul 25 '13 at 17:46
6
$\begingroup$

The weights should equal the counts, because those will be inversely proportional to the variances of the errors. Specifically, the model for the data $(x_i, y_i, n_i)$ is

$$y_i \sim \lambda \Phi((\log(x_i) - \mu)/\sigma + \varepsilon_i$$

with $\mu, \sigma \gt 0,$ and $\lambda \gt 0$ the parameters and $\varepsilon_i$ are independent random variables with zero means and variances

$$\text{Var}(\varepsilon(i)) = \sigma^2 / n_i$$

where $n_i$ are the counts.

The fit to the logarithm of $x$ is visually ok:

Figure

In this figure the x-axis is on a logarithmic scale, the point symbols have areas proportional to the counts (so that large circles will have more influence in the fitting than small ones), and the red line is a least-squares fit. It is clear the model is not really appropriate: the residuals for smaller values of $y$ tend to be small, regardless of the counts. Possibly the sum of squares of relative errors should be minimized to obtain an appropriate fit.

It is evident that the fit is poor for the largest $x$, but those also have small counts.


The R code with (my version of) the data and the fitting and plotting procedures follows.

y <- c(1, 1, 2, 1, 2, 1, 3, 4, 22, 30, 44, 58, 68, 69, 
       71, 72, 75, 72, 80, 78, 87, 86, 80, 82, 92, 90, 85, 61, 38, 36) / 100
x <- ceiling(exp(seq(log(20), log(500), length.out=length(y))))
counts <- c( 10, 3, 17, 20, 38, 31, 44, 55, 58, 68, 77, 
             82, 86, 82, 77, 75, 70, 65, 68, 51, 47, 41, 38, 30, 22, 14, 9, 4, 2, 1)
#
# The least-squares criterion.
# theta[1] is a location, theta[2] an x-scale, and theta[3] a y-scale.
#
f <- function(theta, x=x, y=y, n=counts) 
  sum(n * (y - pnorm(x, theta[1], theta[2]) * theta[3])^2) / sum(n)
#
# Perform a count-weighted least-squares fit.
#
xi = log(x)
fit <- optim(c(median(xi), sd(xi), max(y) * sd(xi)), f, x=xi, y=y, n=counts)
#
# Plot the result.
#
par(mfrow=c(1,1))
plot(x, y, log="x", xlog=TRUE, pch=19, col="Gray", cex=sqrt(counts/12))
points(x, y, cex=sqrt(counts/10))
curve(fit$par[3] * pnorm(log(x), fit$par[1], fit$par[2]), 
          from=10, to=1000, col="Red", add=TRUE)
$\endgroup$
  • 1
    $\begingroup$ Ah, I'm starting to understand the proper way to fit this now. Using the relative errors that are scaled by the counts for the minimization makes much more sense. Previously, I was just passing raw counts rather than variance at each point in as a weight, so I assume this is where my problem arose. Thanks for this explanation, I think I understand it much better now. $\endgroup$ – Thursdays Coming Jul 25 '13 at 23:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.