stata中mata语言学习-《Coding with Mata in Stata》

最新推荐文章于 2024-10-25 11:53:43 发布

肖夕木

最新推荐文章于 2024-10-25 11:53:43 发布

阅读量5.5k

点赞数 4

分类专栏： stata 计量文章标签：经验分享

本文链接：https://blog.csdn.net/weixin_45619112/article/details/105862481

版权

计量同时被 2 个专栏收录

4 篇文章

订阅专栏

stata

3 篇文章

订阅专栏

mata是stata的一种编程语言，类似于c或是c++等语言，作为一门编程语言，其同样有结构，指针，类。但本篇文章主要介绍其mata以及对矩阵的操作。
本篇文章只是mata语言的一个简明教程，可以帮助你认识程序中的各个大的框架，了解mata的一些基础知识。而没有深入的探讨mata中的各个函数。
关于mata的中文资料很少，只是看到连玉君老师分享过一篇文章。希望本文能对他人学习mata提供一些帮助。mata学习的英文资料，可关注我的公众号：“肖夕木的自习室” 回复mata获取

1. 为什么学习mata

开始学习mata是因为在stata中许多程序都是用mata编写的，若是想将这些程序改进以适用于自己的需要，则必然要学习mata

2. 在stata中使用mata

2.1 命令行的使用

mata可以直接在stata的命令窗口中使用，输入命令：

mata

即进入mata工作区，接下来输入的所有命令都将被解释成mata命令。
输入：

end

退出mata工作区，返回到stata平时使用的工作区。

注意：在mata工作区中定义的所有函数和矩阵在你退出mata工作区后，都会保存在mata工作区中，当再输入mata时，可以调用之前定义的所有函数和矩阵。

2.2 在.do文件中使用mata

在.do文件中可以多次使用mata命令，而且每次退出后，再进入mata工作区时，进入的是同一个工作区

2.3 在.ado文件中使用mata

有关使用mata生成新的命令的内容见section 7

3. mata中的矩阵

3.1 create矩阵，向量和标量

3.1.1 创建矩阵

有两种创建矩阵的方式，一种是创建整个矩阵，用\分行，用,分隔列。如创建一个2行2列的矩阵：

A=(1,2\3,4)

或者是创建一个空矩阵，例如创建一个2行三列的空矩阵：

B=J(2,3,.)

然后填充数据：

B[1,1] = 5
B[1,2] = 6
B[1,3] = 7
B[2,1] = 8
B[2,2] = 9
B[2,3] = 10

可以为矩阵命名，命名规则与大部分语言类似，即以字母，数字，下划线组合，但不能以数字为首字母。mata命名区分大小写。
如果矩阵已经被赋值，那么再对此矩阵赋值会覆盖之前的赋值。
输入变量的名字来显示

B

       1   2   3
    +-------------+
  1 |  5   6   7  |
  2 |  8   9   10 |
    +-------------+

3.1.2 创建向量

向量就是 $1\times n$ 或 $\times 1$ 的矩阵，向量的设置，此处不再赘述，但是要指出的是增加'代表转置，如：

f = (1, 2, 3)'

与

f = (  1 \ 2 \ 3  )

等价
也可以使用一系列数字来创建一个列向量或是行向量，如：

g = (3::6)
g
       1
    +-----+
  1 |  3  |
  2 |  4  |
  3 |  5  |
  4 |  6  |
    +-----+
g = (3..6)
g
       1   2   3   4
    +-----------------+
  1 |  3   4   5   6  |
    +-----------------+

3.1.3 创建标量

一个标量就是一个 $1\times1$ 的矩阵，如:
a = 2
但是需要注意的是矩阵命令区分标量和 $1\times1$ 矩阵

3.2 从stata访问数据

首先退出mata，然后获取数据集

webuse auto.dta

3.2.1 创建一个数据的副本

st data(colvector, rowvector)可以通过变量名复制在stata中的数据。如；

X = st_data(.,("mpg", "rep78", "weight"))

创建了一个包含变量mpg,rep78,weight的数据集。，也可以复制其子矩阵，如：

X = st_data((1::5\7::9),("mpg", "rep78", "weight"))

获取1至5与7至9行的子矩阵。
还可以使用函数st_data(matrix, rowvector)将对应编号的变量的值取出，如:

