Deriving least squares coefficients for curve of form $y=a/(x^2+b)$

118 Views Asked by At

I am trying to derive the coefficients $a$ and $b$ for a least squares curve of the form $y=\frac{a}{x^2+b}$

To start I have written $$S=\sum_i(y_i-y)^2=\sum_i\left(y_i-\frac{a}{x_i^2+b}\right)^2=\sum_i\left(y_i^2-\frac{2a y_i}{x_i^2+b}+\frac{a^2}{(x_i^2+b)^2}\right) \\S=\sum_iy_i^2-2a\sum_i\frac{y_i}{x_i^2+b}+a^2\sum_i\frac{1}{(x_i^2+b)^2}$$

Now the way I have seen it done for the line $y=a+bx$ is they minimize $S$ with respect to both $a$ and $b$, applying that for this curve gives the two constraint equations:

$$0=a\sum_i\frac{1}{(x_i^2+b)^2}-\sum_i\frac{y_i}{x_i^2+b}\\ 0=-a\sum_i\frac{1}{(x_i^2+b)^3}+\sum_i\frac{y_i}{x_i^2+b}$$

Summing these together we find $$0=a\left(\sum_i\frac{1}{(x_i^2+b)^2}-\sum_i\frac{1}{(x_i^2+b)^3}\right)$$ Since $a\neq 0$ we must have $$\sum_i\frac{1}{(x_i^2+b)^2}=\sum_i\frac{1}{(x_i^2+b)^3}$$ Which seems to imply $$(x_i^2+b)^2=(x_i^2+b)^3$$ and unfortunately there doesnt appear to be a constant $b$ that solves this equation. Did I make a mistake somewhere or is this the result?

3

There are 3 best solutions below

0
On BEST ANSWER

Starting with a rather silly sentence : since the model is nonlinear with respect to $b$, by the end, you will need a nonlinear regression and, as usual, you will need reasonable starting guesses for the parameters.

You can get these estimates in a preliminary step $$y=\frac{a}{x^2+b} \implies \frac 1y=\frac {x^2} a+\frac ba$$ So, defining for each $(x_i,y_i)$ data points, new variables $t_i=x_i^2$ and $z_i=\frac 1{y_i}$, you have $$z=\alpha\, t+\beta$$ and, from the corresponding linear regression, you will get, as estimates $a=\frac 1 \alpha$ and $b=\frac \beta \alpha$.

If you do not want to use nonlinear regression, you can reduce the problem to a nasty equation in $b$.

Using Rohan's expressions, let us define $$S_1=\sum \frac{y_i}{x_i^2+b} \qquad S_2=\sum \frac{1}{(x_i^2+b)^2}$$ $$S_3=\sum \frac{y_i}{(x_i^2+b)^2} \qquad S_4= \sum \frac{1}{(x_i^2+b)^3}$$ Eliminating $a$, you then need to to find the zero of function $$f(b)=S_2\, S_3-S_1 S_4$$ Since you have the estimate, you could use it as the starting point of Newton method with analytical derivatives since $$f'(b)=S_2S'_3+S'_2S_3-S_1S'_4-S'_1S_4$$ where will appear other summations (I let you the task of writing them - they are not difficult). The iterative process will converge very fast. So $b$ and then $a$.

Edit

For illustration purposes, consider the following data set $$\left( \begin{array}{cc} x & y \\ 0 & 46 \\ 1 & 46 \\ 2 & 45 \\ 3 & 43 \\ 4 & 41 \\ 5 & 38 \\ 6 & 36 \\ 7 & 33 \\ 8 & 30 \\ 9 & 28 \end{array} \right)$$ The preliminary step leads to $z=0.000176905 t+0.0216373$ from which the estimates $a=5652.76$ and $b=122.310$.

Starting with this guess, the iterates of Newton method will be $$\left( \begin{array}{cc} n & b_n \\ 0 & 122.31000 \\ 1 & 121.84307 \\ 2 & 121.85174 \end{array} \right)$$ from which $a=5635.22$. A nonlinear regression would have be providing exactly the same values for the parameters.

4
On

We have $$\frac{\partial S}{\partial a} = -2\sum \frac{y_i}{x_i^2+b}+2a\sum \frac{1}{(x_i^2+b)^2}$$ $$\frac{\partial S}{\partial b} = 2a \sum \frac{y_i}{(x_i^2+b)^2} -2a^2 \sum \frac{1}{(x_i^2+b)^3}$$

3
On

If you actually had to run this regression, you would use this model: $$ ln(y_i) = ln(a) - ln(x_i^2+b)+u_i $$ where $u_i$ is the error term. The sum of squared residuals is $$ S =\sum_i(ln(y_i) - ln(\hat{a}) + ln(x_i^2+\hat{b}))^2 $$ where I am using the usual notation for an estimator. Then $$ \partial S/\partial \hat{a} = 0 \\ \sum_i(ln(y_i) - ln(\hat{a}) + ln(x_i^2+\hat{b})) =0\\ ln(\hat{a}) = \sum_i(ln(y_i)+ln(x_i^2+\hat{b}))/n $$ where $n$ is the number of observations. Also, $$ \partial S/\partial \hat{b} = 0 \\ \sum_i\dfrac{ln(y_i) - ln(\hat{a}) + ln(x_i^2+\hat{b}))}{(x_i^2+\hat{b})} =0\\ $$ The latter expression is not linear in $\hat{b}$. So I do not see an easy algebraic solution for $\hat{a}$ and $\hat{b}$. One would not expect the solution to be linear in the parameters because the original functional form is itself non-linear. Maybe someone else can pick it up from here.