Mata:Mata 编程者指南

Mata:Mata 编程者指南

昨天看完了第二章之后又看了个电影,《正义联盟》,之所以想起来看这部电影,是因为上周五去和丁文亮看了《海王》,昨晚翻豆瓣的时候发现海王之前在正义联盟中出现过,所以就去看了,结果发现正义联盟里面海王的故事似乎和海王里面的不太一样,推测正义联盟的时间早于海底大战。比较好笑的是这部电影里面有个非常好笑的段子,和昨天引用的那个西虹市首富里面的那个有异曲同工之妙。

准备

在进入正文之前,让我们先“四处转转热热身”,和所有的指南一样,我会首先带你看一些你不是很能理解的东西,可能你会无法理解为什么这些很重要。那就先别管它,“Just walk on by”,我们会在稍后分别详述。

Mata 语言中有一些是你熟悉的(类似其它编程语言的),还有一些是你不熟悉的,例如你熟悉的赋值:

Stata
1
a = 0

还有你不熟悉且吃惊的,Mata 允许你同时给多个变量赋值:

Stata
1
a = b = c = 0

Mata 甚至允许你这么做:

Stata
1
a = ((b = c) + 1)

注意上面的(b = c)并不是一个逻辑表达式,而是表示设置b等于c,并且同时设定a = b + 1

另外,不吃惊的是 Mata 可以定义向量和矩阵,但是可能会让你吃惊的是,它们是允许有 0 行、0 列或者 0 行和 0 列的。

下面的部分我们将会交互式的使用 Mata 代码。这或许也是让你吃惊的,因为大多数的编译型语言是不能够交互运行的。交互式的代码使得你能够在使用之前对它们的特性和功能进行实现。

直接显示表达式的结果

Stata
1
2
3
4
5
6
7
8
9
10
11
12
. mata:
------------------------------------------------- mata (type end to exit) ---
: 2 + 2
4

: x = 2 + 2

: x
4

: end
-----------------------------------------------------------------------------

在程序中我们会不经常使用直接表达式(naked expressions),因为printf()函数可以控制更丰富的输出方式,但是直接表达式对于调试代码很有用。例如下面的程序:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
clear all
mata:
mata clear
real scalar nfactorial_over_kfactorial(real scalar n, real scalar k)
{
real scalar result, i
if (n < 0 | n > 1.0x+35 | n != trunc(n)) return(.)
if (k < 0 | k > 1.0x+35 | k != trunc(k)) return(.)
result = 1
for(i = n; i > k; --i) result = result * i
return(result)
}
end
mata: nfactorial_over_kfactorial(5, 3)

知识点:trunc()和int()函数
trunc()和int()函数的功能一样,都是把参数向0的方向进行截断取整,例如int(-5.8) = -5,int(5.2) = 5。需要注意的是int()是ado里面的取整函数,而trunc()是Mata里面的取整函数。

上面函数的运行结果为:

Stata
1
2
. mata: nfactorial_over_kfactorial(5, 3)
20

也就是计算:
$$
\frac{5!}{3!} = 5 \times 4 = 20
$$

为了便于程序的理解,可以在程序里添加提醒:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
mata:
mata clear
real scalar nfactorial_over_kfactorial(real scalar n, real scalar k)
{
real scalar result, i
"计算n!/k!;输入的n和k分别是:"
(n, k)
if(n < 0 | n > 1.0x+35 | n != trunc(n)) return(.)
if(k < 0 | k > 1.0x+35 | k != trunc(k)) return(.)
"开始循环计算:"
result = 1
for(i = n; i > k; --i){
result = result * i
}
"计算完成, 结果为:"
return(result)
}
end
mata: nfactorial_over_kfactorial(5, 3)

运行结果为:

Stata
1
2
3
4
5
6
7
8
9
. mata: nfactorial_over_kfactorial(5, 3)
计算n!/k!;输入的n和k分别是:
1 2
+---------+
1 | 5 3 |
+---------+
开始循环计算:
计算完成, 结果为:
20

赋值

