使用 Stata 寻找分布中密度为最大密度一半的 X

使用 Stata 寻找分布中密度为最大密度一半的 X

这是之前我的一个小伙伴问我的问题,她有一个数据(附件中的 A.csv)这个数据是一个 30x10000 的矩阵(也就是 30 行 10000 列),Stata 的 xpose 命令可以用于数据集的转置,不过这个命令使用的时候对变量数量有限制,经过我测试,变量不能超过 9999,所以这里我们删除最后一个变量:

1
2
3
4
5
6
7
8
9
/* 把 CSV 文件转换成 dta 文件 */
clear all
set maxvar 10000
import delimited using A.csv, clear
* 删除最后一个变量
drop v10000
* 转置
xpose, clear
save A, replace

转置之后得到的是个 9999x30 的数据集。我们的这个小伙伴的问题是这样的,这些变量的分布大致都服从正态分布,然后她想找到变量取何值的时候对应的密度为最大密度的二分之一。

我们先看变量 v1:

1
2
3
use A, clear
* h kdensity
kdensity v1, generate(nv1 nd1) n(10000)

kdendity 命令的 generate 可以生成核密度拟合得到的分布曲线的上的 X 值和密度值。选项 n 表示生成的观测值数量,我们可以用生成的 nv1 和 nd1 自行绘制这幅图:

1
tw line nd1 nv1

这 10000 个 nd1 和 nv1 就是格点数据,下面我们就是利用这个格点数据寻找我们想找到的 “Semi-max X”。

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

这两个点正是我们需要的!下面我们做个循环就可以得到所有的 30 个变量的这四个值了,注意到刚刚我们只有四个值,我使用了四个 local 变量存储,下面我使用的是个 30x6 的矩阵存储,多的两个是最大密度值及其对于的 X:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

下面我们还可以画个图观察下我们得到的这些值是否准确:

X 值:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

密度值:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

那么问题又来了!图例哪儿去了?我用的这个主题似乎会默认把图例 off 掉。

这样就可以出来了:

代码去哪了?

代码可以加入我的知识星球后从知识星球下载附件获取~
要了解如何加入我的知识星球,可以阅读关于界面或者添加我的微信咨询。

通过这个案例大家应该就能理解如何在循环中保存运行结果和使用这些结果了。实际上很多 Stata 命令会把打印在屏幕上的结果也保存着返回值中,这些返回值可以运行 ret list 或者 eret list 查看。

知识星球附件链接:https://t.zsxq.com/Zn2BIuf

#

评论

Your browser is out-of-date!

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

×