听说会Stata的人，数学不会太差？
本文作者：张馨月
文字编辑：王碧琪
技术总编：李朋冲
爬虫俱乐部将于2020年1月5日至11日在湖北武汉举行为期一周的Stata编程技术定制培训，此次采取初级班和高级班分批次培训。课程通过案例教学模式，旨在帮助大家在短期内掌握Stata软件编程、金融计量知识和实证分析方法，使大家熟悉Stata核心的爬虫技术，以及Stata与其他软件交互的高端技术。目前正在火热招生中~详细培训大纲及报名方式，请点击《爬虫俱乐部2020第一期Stata编程训练营开始报名啦！》或点击文末阅读原文呦~
另外，2019年11月29日-12月1日在湖北武汉举办的《2019年Python第四期培训招生》，招生工作已经结束，请大家继续关注我们后续的培训公告哦~
在之前分享的一篇Stata解数学题的推文中，许多小伙伴表示题目太简单没有挑战性。这不，有一个宝妈在给娃辅导作业时，遇到了下面一道数学题，终于难倒了身边的初中生高中生大学生研究生：
那么这道题在Stata中有没有什么好的解法呢？在开始讲解之前，小编再来问一个简单的问题：
大家都知道圆周率的取值是3.1415927，但这个数值究竟是如何计算出来的？
一个简单的方法就是向正方形内投点，因为对于一个正方形，它的内切圆与正方形的面积之比是π/4。那么，在这个正方形内部，随机产生N个点（这些点服从均匀分布），计算各点与原点的距离，判断落在圆的内部点的个数n，可以得到圆内的点占比为n/N= π/4。因此，就可以求出π的值了~
如果我们只随机投放100个点，计算误差会非常大。但是，当投放的数量为1000,10000,20000，计算的π值会不会越来越精准呢？
对于这道题目，我们也可以采用相同的思路。如果对1-6随机排序一次，然后计算绝对值之和，我们无法确定得到的结果就是最小值。但是将这样的过程重复100000次寻找最小值，是不是就可以确定所得到的答案了呢？
为了实现这一目的，我们使用Stata，按照下面的思路求解：
1.生成变量x，取值分别是1，2，……，6；
2.打乱x的排序，分别放入变量x1、x2、x3、x4、x5、x6中并计算题目中要求式子的值；
3.重复以上过程100000次。
下面我们就来动手操作一下吧~
1.生成变量
我们生成一个容量为6的样本，并将每行的行号赋值给变量x，就能完成第一步：
clear allset obs 6 gen x = _n
2.对x1-x6赋值并计算绝对值
首先我们需要完成对x的随机排序，一个思路是生成一组随机数，并将数据集根据该组随机数来排序：
gen ran=uniform()sort ran 运行结果如下：
接下来，如何将x的取值赋给x1-x6呢？如果在同一个数据集中生成变量x1-x6并重复100000次，可以想象工作量是非常庞大的。比较理想的解决方案是生成一个新的数据集，将x的排序结果不断传递到该数据集，这一想法可以通过frame来实现。
frame可以在内存中创建多个框架（helplimits可看到最多是100个），用户可以在它们之间自由切换，并将某个数据链接到其他数据集中。关于它的详细用法，大家可以在推文《Stata16新功能之“框架”——读入多个数据集（1）》《Stata16新功能之“框架”——基础命令大合集（2）》进一步了解。
下面我们来具体应用一下：
frame create abs x1 x2 x3 x4 x5 x6 abs //创建一个框架命名为abs，该数据集中包含的变量为x1 x2 x3 x4 x5 x6 abs，变量定义如下local x1 = x[1]local x2 = x[2]local x3 = x[3]local x4 = x[4]local x5 = x[5]local x6 = x[6]local abs = abs(`x1'-`x2') + abs(`x2'-`x3') + abs(`x3'-`x4') /// + abs(`x4'-`x5') + abs(`x5'-`x6') + abs(`x6'-`x1')frame post abs (`x1') (`x2') (`x3') (`x4') (`x5') (`x6') (`abs') //将每一个“宏”里的信息邮寄到数据集abs中的相应变量通过这一步操作，我们就可以把本次随机模拟的结果储存在数据集abs当中了，结果如下：
3.重复上述过程100000次
为了让程序重复进行，我们可以用forvalue把这段命令放入循环当中，但此时程序会报错：
这是因为我们在循环中定义了变量ran，当循环进行到第二层时，Stata就会提醒我们这个变量已经定义过了，后续程序也就无法进行。这里我们就要用到Stata的一大法宝——起死回生命令preserve和restore，这样程序就可以顺利运行啦：
clear allset obs 6 gen x = _nframe create abs x1 x2 x3 x4 x5 x6 abs forvalue i = 1/100000{preserve gen ran=uniform()sort ran local x1=x[1]local x2=x[2]local x3=x[3]local x4=x[4]local x5=x[5]local x6=x[6]local abs = abs(`x1'-`x2') + abs(`x2'-`x3') + abs(`x3'-`x4') /// + abs(`x4'-`x5') + abs(`x5'-`x6') + abs(`x6'-`x1')frame post abs (`x1') (`x2') (`x3') (`x4') (`x5') (`x6') (`abs')restore }我们将内存中的数据集切换至abs，得到的数据部分截图如下：
cwf abs //切换框架到abs，等同于frame change再运用命令sum，就可以得到100000次模拟中abs的描述性统计结果：
sum abs 结果显示abs的最小值为10，最大值为18。因此，题目的答案为10~
如果我们还想知道abs在10-18之间取每个值的概率为多少，可以再进行如下操作：
local N = `r(N)' //r(N)为sum的返回值forvalues i = `r(min)'/`r(max)'{ count if abs == `i' //计算abs在10-18之间取每个值的次数 di "取值为`i'的概率为:`=`r(N)'/`N''" //r(N)为count if的返回值}从结果可以看到abs的取值均为偶数，为奇数的概率为0，感兴趣的小伙伴可以思考一下背后的原因是什么~
完整程序如下：
clear allset obs 6 gen x = _nframe create abs x1 x2 x3 x4 x5 x6 abs forvalue i = 1/100000{preserve gen ran=uniform()sort ran local x1=x[1]local x2=x[2]local x3=x[3]local x4=x[4]local x5=x[5]local x6=x[6]local abs = abs(`x1'-`x2') + abs(`x2'-`x3') + abs(`x3'-`x4') /// + abs(`x4'-`x5') + abs(`x5'-`x6') + abs(`x6'-`x1')frame post abs (`x1') (`x2') (`x3') (`x4') (`x5') (`x6') (`abs')restore }
cwf abssum abslocal N = `r(N)' forvalues i = `r(min)'/`r(max)'{ count if abs == `i' di "取值为`i'的概率为:`=`r(N)'/`N''" }
事实上，无论是文章开头提到的计算π值，还是这道题目的求解过程，都是蒙特卡洛模拟的初步思想。
蒙特卡罗方法的原理是通过大量随机样本，去了解一个系统，进而得到所要计算的值。具体来讲，蒙特卡洛模拟是以概率和统计的理论、方法为基础的一种数值计算方法，将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解。小到求解小学初中高中的数学问题，大到解决期权定价、人工智能领域以及其他的随机模拟问题，这种方法都有非常广泛的应用。在之后的推文中，我们将继续带领大家探索蒙特卡洛模拟~
关于我们
微信公众号“Stata and Python数据分析”分享实用的stata、python等软件的数据处理知识，欢迎转载、打赏。我们是由李春涛教授领导下的研究生及本科生组成的大数据处理和分析团队。
1）必须原创，禁止抄袭；
2）必须准确，详细，有例子，有截图；
注意事项：
1）所有投稿都会经过本公众号运营团队成员的审核，审核通过才可录用，一经录用，会在该推文里为作者署名，并有赏金分成。
2）邮件请注明投稿，邮件名称为“投稿+推文名称”。
3）应广大读者要求，现开通有偿问答服务，如果大家遇到有关数据处理、分析等问题，可以在公众号中提出，只需支付少量赏金，我们会在后期的推文里给予解答。