Stata
1
2
mata: greeting = "hello world!"
mata: greeting

多变量同时赋值

Stata
1
2
mata: x = y = 1
mata: x; y

这种同时赋值在初始化中是非常方便使用的,例如:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
mata:
mata clear
x = (1, 2, 3, 4, 5)
y = (2, 3, 4, 5, 6)
z = (3, 4, 5, 6, 7)
xsum = ysum = zsum = 0
for(i = 1; i <= 5; i++){
xsum = xsum + x[i]
ysum = ysum + y[i]
zsum = zsum + z[i]
}
(xsum, ysum, zsum)
end

结果:

Stata
1
2
3
4
5
: (xsum, ysum, zsum)
1 2 3
+----------------+
1 | 15 20 25 |
+----------------+

得益于 Mata 的这种允许多变量同时赋值的特性,你可以写出下面的代码:

Stata
1
2
3
4
5
6
if(avg = (zsum = xsum + ysum)/(zn = yn + xn) == 0){
···
}
else{
···
}

而如果你不使用同时赋值,你的代码将会是下面这个样子的:

Stata
1
2
3
4
5
6
7
8
9
zsum = xsum + ysum
zn = xn + yn
avg = zsum / zn
if(avg == 0){
···
}
else{
···
}

原书作者建议少使用花括号,因为他认为那样会让代码难以阅读,显然如果我们把代码写成这个样子,就会方便阅读很多:

Stata
1
2
3
avg = (zsum = xsum + ysum)/(zn = yn + xn)
if(avg == 0) ···
else ···

实数、复数和字符串

实数

Stata
1
2
3
4
5
mata:
x = 2 + 2
x
sqrt(x)
end

Mata 中的 real 相当于 Stata 的 double,实数的大小顺序为:

$$
非缺失值 < . < .a < .b < … < .z
$$

知识点:Stata中的缺失值
Stata中有27个数值型缺失值,其中`.`是系统缺失值(system missing value)或`sysmiss`,而剩下的26个`.a, .b, .c, ···, .z `被称为扩展缺失值(extended missing values),这些拓展缺失值的作用是帮助使用者更好的跟踪缺失值。 而 Stata 中字符串的缺失值只有一种,也就是""。

复数

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
. mata:
------------------------------------------------- mata (type end to exit) --
: z = (2 + 3i) * (5 - 2i)

: z
16 + 11i

: pi()
3.141592654

: exp(1i * pi())
-1

: end
----------------------------------------------------------------------------

值得注意的是这个表达式:

$$
exp(1i * pi()) = e^{πi} = -1
$$

这是欧拉公式得到的:

$$
\begin{align}
e^{πi} & = cos(π) + sin(π) \
& = -1
\end{align}
$$

字符串(ASCII,Unicode 和 binary)

Stata
1
2
3
4
5
6
7
8
9
10
11
12
. mata:
------------------------------------------------- mata (type end to exit) --
: s = "Mary had a lamb"

: s
Mary had a lamb

: substr(s, 6, 3)
had

: end
----------------------------------------------------------------------------

Mata 的字符串可以是 ASCII,Unicode 和 binary,长度范围 0~20 亿字节。

Mata 拥有常见的所有的字符串处理函数,包括substr()strupper()strlower()等等。

Mata 还拥有 Unicode 字符串函数:usubstr()ustrupper()ustrlower()等等。

除此之外,Mata 还拥有一些格式化输出或者文件 I/O 的函数。

在像 C/C++这样的编程语言中,处理二进制文件需要特别的编程技巧,而在 Mata 中,没什么需要特别处理的,二进制数据可以被存储在字符串变量中并可以使用标准的字符串函数进行处理。例如,substr()函数可以从二进制字符串中提取子字符串,正如它在 ASCII 中做的一样。

下面的代码演示了 Mata 读取数据集auto.dta并显示:

Stata
1
2
3
4
5
6
7
mata:
fh = fopen("auto.dta", "r")
"读取200个字符串"
s = fread(fh, 200)
fclose(fh)
s
end

