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] = 5B[1,2] = 6B[1,3] = 7B[2,1] = 8B[2,2] = 9B[2,3] = 10
可以为矩阵命名,命名规则与大部分语言类似,即以字母,数字,下划线组合,但不能以数字为首字母。mata命名区分大小写。
如果矩阵已经被赋值,那么再对此矩阵赋值会覆盖之前的赋值。
输入变量的名字来显示
B1 2 3+-------------+1 | 5 6 7 |2 | 8 9 10 |+-------------+
3.1.2 创建向量
向量就是 1 × n 1\times n 1×n或 n × 1 n \times 1 n×1的矩阵,向量的设置,此处不再赘述,但是要指出的是增加'
代表转置,如:
f = (1, 2, 3)'
与
f = ( 1 \ 2 \ 3 )
等价
也可以使用一系列数字来创建一个列向量或是行向量,如:
g = (3::6)g1+-----+1 | 3 |2 | 4 |3 | 5 |4 | 6 |+-----+g = (3..6)g1 2 3 4+-----------------+1 | 3 4 5 6 |+-----------------+
3.1.3 创建标量
一个标量就是一个 1 × 1 1\times1 1×1的矩阵,如:
a = 2
但是需要注意的是矩阵命令区分标量和 1 × 1 1\times1 1×1矩阵
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))X123+----------------------+1 | 223 2930 |2 | 173 3350 |3 | 22. 2640 |4 | 203 3250 |5 | 154 4080 |6 | 26. 2230 |7 | 203 3280 |8 | 163 3880 |+----------------------+X = st_data((1,5\7,9),("mpg", "weight", "rep78"))X123+----------------------+1 | 22 29303 |2 | 17 33503 |3 | 22 2640. |4 | 20 32503 |5 | 15 40804 |6 | 26 2230. |7 | 20 32803 |8 | 16 38803 |+----------------------+
取出内存中已载入数据集auto
1至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)E1 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]8E[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)A1 2 3+-------------+1 | 1 2 3 |2 | 4 5 6 |3 | 7 8 9 |+-------------+b = diagonal(A)b1+-----+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)X1 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*ad1 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:*y1 2 31 2 3+----------------+1 | 4 10 18 |+----------------+
就是各个元素对应的乘积。
3.6.2 关系与逻辑运算符
等价运算符:
等于:==
不等于!=
此关系运算符用于标量,向量和矩阵。
关系运算符:
<
>
>=
<=
colo relational operator,在之前的运算符前加上:
即可,如:==
例子:
G = (1, 2 \ 3, 4)H = (1, 5 \ 6, 4)G :== H1 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 AX=B AX=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 = 5i = 1while (i<=n) {printf("i=%g\n", i)i++}printf("done\n")
for-loop
for (expr1; expr2; expr3) {stmts}
例如:
n = 5for (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)b1+-----+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 = 5i = 1while (i<=n) {printf("i=%g\n", i)i++}printf("done\n")
for-loop
for (expr1; expr2; expr3) {stmts}
例如:
n = 5for (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)b1+-----+1 | 0 |2 | 0 |3 | 0 |+-----+
7.1 declarations
mata中不需要像其他语言中那样声明创建变量的类型,但是mata建议你声明。
element type:transmorphic, numeric, real,complex, string, pointer
organizational type:matrix, vector, rowvector, colvector, scalar
例如之前的zeros()
程序:
real colvector zeros(real scalar c){real colvector aa = 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 AA = J(c, r, 0)return(A)}
mata中的函数还允许设定可选的参数,可选参数通过’|'隔离来设置,同样对于之前的zero()
函数:
real matrix zeros(real scalar c,| real scalar r){real matrix Aif (args()==1) r = 1A = 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 10mata:real matrix zeros(real scalar c,| real scalar r){real matrix Aif (args()==1) r = 1A = J(c, r, 0)return(A)}mata mosave zeros(), replaceend
运行后,会创建一个zeros.mo
文件存储在工作目录下。如果你关闭stata,再打开stata,该函数同样可以运行。
多个函数文件可以存储在.mlib文件中。
8. 使用mata编码stata命令
8.1 一个例子
创建一个ols回归的命令,文件保存至一个.ado文件,即可调用。
program define ols, eclassversion 10.0syntax varlist(numeric) [if] [in]gettoken depvar indepvar: varlistmarksample tousemata: m_ols("‘varlist’", "‘touse’")tempname b Vmatrix ‘b’ = r(b)’matrix ‘V’ = r(V)local N = r(N)matname ‘b’ ‘indepvar’ _cons, c(.)matname ‘V’ ‘indepvar’ _consereturn post ‘b’ ‘V’, depname(‘depvar’) obs(‘N’) esample(‘touse’)ereturn local cmd = "ols"ereturn displayendcapture mata mata drop m_ols()version 10mata:void m_ols(string scalar varlist, string scalar touse){real matrix M, X, Vreal colvector y, breal scalar n, k, s2M = 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))*bs2 = (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
如果觉得《stata中mata语言学习-《Coding with Mata in Stata》》对你有帮助,请点赞、收藏,并留下你的观点哦!