Statalist The Stata Listserver


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: Two Mata questions


From   wgould@stata.com (William Gould, Stata)
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Two Mata questions
Date   Thu, 30 Mar 2006 09:35:07 -0600

AbdelRahmen El Lahga <uaquap@gnet.tn> asks, 

> I have 2 questions about matrix manipulation on my unbalanced panel data
> sets:  with my old Stata matrix language I can create a sequence of matrix
> by typing for exemple
>
>       local wave=1
>       while `wave'<22 {
>               matrix accum A`wave' = x1 x2 ..xk if t==`wave', noconstant
>               local wave=`wave' +1
>       }
>
>  How can i perform a similar task with Mata?  I don't seek for a complex way
>  to things but in fact i would like to write a more complex expression wich
>  looks like:
>
>       for each time period t the matrix
>            A_{t} = sum_{i=1 to n} (s1_it * s2_it  x_it' * x_it)
>
>   where , s1_it is a scalar =norm(xb) and s2_it= normden(xb) scalars s1_it
>   and s2_it comes froma previeous probit estimation any help?

Okay, hold on to your hat.  This is not difficult, but we are going to use 
some features with compter-scienced technical terms and you might confuse 
the jargon with inherent difficulty.

First off, in Mata, you cannot code things like A`wave'.  A`wave' is a macro
trick, and the value of the macro is substituted at execution time.  Mata is
compiled before execution -- that's what makes Mata fast -- but the side
effect is that compilation prevents trickery like A`wave'.

What AbdelRahmen (am I using the name correctly?) wants is an array of
matrices.  In the next free executable update of Stata scheduled for mid
April, there is a new Mata feature called structures.  They are exactly what
AbdelRahmen needs.  If he can wait, the solution is


        struct amatrix {
                real matrix X
        }

        function ...() 
        {
                struct amatrix vector   A

                ...
                A = amatrix(20)             // <- make 20 of them
                ...

                ... A[1].X ... A[2].X ...   // <- how to refer to them

                ...

                for (i=1; i<=20; i++) { 
                      SUM = SUM + A[i].X       // <- example of use
                }
                ...
         }

So that's the midapril answer.  Structures can do a lot more than solve 
this problem, and I'll be talking about them later.

The answer for today is POINTERS.

A pointer is a memory address value, and said memory address is the 
address of some other Mata object, such as a real matrix.

A pointer might be the value 501,128.  The implication is that, at memory 
address 501,128 is a Mata matrix.

Programmers don't usually write memory addresses in decimal,
although they could.  Rather than write 501,128, they write 0x7a588, which is
the same number in base 16.

Now, when you use pointers, you don't make up address willy nilly.  
Instead, you get addresses and you store them.  If I coded

          function ...(...)
          {
                  real matrix     A
                  pointer vector  p

                  p = J(1, 20, NULL)
                  for (i=1; i<=20; i++) {
                          p[i] = &A
                  }
                  ...
          }

I would create a 20-element vector, each element of which is the address 
of matrix A.  The first assignment line, 

                 p = J(1, 20, NULL)

created a 1x20 vector, each element of which is the special pointer value 
0x0, zero, or in progammer jargon, NULL.  0x0 (NULL) points to nothing.
Then, in 

                  for (i=1; i<=20; i++) {
                          p[i] = &A
                  }

I swept over the vector and filled in each element with the adress of A.
& in this case does not mean "and", it means means "address of".  
& always means "address of" when you use it as a prefix operator instead of 
between two expressions, when it means "and".

Going along with prefix operator & is prefix operator *.  *p[1] means 
"the contents of p[1]" or "what p[1] points to".  What p[1] points to 
is A.  

Typing *p[1] or typing A refer to the same entity.  Change the [1,1] element 
of *p[1], 

                  (*p[1])[1,1] = 5

and you will discover that A[1,1] is now 5.  And vice versa.

Think of the * prefix operator as being analogous to single-quote macro
substitution.  Think of the & prefix operawtor as being analogous to 
macro definition.

Now, in this case, I made each and every element of p[i] point to the same
matrix.  That can be useful in some cases, but that is not what AbdelRahmen El
needs.

Abdel wants to make 20 different matrices and then do some sort of for-loop
calculation across them.  Here's how:

        1.  Create a subroutine that returns the value of one of the 
            created matrices.  Call it repeatedly to obtain the 
            20 matrices.

        2.  Store the addresses of the matrices returned by the subroutine.

        3.  Then use the matrcies returned.

For example, says that mysub(i) calculates and returns the i-th matrix:

         real matrix mysub(real scalar i)
         { 
                 ...
                 return(A)
         }

Then in the calling routine, we do the following, 

         function ...(...)
         {
                 pointer vector   p

                 p = J(1, 20, NULL)

                 for (i=1; i<=20; i++) {
                         p[i] = &mysub(i)
                 }

                 ...

                 for (i=1; i<=20; i++) {
                         SUM = SUM + *p[i]
                 }

                 ...
           }

Note the line 

                         p[i] = &mysub(i)

Rather than saving the matrix returned by mysub(), I saved its address, 
and I stacked the different addresses, one after another, in the pointer 
vector p[].

-- Bill
wgould@stata.com
*
*   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/



© Copyright 1996–2020 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index