糖尿病康复,内容丰富有趣,生活中的好帮手!
糖尿病康复 > stata中mata语言学习-《Coding with Mata in Stata》

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

时间:2023-01-28 21:37:42

相关推荐

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

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 |+----------------------+

取出内存中已载入数据集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_datast_view两个函数的区别。

st_datast_view两个函数所实现的功能相似,但是,st_data运算更加快速,且不改变stata内存中的数据。

注意:对于字符串变量需使用st_sdatast_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》》对你有帮助,请点赞、收藏,并留下你的观点哦!

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。