Adrian Mander <[email protected]> asks
> [...] the first line works (below) but the second line doesn't?
> [...]
> : (select( 0 , 1)) \ 1
> 1
> +-----+
> 1 | 0 |
> 2 | 1 |
> +-----+
>
> : (select( 0 , 0)) \ 1
> <istmt>: 3200 conformability error
> r(3200);
For those who don't know, Mata's select(A, v) function selects rows or columns
of a matrix, which depending on whether v is a row or column vector.
For instance, given the matrix
+- -+
| 1 2 3 |
A = | 4 5 6 |
| 7 8 0 |
+- -+
Then v = (1, 0, 1) produces
+- -+
| 1 3 |
| 4 6 |
| 7 0 |
+- -+
and v = (1, 0, 1)' (note the transpose operator) produces
+- -+
| 1 2 3 |
| 7 8 0 |
+- -+
With that introduction, it's easy enough to understand why the first line,
select(0,1)\1, works,
select(0,1) evaluates to 0
select(0,1)\1 evaluates to 0\1
The tricky part here is thinking about select(0,1), because the arguments
are both scalar. Nevertheless
+- -+
A = | 1 |
+- -+
which is to say, A is a 1x1 matrix, and
v = (1)
which is to say, v is a 1x1 row vector, and so we select the first row of A
and obtain
+- -+
| 1 |
+- -+
Alternatively, we could think of v as
v = (1)'
which is to say, v is a 1x1 column vector, and so we select the first column
of A and obtain
+- -+
| 1 |
+- -+
We get the same result however we think about it.
Next we are going to think about Adrian's second line, select(0,0)\1,
but before we do that, we need to take a detour
Void matrices
-------------
Matrices are an r x c array of elements, and usually we think of matrices
as being r x c subject to the contraint r>=1 and c>=1. Mata, however, allows
r, c, or both to be zero. Such matrices and vectors are called "void".
We can have 1 x 0 row vectors.
We can have 0 x 1 column vectors.
We can have 0 x 0, r x 0, and 0 x c matrices.
For example, let's consider a data matrix X and the regression formula
b = invsym(X'X)*X'y, where
X: n x v (n = # of obs., v = # of vars.)
y: n x 1
and thus
b: v x 1
An example of the usual case is X: 10x2, y: 10x1, and thus b is 2x1.
Let's think about extreme cases.
Case 1: X: 10x0 and y:10x1. In this case, we have 10 observations
and no explanatory variables. Let's work our way through invsym(X'X)*X'y.
X'X amounts to multiplying X' by X. X': 0x10, X: 10x0,
and so, X'X is 0x0
invsym(X'X), X'X: 0x0, results in a 0x0 matrix
X'y amounts to mulitplying X' by y. X': 0x10, y: 10x1.
Thus, X'y is 0x1.
invsym(X'X) * X'y amounts to multiplying invsym(X'X) by X'y.
invsym(X'X): 0x0, X'y: 0x1, and so the result b is 0x1.
In this case we obtain a void column vector for b.
Case 2: X:0x3 and y: 0x1. In this case, we have 0 observations on 3
explanatory variables, and the same number of observations on y.
X'X amounts to multiplying X' by X. X': 3x0, X: 0x3,
and so, X'X is 3x3. X'X contains zeros (yes zero, not missing).
invsym(X'X), X'X: 3x3 zero matrix, results in a 3x3 zero matrix
X'y amounts to mulitplying X' by y. X': 3x0, y: 0x1.
Thus, X'y is 3x1 and contains zeros.
invsym(X'X) * X'y amounts to multiplying invsym(X'X) by X'y.
invsym(X'X): 3x3 zero matrix, X'y: 3x1 zero vector, and so the result
b is a 3x1 zero vector.
In this case we obtain (0,0,0)' for b.
You can try both of these in Mata for yourself:
: X = J(10, 0, .)
: y = uniform(10,1)
: b = invsym(X'X)*X'y
: rows(b), cols(b)
1 2
+---------+
1 | 0 1 |
+---------+
and
: X = J(0, 3, .)
: y = J(0, 1, .)
: b = invsym(X'X)*X'y
: b
1
+-----+
1 | 0 |
2 | 0 |
3 | 0 |
+-----+
It is great fun to think about void matrices and whether, when creating
something from nothing, such as multiplying a 3x0 matrix by a 0x3 matrix to
obtain a 3x3 matrix, that something should be filled with zeros or missings.
That determination is made on the basis of limiting cases.
The purpose of void matrices is to make programming easier. With void matrices,
one usually does not ahve to worry about extreme cases. In some complicated
estimator, for instance, there may be an OPTIONAL second equation that amounts
to linear regression. In most programming languges, one would have to code
if (cols(X)) {
b = invsym(X'X)*X'y)
... continue to make calculation using b...
}
else {
... make limiting calculation when there is no b...
}
In Mata, however, one can usually code
b = invsym(X'X)*X'y)
... continue to make calculation using b...
and that will work even in the limiting case.
Back to Adrian's problem
------------------------
Andrians second line, select(0,0)\1, generated an error.
First let's understand why that happened:
select(0,0) says to select from matrix A = (0) no rows (or
equivalently, no columns). That produced a 0x0 matrix.
select(0,0)\1 says to join select(0,0): 0x0, with 1x1 column
vector 1. That a conformability error.
Adrian also wrote,
> because I need an if statement to handle this [...]
In telling you about void matrices, I told you that the advantages of
them is that one usually does not have to code around extreme cases.
This is an exception for reasons I will explain. But right now, let's
stay with Adrian's problem.
The generic form of what Adrian coded is
select(A, v) \ 1
and that coding assumes that v will select exactly one column from A.
That may not be so, the evidence being that Adrian tried an experiment that
selected 0 columns. Maybe Adrian was confused about what select() does,
but let's assume not.
Secondly, I assume Adrian simplified the problem for us by choosing A: 1x1 and
v: 1x1. Let's assume the general problem Adrian has in mind is A: r x c and
v: 1 x c. And, as I previously mentioned, we're going to assume that v is
such that it can select 0 or more columns, not merely 0 or 1. To handle all
of that, the general form of Adrian's statement would be
select(A, v) \ J(1, rowsum(v), 1)
When rowsum(v==0), I assume Adrian wants a (rows(A)+1) x 0 result. Adrian
needs to code
if (rowsum(v)) c = select(A, v) \ J(1, rowsum(v), 1)
else c = J(rows(A)+1, 0, 1)
It actually interesting to note that the statement
c = select(A, v) \ J(1, rowsum(v), 1)
will work in call case EXCEPT when v is 1x1 and contains 0.
For instance, if v=(0,0,0), the statement works. select(A, v) produces
a rows(A) x 0 result. J(1, rowsum(v), 1) produces a 1 x 0. The two can
be joined to produce a (rows(A)+1)x0 result.
The v: 1x1 and equal to 0 case is special because, in that case, -select()-
cannot tell whether v is a row vector or a column vector. It is both, and so
-select()- does not know whether to return a rows(A) x 0 or a 0 x cols(A)
result. select() returns a 0x0 result in this case. I suspect now we should
have included rowselect() and colselect() along with select(). Adrian could
write the one he needs for himself:
matrix colselect(matrix X, real rowvector v)
{
matrix result
result = select(X, v)
if (rows(result)==0 & col(result)==0) {
return(J(rows(X), 0, missingof(X))
}
return(result)
}
Then Adrian could simply code
c = select(A, v) \ J(1, rowsum(v), 1)
and he would get back what he (should have) expected.
-- Bill
[email protected]
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
