一道非常有意思的stata编程模拟题~求模拟程序

柏拉图有一天问老师苏格拉底什么是爱情。
   
苏格拉底叫他到麦田走一次，要不回头地走，在途中要摘一颗最好的麦穗，但只可以摘一次。柏拉图充满信心地去了。谁知过了半天他仍没有回来，最后，他垂头丧气地出现在老师跟前，诉说空手而回的原因：“很难得看见一颗看似不错的，却不知是不是最好，不得已，因为只可以摘一次，只好放弃，再看看有没有更好的，到发现已经走到尽头时，才发觉手上一颗麦穗也没有。”这时，苏格拉底告诉他：“那就是爱情，爱情是一种理想，而且很容易错过。”


****计算：按照三七[exp(-1)]法则，选到最佳者的成功率
clear all
loc n=30 /*确定试验中的人数*/
loc f=10000 /*确定试验次数*/
loc s=0
set ob `n'
g r=.
loc e=int(_N*exp(-1)) /*预拒绝人数*/
forv i=1/`f'{
  replace r=uniform() /*各次试验中，随机确定各人禀赋及次序*/
  su r in 1/`e'
  loc me=r(max) /*被拒绝者中最优者*/
  forv j=`=`e'+1'/`n'{
    loc t=`j'
    if `me'<r[`j'] continue,break
  }
  loc r=r[`t'] /*所选者*/
  su r
  loc s=`s'+(`r'==r(max)) /*成功次数累计*/
}
n di "运用exp(-1)法则成功的频率：" 100*`s'/`f'  "%"



****找出成功频率最高的预拒绝人数（对于每种人数，试验10000次）
clear all
loc n=30 /*确定试验中的人数*/
loc f=10000 /*确定试验次数*/
set ob `n'
g r=.
forv k=1/`=`n'-1'{
  loc s=0
  forv i=1/`f'{
    replace r=uniform() /*各次试验中，随机确定各人禀赋及次序*/
    su r in 1/`k'
    loc me=r(max) /*被拒绝者中最优者*/
    forv j=`=`k'+1'/`n'{
      loc t=`j'
      if `me'<r[`j'] continue,break
    }
    loc r=r[`t'] /*所选者*/
    su r
    loc s=`s'+(`r'==r(max)) /*成功次数累计*/
  }
  loc bk=cond(`s'>`=`bs'+0',`k',`=`bk'+0')
  loc bs=cond(`s'>`=`bs'+0',`s',`=`bs'+0')
}
n di "成功频率最高的预拒绝人数："`bk'
n di "最高成功频率：" 100*`bs'/`f' "%"


****试验中需要注意的

分别使用原始的rnormal()与uniform()做试验，会得到不同的结果，这是因为数据精度不一致（从而`r'==r(max)结果失真）。

分别使用int(100*rnormal())与int(100*uniform())做试验，从而消除精度的影响，可以看到将得到基本一致的结果。

sungmoo 发表于 2012-8-22 15:03
柏拉图有一天问老师苏格拉底什么是爱情。
   
苏格拉底叫他到麦田走一次，要不回头地走，在途中要摘一颗最 ...

哈。在我们帖子中，哲学变成数学，理想成了优化。喜乎？悲乎？

voodoo 发表于 2012-8-22 14:12
发现模拟竟然和所用Stata版本有关。
我的是用Stata 12.1，你的应该是Stata 12.0。
假如do文档最开始加上一 ...

果然是厉害，连我用的是12.0你都能猜出来，牛！

sungmoo 发表于 2012-8-22 15:03
柏拉图有一天问老师苏格拉底什么是爱情。
   
苏格拉底叫他到麦田走一次，要不回头地走，在途中要摘一颗最 ...

强人，比我那个简单，也是牛人一枚，给好评，加分！

sungmoo大神在21楼直接实现了模拟过程。执行速度比本人8楼通过program sim ... -> simulate ...来得快（在我的电脑上该模拟过程仅需时1.7秒）。
现将sungmoo模拟过程嵌入一个循环中，让其实现寻找最优k的功能。

set more off
clear all
set seed 123456
loc n = 30 /*确定试验中的人数*/
loc f = 10000 /*确定试验次数*/

set ob `n'
g r = .
loc pk = 0

qui forv k = 1/`=`n'-1' {
        loc b = `k'+1
        loc s = 0
        forv i = 1/`f' {
                replace r = uniform() /*各次试验中，随机确定各人禀赋及次序*/
                su r in 1/`k'
                loc me = r(max) /*被拒绝者中最优者*/
                forv j = `b'/`n' {
                        if `me' < r[`j'] {
                                loc r = r[`j']
                                continue, break
                        }
                        loc r = r[`n']
                } /*所选者*/
                su r
                loc s = `s'+(`r'==r(max)) /*成功次数累计*/
        }
         
        if `=`s'/`f'' > `pk' {
                loc kmax = `k'
                loc pk = `s'/`f'
        }
}

di "取得最大概率的k = " `kmax' "；成功找到的概率 = " `pk' "。"

voodoo 发表于 2012-8-22 21:25
sungmoo大神在21楼直接实现了模拟过程。执行速度比本人8楼通过program sim ... -> simulate ...来得快（在我 ...

额，又进行优化了哈，厉害

ipony 发表于 2012-8-22 21:32
额，又进行优化了哈，厉害

厉害的是sungmoo大神！

voodoo 发表于 2012-8-22 21:34
厉害的是sungmoo大神！

额，他的程序比你的短些，不过运行也挺慢的，呵呵，你们都是牛人啊

ipony 发表于 2012-8-22 23:28  不过运行也挺慢的

这与如何设计命令没有太大关系。

设想总人数30，对于每种预拒绝人数都试验10000次，则总试验次数是290000。只要试验次数很大，不可能期望运行时间很短。

另外，个人以为，若想评估每种预拒绝人数的实际效果（从而找出最优人数），应该以以上为基础再做n次（倍）试验，从“平均意义”上看哪种预拒绝人数的“效果”最好。

有点意思，呵呵