X = st_data((1::5\7::9),(3, 4, 7))
X
          1      2      3
    +----------------------+
  1 |    22      3   2930  |
  2 |    17      3   3350  |
  3 |    22      .   2640  |
  4 |    20      3   3250  |
  5 |    15      4   4080  |
  6 |    26      .   2230  |
  7 |    20      3   3280  |
  8 |    16      3   3880  |
    +----------------------+
 X = st_data((1,5\7,9),("mpg", "weight", "rep78"))
 X
           1      2      3
    +----------------------+
  1 |    22   2930      3  |
  2 |    17   3350      3  |
  3 |    22   2640      .  |
  4 |    20   3250      3  |
  5 |    15   4080      4  |
  6 |    26   2230      .  |
  7 |    20   3280      3  |
  8 |    16   3880      3  |
    +----------------------+

取出内存中已载入数据集auto1至5行与7至9行，与3，4，7列的数据存入矩阵X中。
更多请查看help mata:st_data

3.2.2 creating a view on the data

使用st view(X, colvector, rowvector)创建一个当前数据集的view(类似于把变量或某些行内的数据取出)，如：

X = .
st_view(X, ., ("mpg", "weight", "rep78"))

该命令的意思是创建一个包含变量mpg，weight，rep78的值的矩阵X。
使用st subview(Y, X, colvector, rowvector)根据矩阵X创建一个新矩阵Y，如：

st_view(X, ., .)
st_subview(Y, X, (1::5\7::9), (3,4,7))

用X矩阵的1到5行与7到9行的3，4，7列创建一个新矩阵Y。

3.2.3 copy vs. view

这里所说的copy与view是指st_data与st_view两个函数的区别。
st_data与st_view两个函数所实现的功能相似，但是，st_data运算更加快速，且不改变stata内存中的数据。

注意：对于字符串变量需使用st_sdata与st_sview函数。

3.3 管理mata的工作区

使用命令mata describe查看mata中的所有变量。
使用命令mata clear删除mata工作空间的所有变量。
特定矩阵的删除可以使用命令mata drop namelist

3.4 基本的矩阵操作

3.4.1 矩阵的转置

矩阵的转置前文已说

3.4.2 矩阵的分区

创建矩阵E

E = (1,2,5,6,7\3,4,8,9,10\3,4,1,2,3\5,6,4,5,6)
E
        1    2    3    4    5
    +--------------------------+
  1 |   1    2    5    6    7  |
  2 |   3    4    8    9   10  |
  3 |   3    4    1    2    3  |
  4 |   5    6    4    5    6  |
    +--------------------------+

使用[]进行矩阵的分区
如：

E[2, 3]
8
E[2,.]
1 2 3 4 5
+--------------------------+
1  |  3 4 8 9  1  0  |
+--------------------------+

注意可以通过range subscripts提取矩阵的2到3行与2到4列。如下述命令。

E[|2, 2\ 3, 4|]
1 2 3
+-------------+
1 | 4 8 9 |
2 | 4 1 2 |
+-------------+

也可以通过一个连续范围提取矩阵的2到3行与2到4列，如：

E[(2::3), (2..4)]
1 2 3
+-------------+
1 | 4 8 9 |
2 | 4 1 2 |
+-------------+

注意：第一种方法的运算速度快于第二种。

提取单个行和单个列组成矩阵：

E[(1\ 4), (3, 5, 2)]
1 2 3
+-------------+
1 | 5 7 2 |
2 | 4 6 6 |
+-------------+

更加一般化的是，可以用两个定义好的变量来提取行和列，如；

r = (1\ 4)
c = (3, 5, 2)
E[r,c]
1 2 3
+-------------+
1 | 5 7 2 |
2 | 4 6 6 |
+-------------+

注意；尽管stata认为将第一个下标向量指定为列向量，将第二个下标向量指定为行向量是一种好方法，但是这并不是必须的。

3.4.3 提取对角线并创建对角矩阵

diagonal(A)可以提取矩阵的对角线上的元素，并以列向量的形式表示，如：

A = (1,2,3\4,5,6\7,8,9)
A
       1   2   3
    +-------------+
  1 |  1   2   3  |
  2 |  4   5   6  |
  3 |  7   8   9  |
    +-------------+
b = diagonal(A)
b
       1
    +-----+
  1 |  1  |
  2 |  5  |
  3 |  9  |
    +-----+

diag(b)创建一个对角线为b向量的矩阵，如:

 diag(b)
[symmetric]
       1   2   3
    +-------------+
  1 |  1          |
  2 |  0   5      |
  3 |  0   0   9  |
    +-------------+

3.4.4 提取上下三角矩阵

lowertriangle(A)提取下三角矩阵，uppertriangle(A)提取上三角矩阵。

3.4.5 对矩阵进行排序

sort（X,idx）返回一个以idx列排序的矩阵X，如sort（X,1）即以第一列对X进行排序

