名师讲堂｜专利的前向相似度与后向相似度：使用 Stata 计算专利关键数字技术突破创新指标（二）

上次课程主页：https://rstata.duanshu.com/#/course/bcb7f4d90a8445e48c18e4ed769ef071

继续上次课的内容，今天我们继续讲解第二部分：

专利的前向相似度与后向相似度：使用 Stata 计算专利关键数字技术突破创新指标（二）：基于 Mata 程序的专利间前瞻、回溯相似度计算；

tempdata7.dta 是上次课的处理结果：

use tempdata7.dta, clear

一个简单的例子

在开始计算之前，我们先考虑一个简单的例子。

data 数据包含 9 个变量，第一列是专利的编号，第二列是专利的申请年份，后面 7 列就是专利的 7 个属性维度，我们来计算每个专利的前瞻和回溯相似度：

mata:
mata clear
// 假设数据已经加载到 Mata 中
data = (1, 2000, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7
        2, 2001, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8
        3, 2002, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9
        4, 2001, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0)

// 提取专利编号、申请日和属性
patent_id = data[., 1]
application_date = data[., 2]
attributes = data[., 3..cols(data)]

// 第 3 个专利
// 早于该专利的
index_p = selectindex(application_date :< application_date[3])
index_p

// 这些专利的属性矩阵
M = attributes[index_p, .]

M

