1
$\begingroup$

Apologies in advance for the tedious beginner question.

I'm trying to translate a least-squares problem from a manual process (using Excel for matrix transposition and multiplication) to using the Python statsmodels package. In this case, I'm performing an affine transform from a set of observed coordinates to a set of ground coordinates in eastings (E) and northings (N). I've used the following formula to form the A (design) matrix: $$ \begin{equation} f_i(a_{0}, a_{1}, a_{2}) = a_{0} + a_{1}x_{i} + a_{2}y_{i} \\ f_i(b_{0}, b_{1}, b_{2}) = b_{0} + b_{1}x_{i} + b_{2}y_{i} \end{equation} $$ Which gives me a matrix that looks like:

$$ \begin{bmatrix} 1 & E_1 & N_1 & 0 & 0 & 0 \\ 1 & E_n & N_n & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & E_1 & N_1 \\ 0 & 0 & 0 & 1 & E_n & N_n \\ \end{bmatrix} $$

and a b vector, which is all x co-ordinates, followed by all y co-ordinates:

$$ \begin{bmatrix} x_1 \\ x_n \\ y_1 \\ y_n \end{bmatrix} $$

I also construct a square covariance matrix, W, using the square of the standard errors of the eastings and northings, arranged in a diagonal. This is used as the weight matrix during the least-squares process (the standard errors are assumed to be independent)

I then calculate:

$A^TWA$, $A^TWb$, and $(A^TWA)^{-1}$, then multiply $(A^TWA)^{-1}$ by $A^TWb$ to determine a vector x, which contains values for $a_0, a_1, a_2$ and $b_0, b_1, b_2$ 5. Multiply A by x, and subtract b from the result, to determine a residuals vector, v.

I can calculate the unit variance ($\sigma0$) by obtaining the square root of $\frac{v^TWv}{observations - unknowns}$, and multiplying it by the a priori standard error of each co-ordinate, in order to assess the quality (a posteriori standard error) of the transform. I can also calculate the standard error of my x vector (the diagonal values of Cx) by multiplying $(A^TWA)^{-1}$ by $\sigma0^2$

Now that the tedious step-by-step manual explanation is out of the way, let's say I have Pandas DataFrames for A, b and W:

In [124]: A_matrix
Out[124]: 
<class 'pandas.core.frame.DataFrame'>
Int64Index: 108 entries, 0 to 107
Data columns (total 6 columns):
As             108  non-null values
Eastings_a     108  non-null values
Northings_a    108  non-null values
Bs             108  non-null values
Eastings_b     108  non-null values
Northings_b    108  non-null values
dtypes: float64(4), int64(2)

In [125]: b_vector
Out[125]:
<class 'pandas.core.frame.DataFrame'>
Int64Index: 108 entries, 0 to 107
Data columns (total 1 columns):
coordinates    108  non-null values
dtypes: float64(1)

In [162]: Weight_matrix
Out[162]: 
<class 'pandas.core.frame.DataFrame'>
Int64Index: 108 entries, 0 to 107
Columns: 108 entries, 0 to 107
dtypes: float64(108)

How do I use statsmodels.ols_regression to easily calculate my residuals and $\sigma 0$?

$\endgroup$
  • 1
    $\begingroup$ I'm wondering why you don't use numpy (which is simpler to understand if you are not coming from an R background) docs.scipy.org/doc/numpy/reference/generated/… $\endgroup$ – seanv507 Aug 4 '13 at 23:12
  • $\begingroup$ @seanv507 That works, but I have to work out the residuals (Ax-b) manually. $\endgroup$ – urschrei Aug 5 '13 at 0:12
  • 1
    $\begingroup$ Statsmodels performs all the matrix manipulations for you. Once you fit the model, the results instance has all the information you are looking for. For example, you can get the residuals by calling .resid(). $\endgroup$ – TomAugspurger Aug 5 '13 at 2:33
  • $\begingroup$ @urschrei... numpy returns residuals too? as an extra output parameter (unless i have misunderstood your whole setup) $\endgroup$ – seanv507 Aug 5 '13 at 13:07
  • $\begingroup$ @seanv I've put my progress so far on Github: github.com/urschrei/linalg. When I ran np.linalg.lstsq, using my design matrix and outcome vector, I got back parameters, but not residuals. It's quite likely I did something wrong, though. $\endgroup$ – urschrei Aug 5 '13 at 14:01
0
$\begingroup$

Calculating σ0 is simple: take the square root of the returned fit object's mse_resid property.

Calculating the residuals: It's just the resid property of the returned fit object.

| cite | improve this answer | |
$\endgroup$

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.