这是这个星期（准确来说是这个星期和下个星期）的机械工学实验的题目——嘛因为是机械工学不是情报工学（信息工程学）所以本来来说C语言神马的本来并不是偶们的必修，8过这一次的主题是情报处理，于是考虑到绝大部分同学应该根本就从来没编过程序所以才要用两个星期来处理这个的吧w（于是好歹还是会那么一点C的偶表示优越感233
好了首先介绍一下这次的实验课题：使用蒙地卡罗方法来求得近似的圆周率。什么是蒙地卡罗方法呢？简单来说就是通过无数的随机函数来确定一个我们想要的东西的近似值的方法。因为是通过随机函数的方式，所以自然做的次数越多则近似值越精确
这样说估计大家还是不太听得懂——没关系，这里偶会使用下面的图片来解释这次的实验课题，顺便大家也可以通过这个课题的做法来大致了解蒙地卡罗方法w
嗯，题目就是这样——8过难得这次的题目（对于偶来说）这么方便（最重要的是可以用电脑不用手写实验报告XD），所以今天的博客偶稍微把题目改动一下——改为求证实验次数的增加对于蒙地卡罗法来说是否能提高精确度
实验方法非常简单，根据上述的求圆周率的方法来求得圆周率，每一次实验得出的圆周率都将减去实际的圆周率（C++提供了圆周率函数），并取其绝对值就能得出每一次实验的计算误差值；再将这个误差值作为Y轴、实验的次数作为X轴，两轴均取对数描图，就能得出随实验次数的增加，误差值会如何变化的图像了
这里作为例子，仅设置一共实验10000次；需要用到的试验工具有：
- 能够写代码的工具（比如文本编辑或者记事本之类）
- gcc（或者其他的任何可以编译C++的环境）
- GNUPLOT（描图工具，跨平台8过偶没用过Windows版不知道怎么用w总之先确保自己的机子内安装好GNUPLOT吧）
以上；8过因为例题也好偶自己的环境也好都是UNIX/Linux环境（准确来说偶的是Mac OS X内的UNIX环境w），所以不确定Windows环境的实验进行方式，有兴趣的童鞋欢迎自行捣鼓w
好了，试验开始：首先使用C++编写一个利用蒙地卡罗法球圆周率的程序（虽然使用文本编辑也可以，8过为了大家看着方便这里用Xcode进行编写，好处是根据代码的不同含义能自动有不同的配色，大家看起来可能比较方便理解=w=a），将这个程序保存在自己喜欢的地方——比如偶保存在了桌面；文件名也可以用大家喜欢的名称，只要后缀名是.cpp就行（其实如果只在意是否能编译成功的话后缀名不用太在意也没太大关系的本来，只不过为了突出显示这是一个C++代码文件所以用.cpp的后缀名而已w），比如偶这里用的是pi.cpp
接下来在Mac上打开Terminal（终端），前往刚才保存到的文件夹的路径（比如偶的就在桌面），使用C++编译代码文件（编译后程序名称也可以随便弄，比如偶这里用的pi），执行编译好的程序并将执行结果写入到新文本（文本名称也随便弄，比如偶这里是pi.dat），到这里文件的编译就算完成了，接下来是如何将写出的文本制成图像
仍然还是在终端里面，使用gnuplot命令（当然前提是GNUPLOT已经安装好了）进入GNUPLOT
在GNUPLOT里面，首先设置X轴的标签为i（次数），然后设置Y轴的标签为error（误差），接下来是两个轴均使用对数坐标轴，并设置Y轴的范围是从1到0.00001（也就是Y轴下面是1上面是0.00001，因为这样设置出来的话图像会随着误差的降低而升高，看上去有种精度提高了的感觉www当然这只是偶的个人喜好，不喜欢的完全可以不设置Y轴范围，这样Y坐标轴就是默认的越往上走数值越大），这些都设置好了之后就可以载入刚刚写出的结果文本描图了；如果不指定的话默认是画点上去，8过如果在文本后面指定with lines [1]的话就可以使用直线描图了
比如，这里这个执行了10000次的结果描出来的图就是这样的感觉——虽然中途的确有几次距离真实的圆周率非常接近了，8过除去那样的低概率事件，大体上来说的确是实验次数越多则精确度也越高
顺便一说，如果有兴趣的话大家也可以点击编译好的程序（比如偶这里是pi）运行看看结果是什么——其实就是每一次的误差值了啦w
最后——这个只建议家里电脑性能超好而且自己超级闲得找不到事干了的人做——如果在编写C++代码的时候将n值设的很高的话，那么编译并执行程序需要花很长的时间不说，使用GNUPLOT载入如此庞大的数据量制图也会非常花时间——8过结果倒是的确非常有趣，比如这里偶设置了一次n值为100000000，也就是做了100000000次实验的代码，于是从开始编译执行到最终制成图像似乎花了一共有将近半个小时的时间（可怜了偶的女仆本，执行这些命令的时候偶还是一如既往开着RSS、邮件、Twitter、浏览器、IM聊天工具甚至Xcode等东西wwwwww），于是制出来的图像的确更能体现出随着实验次数的增加精确度也在增加的倾向——注意这图和前面那个10000次相比Y轴在小数点后面更多了1位哦w
嘛，以上就是本次实验的课题外加本人的一点追加wwww
DT..
完全可以使用Xcode写个Mac程序出来还带GUI比那些命令行都好（前提是有时间的话..
各种天书。。。o(╯□╰)o囧
———-
置百丈玄冰而崩裂，掷须臾池水而漂摇。
没看懂由计数器相除乘4得到pi值？
pi = 4*(double)a/(double)i;
这里！求讲解，我几何肯定不行，C++略懂皮毛
吐槽点好多
比如同时用了C++和C的头文件声明格式
再比如用了一个常数来初始化随机数发生器
不过这都是小事
关键是算法上，主要是两点
首先，数据不是独立的，n=k的情况和n=k+1有k个点是相同的
其次，样本容量为1的样本实在没什么说服力耶。我一般做实验都要多测几次算均值算方差算置信区间，虽然如果blog的一篇文章分析这么多我也会觉得很蛋疼，不过至少也多测几次算个均值吧orz
太复杂了wwwwwww
是不是和当初的fb及后来的vb一样道理？
不擅长理科和编程，完全看不懂。
这个GNUPLOT第一次见呢… 话说是基于X11运行的？
技术文啊。。算法啊。。啊啊啊啊啊
@conanshang:
不是，基于Terminal运行的，X11只是提供窗口而已
@ciasoMs:
蒙地卡罗法就不存在“样本容量”这个概念……嘛硬要说的话那个n值本身就是“样本容量”了吧
header文件偶是照抄老师发的单子上面的，反正偶没学过C++
8过算法的吐槽明显吐错了吧
@Phoenix:
i是目前为止生成随机函数的次数，也就是生成的坐标点的数量，而a是目前为止这些随机生成的坐标点里面在扇形范围内的坐标点的数量，因为可以把这些点看做面积相同的比如圆圈之类的有大小的图像，所以扇形内的点的数量除以整个方框内的点的数量大致就是这个扇形占这个方框的面积的几分之几哦 
然后因为整个方框的面积是1*1=1，所以这个几分之几实质上就是扇形的面积；扇形是1/4圆，所以乘以4就是圆形的面积
和C++以及几何无关
@Bill gates hxk:
你这句话直接把全世界的超级计算机的存在意义都给抹杀了
@Lovee:
这样啊，
Mac下Terminal的指令和Unix一样的么？偶对Unix还一窍不通呢…
@conanshang:
8过偶也没系统学过UNIX不敢完全肯定，只知道Leopard开始Mac OS X就通过了UNIX认证，也就是说任何针对UNIX编写的程序都不需要经过修改就能原生在Mac OS X上运行了
内核是一样的所以基本上命令应该绝大部分都是一样的吧
@Lovee:
嗯嗯，这样貌似方便不少呢，
上周有个HP来偶们学校找实习生，笔试题目里关于操作系统的都是Unix的…偶就爆了
话说，偶的博客的评论邮件回复正常么》上次无意被偶关掉后忘记开了…后来就没测试过
@conanshang:
回复通知很正常哦
@Lovee:
谢谢啦
@Lovee:
据我了解超级计算机在往图形化发展，比如nVidia的CUDA（官方称：个人超级计算机）
@Lovee:
有的，只要是测量就有样本容量
具体到这里，执行一次N点计算程序就是一次实验/一次测量，执行n次并记录就得到一个容量为n的样本
顺便UNIX指的是接口不是实现，内核模型不一样的话编译成二进制文件的软件肯定不能直接用的
针对UNIX编写的程序不经过修改就能在Mac上编译通过，但有的软件的编译就是件困难的事情
@ciasoMs:
前后两条评论连起来看了足足5分钟之后终于理解你要说的意思了orz
于是：
1、电脑本来就只能生成有限精度的小数，而且随机函数准确来说也只是“伪随机”
2、这里只是求随机生成点数的数量的增加对计算精度的影响的推移（或者说精度的变化趋势），而非准确地通过蒙地卡罗法计算圆周率
所以使用i+1个随机点计算出来的圆周率比使用i个随机点计算出来的有什么差别这一点完全不需要重置前面i个点——因为从根本上来说偶就不是要求证这个关系
事实上图上的两个箭头都只是随意画出来而非严格通过统计学的最小平方偏差法严谨地求出来的，所以不论偶是否重置前i个点最后生成的结果对于偶画这条线不会有任何影响
参考
8过这倒是让偶想到了另一个或许会比较有趣的实验——执行100次10000个随机点得出的结果的平均，和执行1000次1000个随机点得出的平均，会是哪一个更接近准确值呢
@Bill gates hxk:
图形化和GUI化是两个概念——GNUPLOT还是把计算结果图形化呢
@conanshang:
8用
@Lovee:
不一样哦，CUDA是直接抛弃了CPU使用的是GPU来进行图形运算
所以图形运算几乎是瞬时的不像一般的电脑还需要等
我要是吐槽用Super PI 来计算圆周率会不会有人砸我______WWwwwwWWWww
@Bill gates hxk:
于是你不知道这个世界上有一种东西叫做GPGPU还有一种东西叫做OpenCL么
另外偶已经说过了GUI化和图形化是两个概念
@Lovee:
好吧..我只是想说图形化操作是一种趋势而已..
@Bill gates hxk:
用GPU不代表就是GUI化，说到底类似于Earth Simulater这类的超级计算机的用途上只用Source Code做复杂科学处理不论是开发成本还是维护成本抑或者程序修改方便程度都比还要弄个GUI要好多了——不就是模拟个地球环境变化么还专门搞个GUI来干嘛用起来反而麻烦
所以偶说了图形化（准确来说是可视化）和GUI化是两个概念
@Lovee:
不扯那么远…有个GUI好过去输那些麻烦的命令行..
如果是重复大量调用的话还是用图形界面更好些..
说起来这个好像可以直接使用Number或者Excal之类的工具完成而不需要专门写个程序…
@Bill gates hxk:
刚好相反——大量的重复调用反而用命令更方便，尤其是针对那些“真正需要用到这种大型复杂科学计算”的人来说，比如用命令行和C写这个东西只需要最多不超过10分钟，加上整个思考程序架构的时间也不会超过30分钟；要用GUI写这个玩意儿光是码代码和窗口环境关联以及美化等就不止10分钟，还不算最开始的架构思考和调试，而且更重要的是要维护这玩意儿纯CLI的话只要看代码就足够了GUI的话还要看其他的各种东西
至于Numbers或者Excel——你试试看用Excel随机生成100000000个点球圆周率的误差然后把这些误差像上面的GNUPLOT一样画个图出来看看=_,=不说100000000个点，光是1000个点的误差要用Numbers生成一个折线图耗时就不止15分钟（反正偶是等了15分钟就受不了强制退出了），Excel虽然（在Windows上）的执行效率比Numbers要高点但是恐怕也不可能迈过1分钟这个坎——但是1000个点用GNUPLOT的话只需要大概1秒钟的时间=_,=这也是为什么有时候GUI反而不一定有CLI好用的因素之一
@Lovee:
GUI党表示看到代码很长就心烦…
特别是便携设备上如果木有图形界面，那绝对是反人类…
使用那些办公工具自然效率低，但是至少是可行的特别是在ipad等设备上…
@Bill gates hxk:
好吧你终于说到点子上了——于是你见过在iPad上模拟地球环境变化的么
顺便一说要在iPad上的Numbers描1000个点的折线图绝对中途Crash，管你用1还是2
虽然不是Geek但是不是所有程序都适合用GUI
@Lovee:
表示那些后台的东西交给云处理就好了…（我不是指iwork
mac有个好处就是可以用Xcode写GUI程序很方便的调用插件，当然也可以用auto那啥的创建个GUI类似image2pdf那样，多次操作的话还是有个GUI比较舒服，如果只是偶尔用用的话还是只用命令行好了
表示…………kita只能看得懂圖形
其实要用到随机数（乱数）的地方 特别是对随机数的随机性要求高的地方 我都不建议用c/c++来做
嘛 不过你们专业如此。。。
不过一様乱数已经算好了 如果再做变换函数 那误差更是伤不起啊
话说抛针100000000次才花掉半个小时 我要想要女仆本啊
@hicrokee:
（C++只是照着学校的单子抄的
小小表示只会一点C和BASIC
@ciasoMs:
1，这不是混用，本来就是继承过来的语法，不然我们为何黑c++
2，srand一次就够，因为只运行一次
3，这是情报工很基础的试验，完全没问题。不是每个蒙特卡罗都用来做时序列模拟的。
@hicrokee:
1. You’ll usually have problems if you try to intermix the two forms in a single program. –Thinking in C++ Page.87
当然我没碰到过问题所以到底有什么问题我也不知道
2.我的意思是srand使用了常数作参数，这样多次运行程序会有同样的结果，起不到随机的作用
@Lovee:
设第n次时圆弧内的点数是a[n]
如果第n+1次的结果是r（点在圆内r=1，不在圆内r=0），那总和就变成a[n+1]=a[n]+r，期望E(a[n+1])=E(a[n])+E(r)=a[n]+pi/4，此时算出的\hat{pi}=a[n+1]*4/(n+1)=a[n]*4/(n+1)+r/(n+1)。换句话说第n+1次实验的结果是依赖于前一次的
把这样算出来的误差放在一张图里画出来，说是蒙氏方法随着数据量增加的误差变化趋势，我觉得不妥
@ciasoMs:
1，这问题我也没遇到过，一般不用c++做实验，double的精度放在那里。
2，多次运行当然用时间作变量srand比较好。
3，我理解你的意思。这种做法下a(n)跟a(n+k)存在自相关，直接比较的确不妥。只是展示a(n)×4/n随n增加而收束于某值的话，这个实验的做法是可以接受的。至于提高重复试验求均值方差就233了，跟增大n没什么两样。
@hicrokee:
重复实验和增大n是两回事
重复实验是为了增加实验结果的可信度，增大n是为了研究准确度和n的关系，没有前者的话后者是不成立的
如果让我做的话应该是以下的步骤
1. 对于n=1的情况重复实验10次，记下实验结果
2. 用蒙氏法求pi的点估计（其实就是矩估计）和某置信度下的置信区间
3. n++，回到1，直到n足够大
4. 画出pi的置信区间长度对n的曲线
好像表示实验误差不应该用点估计值而应该用置信区间？不是很确定
ps.写上面的东西顺便复习了一下数理统计，发现有的东西全忘了
@ciasoMs:
@hicrokee:
结论：果然工科的试验思考方式和理科的试验思考方式不同么
p.s. 球更好的随机数语言
@Lovee:
其实我是工科生
买个硬件随机数发生器吧
@ciasoMs:
为了一个区区耗时两节课的实验去买个硬件随机数发生器除非学校啥时候那么大方给偶报销否则偶才不要去买
2020年前来观摩。
楼主的网站内容充实啊。
我是觉得蒙特卡罗方法精度太低，但为什么精度这么低呢？不够随机？
现在在用py，楼主毕设竟然是cuda相关的，13年竟然就开始cuda编程了。我用numba。
很有趣。