sqrt(attributes[3, .] * attributes[3, .]')

// 余弦相似度

attributes[3, .] * M'
attributes[3, .] * attributes[3, .]'
M * M'

// 提取对角线元素再开根号
sqrt(diagonal(M * M')')

// 计算余弦相似度
cos = (attributes[3, .] * M') :/ (sqrt(attributes[3, .] * attributes[3, .]') :* sqrt(diagonal(M * M')'))
rowsum(cos)

// 晚于该专利的
index_l = selectindex(application_date :> application_date[1])
index_l

// 这些专利的属性矩阵
M = attributes[index_l, .]

M

sqrt(attributes[1, .] * attributes[1, .]')

// 余弦相似度

attributes[1, .] * M'
attributes[1, .] * attributes[1, .]'
M * M'

// 提取对角线元素再开根号
sqrt(diagonal(M * M')')

// 计算余弦相似度
cos2 = (attributes[1, .] * M') :/ (sqrt(attributes[1, .] * attributes[1, .]') :* sqrt(diagonal(M * M')'))
rowsum(cos2)
end

计算原理非常简单，就是分别提取申请时间早于和晚于该专利的，使用余弦相似度公式进行计算即可。

然后我们就可以循环计算简单示例中所有专利的了：

mata:
mata clear
// 假设数据已经加载到 Mata 中
data = (1, 2000, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7
        2, 2001, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8
        3, 2002, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9
        4, 2001, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0)

// 提取专利编号、申请日和属性
patent_id = data[., 1]
application_date = data[., 2]
attributes = data[., 3..cols(data)]

// 创建文件保存结果
stata("cap erase res1.csv")
filename = "res1.csv"
fh = fopen(filename, "rw")
for(i = 1; i <= rows(data); i++) {
    index_p = selectindex(application_date :< application_date[i])
    index_l = selectindex(application_date :> application_date[i])

    M_p = attributes[index_p, .]
    M_l = attributes[index_l, .]

    cos_p = (attributes[i, .] * M_p') :/ (sqrt(attributes[i, .] * attributes[i, .]') :* sqrt(diagonal(M_p * M_p')'))
    cos_l = (attributes[i, .] * M_l') :/ (sqrt(attributes[i, .] * attributes[i, .]') :* sqrt(diagonal(M_l * M_l')'))
    fput(fh, invtokens(strofreal((patent_id[i], application_date[i], rowsum(cos_p), rowsum(cos_l)), "%16.12g"), ","))
}
end

处理结果如下：

import delimited using res1.csv, clear
*> (encoding automatically selected: UTF-8)
*> (4 vars, 4 obs)
list
*>      +---------------------------------+
*>      | v1     v2         v3         v4 |
*>      |---------------------------------|
*>   1. |  1   2000          0   2.969368 |
*>   2. |  2   2001   .9965457   .9982744 |
*>   3. |  3   2002    2.98728          0 |
*>   4. |  4   2001   .9828722   .9990562 |
*>      +---------------------------------+

基于真实的专利数据

然后我们就可以把这个代码应用到真实的专利数据了，首先尝试 100 条专利的计算：

use tempdata7, clear
*> (数据处理：微信公众号 RStata)
keep in 1/100
*> (586,192 observations deleted)
destring newipzlid, replace
*> newipzlid: all characters numeric; replaced as long

*- 此处代码需下载讲义材料查看～

import delimited using res2.csv, clear
list in 1/10

这里仅仅演示了 100 条的计算结果，不过想要靠这段代码计算所有的是不太可能的（一个几十万 x 几十万维度的矩阵计算是非常耗时的）。所以下面我们得思考下如何提升计算效率。

观察 tempdata7.dta 的数据可以发现，其实大多数专利的 7 个指标都是相同的，互不相同的并不多：

use tempdata7, clear
*> (数据处理：微信公众号 RStata)
duplicates drop 人工智能技术 - 高端芯片技术, force
*> Duplicates in terms of 人工智能技术 元宇宙技术 区块链技术 工业互联网技术 物联网技术 量子信息技术 高端芯片技术
*> (585,755 observations deleted)
di `=_N'
*> 537

那么在计算余弦相似度的时候，我们实际上并不需要计算所有的，只需要计算互不相同的，然后再根据数量加权求和即可。下面我们再来改进下代码。

为了更好的分辨哪些是相等的，我们再生成一个 group 变量，group 值相等的就是完全一样的：

use tempdata7, clear
*- 如果 group 相等，那么就是同样的属性
egen group = group(人工智能技术 元宇宙技术 区块链技术 工业互联网技术 物联网技术 量子信息技术 高端芯片技术)

keep in 1/100
destring newipzlid, replace
mata:
data = st_data(., .)
// 提取专利编号、申请日和属性
patent_id = data[., 1]
application_date = data[., 2]
attributes = data[., 3..cols(data)]

i = 2
index_p = selectindex(application_date :< application_date[i])
index_l = selectindex(application_date :> application_date[i])

M_p = attributes[index_p, .]
M_l = attributes[index_l, .]

// 互不相同的行
M_p2 = uniqrows(M_p)
M_l2 = uniqrows(M_l)

// 各行的数量（需要循环统计）
freq_p = J(rows(M_p2), 1, .)
for (j=1; j<=rows(M_p2); j++) {
    freq_p[j] = colsum(M_p[., 8] :== M_p2[j, 8])
}

freq_l = J(rows(M_l2), 1, .)
for (j=1; j<=rows(M_l2); j++) {
    freq_l[j] = colsum(M_l[., 8] :== M_l2[j, 8])
}

// 去除最后一列
M_p2 = M_p2[., 1..7]
M_l2 = M_l2[., 1..7]
attributes = attributes[., 1..7]

cos_p = (attributes[i, .] * M_p2') :/ (sqrt(attributes[i, .] * attributes[i, .]') :* sqrt(diagonal(M_p2 * M_p2')'))
cos_l = (attributes[i, .] * M_l2') :/ (sqrt(attributes[i, .] * attributes[i, .]') :* sqrt(diagonal(M_l2 * M_l2')'))
rowsum(cos_p :* freq_p')
rowsum(cos_l :* freq_l')
end

可以看到这里的计算结果和前面的结果是一样的。由于时间关系，我们还是只选择前 100 个运行：

然后就可以循环所有的：

*- 此处代码需下载讲义材料查看～

import delimited using res2.csv, clear
*> (encoding automatically selected: UTF-8)
*> (4 vars, 100 obs)
list in 1/10
*>      +------------------------------------------+
*>      |         v1      v2         v3         v4 |
*>      |------------------------------------------|
*>   1. | 2010000022   18406          1   7.352787 |
*>   2. | 2010000023   18414   4.541115   3.811672 |
*>   3. | 2010000028   18443   24.91969          1 |
*>   4. | 2010000029   18408    12.7567   12.94932 |
*>   5. | 2010000030   18426   21.95129    4.02529 |
*>      |------------------------------------------|
*>   6. | 2010000031   18409   14.02726   11.94932 |
*>   7. | 2010000033   18452   25.97658          0 |
*>   8. | 2010000042   18408          2    6.08223 |
*>   9. | 2010000043   18339          0   8.352787 |
*>  10. | 2010000044   18438   7.811672   .5411147 |
*>      +------------------------------------------+

不过如果想要运行所有的还是比较耗时的，下次课我们再讲解如何改用并行的方法。

上面的所有计算都是针对之前和之后所有专利的，按照 Kelly et al.（2021）中的介绍，我们可以限制计算之前 5 年的，例如这里我们可以仅仅计算上个日历年的。这里需要注意，如果选择计算申请日前 365 天的，可以在 Mata 里面使用下面的代码：

... Other codes ...
index_p = selectindex(application_date[i] :- application_date :<= 365)
... Other codes ...

如果想要计算上个日历年的，就需要在 Stata 数据里面进行调整下了：

use tempdata7, clear
*- 抽取 2% 的样本
set seed 1234
sample 0.02
*- 把申请日变成年份
replace 申请日 = yofd(申请日)
format 申请日 %6.0f
*- 如果 group 相等，那么就是同样的属性
egen group = group(人工智能技术 元宇宙技术 区块链技术 工业互联网技术 物联网技术 量子信息技术 高端芯片技术)

destring newipzlid, replace
mata:
data = st_data(., .)
// 提取专利编号、申请日和属性
patent_id = data[., 1]
application_date = data[., 2]
attributes = data[., 3..cols(data)]

stata("cap erase res3.csv")
filename = "res3.csv"
fh = fopen(filename, "rw")
for(i = 1; i <= rows(data); i++) {
    i
    index_p = selectindex(application_date[i] :- application_date :>= 1)
    index_l = selectindex(application_date :- application_date[i] :>= 1)

    M_p = attributes[index_p, .]
    M_l = attributes[index_l, .]

    // 互不相同的行
    M_p2 = uniqrows(M_p)
    M_l2 = uniqrows(M_l)

    // 各行的数量（需要循环统计）
    freq_p = J(rows(M_p2), 1, .)
    for (j=1; j<=rows(M_p2); j++) {
        freq_p[j] = colsum(M_p[., 8] :== M_p2[j, 8])
    }

    freq_l = J(rows(M_l2), 1, .)
    for (j=1; j<=rows(M_l2); j++) {
        freq_l[j] = colsum(M_l[., 8] :== M_l2[j, 8])
    }

    // 去除最后一列
    M_p3 = M_p2[., 1..7]
    M_l3 = M_l2[., 1..7]
    attributes3 = attributes[., 1..7]

    cos_p = (attributes3[i, .] * M_p3') :/ (sqrt(attributes3[i, .] * attributes3[i, .]') :* sqrt(diagonal(M_p3 * M_p3')'))
    cos_l = (attributes3[i, .] * M_l3') :/ (sqrt(attributes3[i, .] * attributes3[i, .]') :* sqrt(diagonal(M_l3 * M_l3')'))
    fput(fh, invtokens(strofreal((patent_id[i], application_date[i], rowsum(cos_p :* freq_p'), rowsum(cos_l :* freq_l')), "%16.12g"), ","))
}
end

import delimited using res3.csv, clear

这样我们就实现了基于指定窗口范围的计算。

下次课程中我们再讲解如何使用并行运算提升速度。

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜专利的前向相似度与后向相似度：使用 Stata 计算专利关键数字技术突破创新指标（二）