X = (2, 3, 1\ 2, 2, 2\ 1, 1, 3)
X
1 2 3
+-------------+
1 | 2 3 1 |
2 | 2 2 2 |
3 | 1 1 3 |
+-------------+
sort(X,1)
1 2 3
+-------------+
1 | 1 1 3 |
2 | 2 3 1 |
3 | 2 2 2 |
+-------------+

使用help mata sort()可获取更多信息。

3.5 基本的矩阵运算符号

+， $\times$ ，-，/，^

接下来的运算符号，以例子的形式讲解。
设定矩阵：

a = (1..5)
b = (6::10)
c = 3

矩阵的 + 与 $\times$ 与数学上的定义相同

d=a+c*a
d
1 2 3 4 5
+--------------------------+
1 |  4  8  12  16  20 |
+--------------------------+

3.6 矩阵乘法

函数cross(X, Z)可以计算X’Z。cross(X, Z)函数计算有以下三点优势

自动忽略缺省值
运算速度快
占用空间小

3.6.1 Element-by-element operators

即矩阵间，各个内部元素（后文称之为colon）的相互运算，而不是两个矩阵的运算。运算符： :* :/ :^，例如：

x = (1, 2, 3)
y = (4, 5, 6)
x:*y
1 2 3
        1    2    3
    +----------------+
  1 |   4   10   18  |
    +----------------+

就是各个元素对应的乘积。

3.6.2 关系与逻辑运算符

等价运算符：
等于： ==
不等于 !=
此关系运算符用于标量，向量和矩阵。
关系运算符：
< > >= <=
colo relational operator，在之前的运算符前加上:即可，如:== 例子：

G = (1, 2 \ 3, 4)
H = (1, 5 \ 6, 4)
G :== H
       1   2
    +---------+
  1 |  1      |
  2 |  0   1  |
    +---------+

逻辑运算符：
& && 和
| || 或
同样可以加:

3.7 一些特殊矩阵的创建

3.7.1 创建一个含相同元素的矩阵

 J(2, 3, 0)
       1   2   3
    +-------------+
  1 |  0   0   0  |
  2 |  0   0   0  |
    +-------------+

3.7.2 创建一个单位阵

 I(3)
[symmetric]
       1   2   3
    +-------------+
  1 |  1          |
  2 |  0   1      |
  3 |  0   0   1  |
    +-------------+

3.7.3 创建单位向量

 e(2, 5)
       1   2   3   4   5
    +---------------------+
  1 |  0   1   0   0   0  |
    +---------------------+

3.7.4 创建随机矩阵

 uniform(2, 3)
                 1             2             3
    +-------------------------------------------+
  1 |  .3488717046   .2668857098   .1366462945  |
  2 |  .0285568673   .8689332692    .350854896  |
    +-------------------------------------------+

uniform()创建的是一个0，1之间的随机数。
uniformseed(newseed)与stata中一样可以设立一个随机种子，以使结果具有可重复性。
mata无法直接生成随机正态分布的矩阵，但可以通过嵌套的方法来实现，如：

 invnormal(uniform(2,3))
                  1              2              3
    +----------------------------------------------+
  1 |  -1.467610024   -.4583014284    .1385652864  |
  2 |   1.155176934   -.8249166005    1.241333377  |
    +----------------------------------------------+

3.8 逆与线性代数

mata有多种计算矩阵逆的函数：
luinv(A) A为满秩，方阵
cholinv(A) A为正定对称矩阵
invsym(A) generalized inverse of positive-definite, symmetric matrix A
mata也提供了一些方法来求解线性方程组 $A X = B$ ：
lusolve(A,B) A为满秩，方阵。
cholinv(A) A为正定对称矩阵。
help m4 solvers获得更多相关函数
一个使用auto.data数据进行线性回归的例子：

y = st_data(.,"price")
X = st_data(.,("mpg", "weight"))
X = X, J(rows(X),1,1)
b = invsym(X’*X)*X’*y

4. Controlling the flow

此处的流程控制语句与c中类似

4.1 循环语句

while-loop：

while (expr) {
stmt
}

例如:

n = 5
i = 1
while (i<=n) {
printf("i=%g\n", i)
i++
}
printf("done\n")

for-loop

for (expr1; expr2; expr3) {
stmts
}

例如:

n = 5
for (i=1; i<=n; i++) {
printf("i=%g\n", i)
}
printf("done\n")

do循环

do {
stmt
} while (exp)

4.2 条件语句

if condition

 if (expr) {
stmt
}

5. 编码 mata 函数

mata允许创建一些新函数，例如:

function zeros(c)
{
a = J(c, 1, 0)
return(a)
}

输出结果

b = zeros(3)
 b
       1
    +-----+
  1 |  0  |
  2 |  0  |
  3 |  0  |
    +-----+

5.1 declarations