结果:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
. mata:
------------------------------------------------- mata (type end to exit) --
: fh = fopen("auto.dta", "r")

: "读取200个字符串"
读取200个字符串

: s = fread(fh, 200)

: fclose(fh)

: s
<stata_dta><header><release>118</release><byteorder>LSF</byteorder><K>\uc\
> u0</K><N>J\u0\u0\u0\u0\u0\u0\u0</N><label>\u14\u01978 Automobile Data</lab
> el><timestamp>\u1110 Dec 2018 16:56</timestamp></header><map>\u0\u0\u0\u0\
> u0\u0\u0\u0\xb2\u0\u0\u0\u0\u0\u0\u0-

: end
----------------------------------------------------------------------------

\u0是 Mata 显示二进制 0 的方式。我们还可以发现,s 的长度确实是 200:

Stata
1
2
. mata: strlen(s)
200

而如果我们这样:

Stata
1
2
3
4
5
6
7
. mata: sp = " <stata_dta><header><release>118</release><byteorder>LSF</byte
> order><K>\uc\u0</K><N>J\u0\u0\u0\u0\u0\u0\u0</N><label>\u14\u01978 Automob
> ile Data</label><timestamp>\u1110 Dec 2018 16:56</timestamp></header><map>
> \u0\u0\u0\u0\u0\u0\u0\u0\xb2\u0\u0\u0\u0\u0\u0\u0-"

. mata: strlen(sp)
260

你会发现字符串 sp 的长度是 260!这是因为 Mata 中\u0被视为一个字符,像 abcd 这些一样。

Mata 不把二进制的 0 视为字符串的结尾并且允许字符串的最大长度为 20 亿字节,这些特性是的使用使用 Mata 处理二进制文件变得更加容易。

不管是 ASCII 字符串还是二进制字符串,都能相加或被乘:

Stata
1
2
3
4
5
. mata: "abc" + "def"
abcdef

. mata: 2 * "abc"
abcabc