mata中不需要像其他语言中那样声明创建变量的类型，但是mata建议你声明。
element type：transmorphic, numeric, real,complex, string, pointer
organizational type：matrix, vector, rowvector, colvector, scalar
但是，mata中可以强制要求声明：

set matastrict on

6. Controlling the flow

此处的流程控制语句与c中类似

6.1 循环语句

while-loop：

while (expr) {
stmt
}

例如:

n = 5
i = 1
while (i<=n) {
printf("i=%g\n", i)
i++
}
printf("done\n")

for-loop

for (expr1; expr2; expr3) {
stmts
}

例如：

n = 5
for (i=1; i<=n; i++) {
printf("i=%g\n", i)
}
printf("done\n")

do-循环

do {
stmt
} while (exp)

6.2 条件语句

 if condition  
 if (expr) {
stmt
}

7. 编码 mata 函数

mata允许创建一些新函数，例如：

function zeros(c)
{
a = J(c, 1, 0)
return(a)
}
输出结果

b = zeros(3)
 b
       1
    +-----+
  1 |  0  |
  2 |  0  |
  3 |  0  |
    +-----+

7.1 declarations

real colvector zeros(real scalar c)
{
real colvector a
a = J(c,1,0)
return(a)
}

但是，mata中可以强制要求声明：

set matastrict on

7.2 传递参数（stata中的argument）

function中如果有多个需要传递的参数，则这些参数可以用,隔开，对于之前的zero()如：

real matrix zeros(real scalar c, real scalar r)
{
  real matrix A
  A = J(c, r, 0)
  return(A)
}

mata中的函数还允许设定可选的参数，可选参数通过’|'隔离来设置，同样对于之前的zero()函数：

real matrix zeros(real scalar c,| real scalar r)
{
real matrix A
if (args()==1) r = 1
A = J(c, r, 0)
return(A)
}

函数args()决定参数的数量。help m2 syntax获取更多信息。

注意： mata传递的是参数的地址，而不是参数的值。

7.3 返回值

使用函数return(expr)会结束当前函数进程，且返回相应的值。如：

real matrix zeros(real scalar c,| real scalar r)
{
if (args()==1) {
    return(J(c, 1, 0))
}
else {
    return(J(c, r, 0))
}
}

此时会返回一个矩阵
如果一个函数不返回任何东西，则称此函数为void，例如：

void zeros(real matrix A, real scalar c,| real scalar r)
{
if (args()==1) {
A = J(c, 1, 0)
}
else {
A = J(c, r, 0)
}
}

此时函数不会返回任何值。

7.4 .mo和.mlib文件

使用mata mosave可以将函数保存至文件中，这样在stata中可以随时调用，而不用重新定义，例如创建一个.do文件：

version 10
mata:
real matrix zeros(real scalar c,| real scalar r)
{
real matrix A
if (args()==1) r = 1
A = J(c, r, 0)
return(A)
}
mata mosave zeros(), replace
end

运行后，会创建一个zeros.mo文件存储在工作目录下。如果你关闭stata，再打开stata，该函数同样可以运行。
多个函数文件可以存储在.mlib文件中。

8. 使用mata编码stata命令

8.1 一个例子

创建一个ols回归的命令,文件保存至一个.ado文件，即可调用。

program define ols, eclass
version 10.0
syntax varlist(numeric) [if] [in]
gettoken depvar indepvar: varlist
marksample touse
mata: m_ols("‘varlist’", "‘touse’")
tempname b V
matrix ‘b’ = r(b)’
matrix ‘V’ = r(V)
local N = r(N)
matname ‘b’ ‘indepvar’ _cons, c(.)
matname ‘V’ ‘indepvar’ _cons
ereturn post ‘b’ ‘V’, depname(‘depvar’) obs(‘N’) esample(‘touse’)
ereturn local cmd = "ols"
ereturn display
end



capture mata mata drop m_ols()
version 10
mata:
void m_ols(string scalar varlist, string scalar touse)
{
real matrix  M, X, V
real colvector  y, b
real scalar  n, k, s2
M = X = y = .
st_view(M, ., tokens(varlist), touse)
st_subview(y, M, ., 1)
st_subview(X, M, ., (2\.))
n = rows(X)
k = cols(X)
XX = cross(X,1,X,1)
if (rank(XX) < k+1) {
errprintf("near singular matrix\n")
exit(499)
}
Xy = cross(X,1,y,0)
b = cholsolve(XX,Xy)
e = y - (X, J(n,1,1))*b
s2 = (e’e)/(n-k)
V = s2*cholinv(XX)
st_eclear()
st_matrix("r(b)", b)
st_matrix("r(V)", V)
st_numscalar("r(N)", n)
}
end