标量、矢量和矩阵

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
mata:
"行向量"
x = (1, 2, 3)
x
"列向量"
y = (1\2\3)
y
"矩阵"
X = (22, 2.9 \ 17, 3.35 \ 22, 2.64)
X
"矩阵的逆"
invsym(X'X)
end

运行结果:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
. mata:
------------------------------------------------- mata (type end to exit) --
: "行向量"
行向量

: x = (1, 2, 3)

: x
1 2 3
+-------------+
1 | 1 2 3 |
+-------------+

: "列向量"
列向量

: y = (1\2\3)

: y
1
+-----+
1 | 1 |
2 | 2 |
3 | 3 |
+-----+

: "矩阵"
矩阵

: X = (22, 2.9 \ 17, 3.35 \ 22, 2.64)

: X
1 2
+---------------+
1 | 22 2.9 |
2 | 17 3.35 |
3 | 22 2.64 |
+---------------+

: "矩阵的逆"
矩阵的逆

: invsym(X'X)
[symmetric]
1 2
+-------------------------------+
1 | .0182372198 |
2 | -.1225979159 .8617434448 |
+-------------------------------+

: end
----------------------------------------------------------------------------

向量和矩阵被允许拥有 2.81 亿行和列(如果你的电脑内存允许的话)。

函数 rows()、cols()和 length()

Stata
1
2
3
4
5
6
7
8
mata:
"行数"
rows(X)
"列数"
cols(X)
"行数 x 列数(对于向量来说就是长度)"
length(X)
end

这三个函数也能被用在 scalar 上,scalar 是$1×1$的矩阵。

函数 I()

函数 I(n)可以返回一个$n×n$的单位矩阵:

Stata
1
2
3
4
5
6
7
8
. mata: I(3)
[symmetric]
1 2 3
+-------------+
1 | 1 |
2 | 0 1 |
3 | 0 0 1 |
+-------------+

函数 J()

函数 J(r, c, value)可以返回一个$r × c$的矩阵,矩阵元素为 value:

Stata
1
2
3
4
5
6
. mata: J(2, 3, 0)
1 2 3
+-------------+
1 | 0 0 0 |
2 | 0 0 0 |
+-------------+

第三个参数 value 并不一定非要是标量,也可以是个向量或矩阵:

Stata
1
2
3
4
5
6
7
8
. mata: J(4, 2, (3, 4))
1 2 3 4
+-----------------+
1 | 3 4 3 4 |
2 | 3 4 3 4 |
3 | 3 4 3 4 |
4 | 3 4 3 4 |
+-----------------+

Stata
1
2
3
4
5
6
7
8
9
10
11
12
. mata: J(4, 2, (3, 4\2, 1))
1 2 3 4
+-----------------+
1 | 3 4 3 4 |
2 | 2 1 2 1 |
3 | 3 4 3 4 |
4 | 2 1 2 1 |
5 | 3 4 3 4 |
6 | 2 1 2 1 |
7 | 3 4 3 4 |
8 | 2 1 2 1 |
+-----------------+

同样,value 参数也可以是字符串、复数或者其它类型的矩阵:

Stata
1
2
3
4
5
6
. mata: J(2, 2, ("A", "B"))
1 2 3 4
+-----------------+
1 | A B A B |
2 | A B A B |
+-----------------+

行连接和列连接运算符

我们已经知道了逗号用于列连接而反斜线用于行连接,与加减号一样,逗号和反斜线实际上都是一种运算符。

逗号运算被称为列连接(连接列);
反斜线运算被称为行连接(连接行);

在连接运算中逗号运算优先于反斜线运算。例如,a\b,c意味着a\(b, c),而a\(b, c)意味着我要把a作为第一行,而(b, c)作为第二行。

考虑这样的一个表达式:

Stata
1
X = (86, 13 \ 13, 22)

这个表达式中的括号其实是可选择的。加上仅仅是为了易读性。按照上面所说的运算规则,X 中第一行是(86, 13),第二行是(13, 22),所以这是个$2×2$的矩阵。
更多示例:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

. mata:
------------------------------------------------- mata (type end to exit) --
: mata clear

: a = 1, 2

: b = 3, 4, 5

: a, b
1 2 3 4 5
+---------------------+
1 | 1 2 3 4 5 |
+---------------------+

: a, b, a
1 2 3 4 5 6 7
+-----------------------------+
1 | 1 2 3 4 5 1 2 |
+-----------------------------+

: "行连接"
行连接

: a = 1, 2

: b = 3, 4

: A = a\b

: A
1 2
+---------+
1 | 1 2 |
2 | 3 4 |
+---------+

: A \ A
1 2
+---------+
1 | 1 2 |
2 | 3 4 |
3 | 1 2 |
4 | 3 4 |
+---------+

:
: "列连接"
列连接

: B = A, A

: B
1 2 3 4
+-----------------+
1 | 1 2 1 2 |
2 | 3 4 3 4 |
+-----------------+

: end
----------------------------------------------------------------------------

空向量与空矩阵

正如前文所述,Mata 中的向量和矩阵都是可以为空的:

  • J(0, 0, .)可以创建一个实数空矩阵,J(0, 0, 3)也一样;
  • J(0, 0, C(.))可以创建一个复数空矩阵,J(0, 0, 1+2i)也一样;
  • J(0, 0, “”)可以创建一个字符串空矩阵,J(0, 0, “text”)也一样;

需要注意的是,不同类型的空矩阵是不相等的:

Stata
1
2
3
4
5
6
7
8
9
10
11
. mata:
------------------------------------------------- mata (type end to exit) --
: x = J(0, 0, .)

: s = J(0, 0, "")

: x == s
0

: end
----------------------------------------------------------------------------

但是同类型的空矩阵是相等的:

Stata
1
2
3
4
5
6
7
8
9
10
11
. mata:
------------------------------------------------- mata (type end to exit) --
: x = J(0, 0, 1)

: y = J(0, 0, 2)

: x == y
1

: end
----------------------------------------------------------------------------

空向量和空矩阵在处理极值问题中非常有用,例如考虑线性回归的例子,线性回归的逻辑如下:

  • 考虑$y = X \times b$,近似的来说,就是$y = X \times b + e$,在这里$e$表示误差项。
  • X 是一个$r × c$的矩阵,包含 c 个变量的 r 个观测值;
  • y 是一个$r × 1$的向量,包含对应 X 的 r 个因变量的观测值;
  • 那么最小化 e 的平方和的 b 是invsym(X'X)*X'y;

也就是b = invsym(X'X)*X'y会出现在我们的程序中,我们对针对特定的 X 和 y 运行这个程序来产生 b。在前面我们还没有讨论过如何编写 Mata 程序,但是我们可以想象一个包含b = invsym(X'X)*X'y的程序,这个程序大概是下面这个样子的:

Stata
1
2
3
4
5
mata:
X = ···
y = ···
invsym(X'X)*X'y
end

如果我们输入 X 和 y 我们就能计算出b = invsym(X'X)*X'y

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
. mata:
------------------------------------------------- mata (type end to exit) --
: X = (1, 2, 4 \ 3, 4, 6 \ 5, 6, 7)

: y = (3 \ 4 \ 2)

: invsym(X'X)*X'y
1
+--------+
1 | 4 |
2 | -6.5 |
3 | 3 |
+--------+

: end
----------------------------------------------------------------------------

这个结果意味着拟合的模型为:

$$
y = 4x_1 - 6.5x_2 + 3x_3 + ε
$$

当我们真的要去写一个程序去计算这个invsym(X'X)*X'y,我们就要确保这个程序能执行处理一些特殊情况,例如常见的数据缺失问题;在大多数的编程语言中,我们需要写一些额外的代码确保程序自动处理错误,但是 Mata 中我们不用担心这个问题:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
. mata:
------------------------------------------------- mata (type end to exit) --
: X = J(0, 3, .)

: y = J(0, 1, .)

: invsym(X'X)*X'y
1
+-----+
1 | 0 |
2 | 0 |
3 | 0 |
+-----+

: end
----------------------------------------------------------------------------

虽然 X 和 y 都是缺失值,但是程序依然正确处理了,拟合的结果是:

$$
y = 0x_1 + 0x_2 + 0x_3 + ε
$$

显然这个结果是可以接受的,也就是说在 Mata 中我们不需要写程序处理rows(X)=0的问题。
再考虑另外一种特殊情况,cols(X)=0的情况:

Stata
1
2
3
4
5
6
7
8
9
10
. mata:
------------------------------------------------- mata (type end to exit) --
: X = J(5000, 0, .)

: y = J(5000, 1, .)

: invsym(X'X)*X'y

: end
----------------------------------------------------------------------------

可以看出,没有计算出来东西,也就是得到了一个空矩阵,这个时候不会打印出来结果的,你可以按照下面的方法看看这个空矩阵是啥样的:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
. mata:
------------------------------------------------- mata (type end to exit) --
: X = J(5000, 0, .)

: y = J(5000, 1, .)

: b = invsym(X'X)*X'y

: rows(b)
0

: cols(b)
1

: length(b)
0

: end
----------------------------------------------------------------------------

这样就能看到,结果得到的是一个$1×0$的矩阵。

空向量的另外一个用处就是存储循环中产生的结果了,例如下面的 Mata 程序:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
mata:
zvals = J(1, 0, .)
a = b = z = lastz = 1
while(z - lastz < 100){
olda = a; oldb = b
a = olda * oldb
b = olda + oldb
z = a + b
zvals = (zvals, z)
}
zvals
end

运算结果:

Stata
1
2
3
4
5
: zvals
1 2 3 4 5
+-------------------------------+
1 | 3 5 11 41 371 |
+-------------------------------+

空向量的好处就是不会有无用的东西。

Mata 的高级特性

结构体指针就是 Mata 的高级特性,称呼它们为“高级”就意味着它们不是必须的内容。

变量类型

Mata 的变量类型可以是:

  • 实数标量;
  • 字符串向量;
  • 复数矩阵;
  • ···

变量类型由元素类型(elements type)和组织类型(organizational type)组成,这两种东西也被成为eltypeortype,两者之间的任何组合都是可以的,所以下面的也是 OK 的:

  • 字符串标量;
  • 复数向量;
  • 实数矩阵;
  • ···

你应该在你的程序中声明参数、变量和返回值的类型, 例如:

Stata
1
2
3
4
5
6
real colvector bofX(real matrix X, real scalar i)
{
real scalar j, k
real colvector b
···
}

程序中的real colvectorreal matrixreal scalar就是声明,注意在交互式代码运行中声明是不被允许的。

如果变量类型没有被显式地声明,那么它会被被赋予一个默认的格式,也就是被称为:

transmorphic matrix
跨形态矩阵

跨形态 意味着它的eltype是可以被改变的:

Stata
1
2
3
4
5
mata:
x = 4
x = 1i
x = "a"
end

matrix 则意味着它的纬度可以变化:

Stata
1
2
3
4
5
mata:
x = 4
x = (1i, 2i)
x = ("a", "b" \ "c", "d")
end

变量类型这是eltypeortype的组合。

八种eltype是:

Eltype Comment
transmorphic 可以转换成下面的任何一种形态
numeric 可以在实数和复数之间转换
real  
complex  
string  
struct name  
class name  
pointer  

五种ortype是:

Ortype Dimension
matrix $r×c$
vector $1×c$ 或 $r×1$
rowvector $1×c$
colvector $r×1$
scalar $1×1$

结构体

结构体(structures)是自身包含其它变量的对象,结构是通过程序来定义的,例如struct coord

1
2
3
4
5
struct coord
{
real scalar x
real scalar y
}

这个定义就产生了一种新的eltype。一旦定义,你就可以在你的程序中使用这种新的元素类型:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
real scalar foo()
{
struct coord scalar c
real scalar d
···
c.x = ···
c.y = ···
···
d = sqrt(c.x^2 + c.y^2)
···
d = distance_from_origin(c)
···
}

c 是一个struct coord scalar,因为我们就是这样声明它的。实际使用中,很多结构体变量被声明为标量,但是向量和矩阵也是被允许的。

一旦 c 被定义了,你就可以使用 c.x 和 c.y 引用 c 的成员变量。c.x 和 c.y 都是实数标量,因为我们就是这样声明它们的,你可以像使用其它实数变量那样使用 c.x 和 c.y。

在上面的程序中,distance_from_origin()函数的参数是 c,那么这个函数可能是这样的:

Stata
1
2
3
4
real scalar distance_from_origin(struct coord scalar c)
{
return(sqrt(c.x^2 + c.y^2))
}

值得注意的是,当我们把 c 作为参数传递给一个函数的时候,函数接收的是整个 c,c 有两个成员变量,但是它还可以有成百上千个成员,函数可以使用他们中的任何一个。

同样值得注意的是,你也可以写这样的函数:

Stata
1
2
3
4
5
6
7
8
struct coord  scalar unitlength(struct coord scalar c)
{
struct coord scalar toret
real scalar distance
distance = distance_from_origin(c)
toret.c1 = c.c1 / distance
toret.c2 = c.c2 / distance
}

结构体能够被传递和返回正是为什么你要在程序中使用它们的原因,让我们考虑一个包含很多中变量的结构体:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
struct regression_results
{
real vector b
real matrix VCE
real scalar N
real scalar ndf
real scalar ddf
real scalar MSS
real scalar RSS
real scalar R_squared
real scalar s2
}

定义好之后我们就能将这个结构体传递给下面的子程序使用了。如果 r 是一个struct regression_results的话,你就可以传递 r 给子程序,然后子程序就可以获取 r 的所有成员变量了。这些会在后面再进行讨论。

类是广义的结构体,Mata 的类允许共有变量和私有变量、继承、屏蔽、虚函数(virtual functions)和多态(polymorphisms)。

屏蔽(shadowing)
在对象继承过程中,在子类中可以直接使用由父类继承下来的方法和属性;但是如果子类中又声明了相同名称的属性的话,那么当你直接使用的时候就是在使用自己的属性,而不是继承自父类的属性了,这种情形我们称为shadow;如果还是要使用父类的属性,那就需要用super关键字。
Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
class regression
{
public:
real vector b
real matrix VCE
real scalar N
real scalar ndf
real scalar ddf
real scalar R_squared
real scalar s2
private:
real scalar MSS
real scalar RSS
public:
void calc()
private:
real matrix XX()
real colvector Xy()
}

公有变量的使用和结构体中的成员变量是一样的,类的成员函数的调用也是使用点运算符,如果 c 类中有公有函数 calc(),那么调用这个函数的方法就是c.calc()

而私有变量和私有函数只能被类中的成员函数使用。

要定义一个一个类,我们需要先声明它,然后再编写它的成员函数:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
class regression
{
public:
real vector b
real matrix VCE
real scalar N
real scalar ndf
real scalar ddf
real scalar R_squared
real scalar s2
private:
real scalar MSS
real scalar RSS
public:
void calc()
private:
real matrix XX()
real colvector Xy()
}

void regression::calc(real colvector y, real matrix X)
{
real matrix XXinv
XXinv = invsym(XX(X))
b = XXinv * Xy(X, y)
RSS = sum((y - X * b):^2)
s2 = RSS / (rows(X) - cols(X))
···
···
}
real matrix regression::XX(real matrix X)
{
return(X'X)
}
real matrix regression::Xy(real colvector y, real matrix X)
{
return(X'y)
}

注意这里面用到了分号运算符:^,分号运算符会分别作用于一个向量或矩阵中的每一个元素,例如:

$$
e : *
\left[
\begin{matrix}
a & b \
c & d
\end{matrix}
\right] =
\left[
\begin{matrix}
a * e & b * e \
c * e & d * e
\end{matrix}
\right]
$$

分号运算符的详细介绍可以参考:help [M-2] op_colon

需要注意的是,成员函数calc()使用成员bRSSs2XX()Xy()并没有使用c.前缀,这是因为在类函数的内部,成员是可以直接调用的。而在成员函数的外部,你必须使用varname.membername的方式进行调用:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
void myreg(real colvector y, real matrix X)
{
class regression scalar r
real scalar i
r.calc(y, X)
printf("s2 = %11.0g\n", r.s2)
"系数向量"
for(i = 1; i <= length(r.b); i ++){
printf("x%g %10.0g\n", i, r.b[i])
}
···
···
}

关于类的内容会在后面的章节详述。

指针

指针(pointers)是非常邪门的,它们似乎总是出现在复杂的程序中。但是这些复杂的程序的复杂并不是因为它们使用了指针,指针本身只非常容易理解的,它们可以解决一些其它东西解决不了的问题。

假设你使用线性回归类写了一个线性回归方程,如果你完全开发好了这个类,你会很快发现把矩阵 X 的一个拷贝放在这个类里面是很方便的,但是因为内存问题,你会犹豫这么做。如果 X 包含 100 个变量的 10 万个观测值,那么这个矩阵会占用 3.8 亿字节。那么你会觉得把一个 3.8 亿字节的东西创建一个拷贝放入你的类中是一件方便的事情么?这个时候使用指针就会很方便,它只会在内存中占有 8 个字节。

指针是一个包含了另一个变量的内存地址的变量。如果你输入下面的代码:

Stata
1
2
3
mata:
p = &a
end

那么 p 就会包含 a 的内存地址,&前缀是 Mata 中的地址运算符号,例如:

Stata
1
2
3
4
5
mata:
a = 2
p = &a
p
end

结果:

Stata
1
2
3
4
5
6
7
8
9
10
11
. mata:
------------------------------------------------- mata (type end to exit) --
: a = 2

: p = &a

: p
0x7f8503db8098

: end
----------------------------------------------------------------------------

0x7f8503db8098就是变量 a 存储的十六进制的地址。我们在乎的不是这个地址,而是 p 指向&a。

&a 的含义就是 a 的内存地址。

*p 的含义是 p 中的地址上存储的内容:

Stata
1
2
. mata: *p
2

是的,的确,*p就是 2,也就是说*p是 a 的同义词,如果改变*p的值,a 的值同时也会改变:

Stata
1
2
. mata: *p = 4; a
4

可以看到,a 的值也变成了 4。

最后一个需要注意的事情是,指针不能指向任何包含NULL的东西,这里的NULLnull vectorsnull matries中的null是不一样的,大写的NULLnull memory address,即空内存地址。

现在 p 指向 a,我们做下面的改变:

Stata
1
2
3
4
5
6
mata:
p
p = NULL
p
*p
end

运行结果:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
. mata:
------------------------------------------------- mata (type end to exit) --
: p
0x7f8503db8098

: p = NULL

: p
0x0

: *p
<istmt>: 3010 attempt to dereference NULL pointer
----------------------------------------------------------------------------
r(3010);

. end
command end is unrecognized
r(199);

你会发现,结果出错了。再后面会对指针的使用有更详细的介绍。

注意事项

如何交互式的编写 Mata 代码

运行mata:进入 Mata 模式,然后即可逐行运行代码。
让我们试着求一个满秩矩阵的逆:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
. mata:
------------------------------------------------- mata (type end to exit) --
: X = (86, 13 \ 13, 22)

: Xinv = invsym(X)

: Xinv * X
[symmetric]
1 2
+---------+
1 | 1 |
2 | 0 1 |
+---------+

: end
----------------------------------------------------------------------------

invsym()函数也能作用于奇异矩阵(非满秩矩阵,但是必须是方阵)。

Stata
1
2
3
4
5
mata:
X = (86, 13, 0 \ 13, 22, 0 \ 0, 0, 0)
Xinv = invsym(X)
Xinv * X
end

运算结果:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
. mata:
------------------------------------------------- mata (type end to exit) --
: X = (86, 13, 0 \ 13, 22, 0 \ 0, 0, 0)

: Xinv = invsym(X)

: Xinv * X
[symmetric]
1 2 3
+-------------+
1 | 1 |
2 | 0 1 |
3 | 0 0 0 |
+-------------+

: end
----------------------------------------------------------------------------

需要注意的是invsym()函数只能求对称矩阵的逆,如果求一个非对称方阵的逆:

Stata
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
. mata:
------------------------------------------------- mata (type end to exit) --
: X = (1, 2 \ 3, 5)

: Xinv = invsym(X)

: X * Xinv
1 2
+-----------+
1 | 1 0 |
2 | 5 -1 |
+-----------+

: end
----------------------------------------------------------------------------

你会发现你得到了一个错误的结果。这就提醒我们,当我们想要求一个矩阵的逆的时候,我们应该首先检查这个矩阵的对称性,Mata 中的issymmetric()函数可以完成这一操作:

Stata
1
2
. mata: issymmetric(X)
0

结果为 0 表示不对称。但是这种检查并不总是必要的,例如你想求X'X的逆,这一操作就是没有意义的,因为X'X一定是对称的。

_error():中止函数

_error()是 Mata 中产生回溯日志和中止程序的执行。我们可以在使用该函数进行手动抛出异常并中止程序:

Stata
1
2
3
mata:
if(situation irretrievable) _error(...)
end

大部分时候我们是不需要这么做的,因为 Mata 会自行抛出异常。

_error()有三种用法:

  • _error(#)
  • _error(“string”)
  • _error(#, “string”)

#表示错误代码,"string"表示错误日志。你可以使用下面的方法检查错误代码的内容:

Stata
1
2
3
4
5
6
7
8
9
10
11
. mata:_error(3300)
<istmt>: 3300 argument out of range
r(3300);

. mata:_error("argument N out of range")
<istmt>: 3498 argument N out of range
r(3498);

. mata:_error(3300, "argument N out of range")
<istmt>: 3300 argument N out of range
r(3300);

# Mata

评论

程振兴

程振兴 @czxa.top
截止今天,我已经在本博客上写了165.1k个字了!

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×