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

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

由于没人能知道真正的缘分何时到来，没人能知道下一个来求爱的男生会是什么样子，接受表白的时机早晚实在很难决定。怎么办？去向《非诚勿扰》的黄菡老师和乐嘉老师请教一下？其实你还可以向欧拉老师请教一下。你没听错。大数学家欧拉对一个神秘的数学常数 e ≈ 2.718 深有研究，这个数字和“拒人问题”竟然有着直接的联系。

“拒人问题”的数学模型

为了便于我们分析，让我们把生活中各种复杂纠纷的恋爱故事抽象成一个简单的数学过程。假设根据过去的经验，MM 可以确定出今后将会遇到的男生个数，比如说 15 个、30 个或者 50 个。不妨把男生的总人数设为 n。这 n 个男生将会以一个随机的顺序排着队依次前来表白。每次被表白后，MM 都只有两种选择：接受这个男生，结束这场“征婚游戏”，和他永远幸福地生活在一起；或者拒绝这个男生，继续考虑下一个表白者。我们不考虑 MM 脚踏两只船的情况，也不考虑和被拒男生破镜重圆的可能。最后，男人有好有坏，我们不妨假设 MM 心里会给男生们的优劣排出个名次来。

聪明的 MM 会想到一个好办法：先和前面几个男生玩玩，试试水深；大致摸清了男生们的底细后，再开始认真考虑，和第一个比之前所有人都要好的男生发展关系。从数学模型上说，就是先拒掉前面 k 个人，不管这些人有多好；然后从第 k+1 个人开始，一旦看到比之前所有人都要好的人，就毫不犹豫地选择他。不难看出，k 的取值很讲究，太小了达不到试的效果，太大了又会导致真正可选的余地不多了。这就变成了一个纯数学问题：在男生总数 n 已知的情况下，当 k 等于何值时，按上述策略选中最佳男生的概率最大？

如何求出最优的 k 值？

对于某个固定的 k，如果最适合的人出现在了第 i 个位置（k < i ≤ n），要想让他有幸正好被 MM 选中，就必须得满足前 i-1 个人中的最好的人在前 k 个人里，这有 k/(i-1) 的可能。考虑所有可能的 i，我们便得到了试探前 k 个男生之后能选中最佳男生的总概率 P(k)：
                              

用 x 来表示 k/n 的值，并且假设 n 充分大，则上述公式可以写成：


对 -x · ln x 求导，并令这个导数为 0，可以解出 x 的最优值，它就是欧拉研究的神秘常数的倒数—— 1/e ！

也就是说，如果你预计求爱者有 n 个人，你应该先拒绝掉前n/e 个人，静候下一个比这些人都好的人。假设你一共会遇到大概 30 个求爱者，就应该拒绝掉前 30/e ≈ 30/2.718 ≈11 个求爱者，然后从第 12 个求爱者开始，一旦发现比前面11 个求爱者都好的人，就果断接受他。由于 1/e 大约等于37%，因此这条爱情大法也叫做 37% 法则。

不过，37% 法则有一个小问题：如果最佳人选本来就在这37% 的人里面，错过这 37% 的人之后，她就再也碰不上更好的了。但在游戏过程中，她并不知道最佳人选已经被拒，因此她会一直痴痴地等待。也就是说，MM 将会有 37% 的概率“失败退场”，或者以被迫选择最后一名求爱者的结局而告终。

37% 法则“实测”！

37% 法则的效果究竟如何呢？我们在计算机上编写程序模拟了当 n = 30 时利用 37% 法则进行选择的过程（如果 MM 始终未接受求爱者，则自动选择最后一名求爱者）。编号越小的男生越次，编号为 30 的男生则表示最佳选择。程序运行 10000 次之后，竟然有大约 4000 次选中最佳男生，可见 37% 法则确实有效啊。

计算机模拟 10000 次后得到的结果


这个问题由数学家 Merrill M. Flood 在1949 首次提出，这个问题被他取名为“未婚妻问题”。这个问题的精妙之处在于，在微积分界叱咤风云的自然底数 e，竟也出人意料地出现在了这个看似与它毫不相关的问题中。不知道此问题在果壳网上发表后，Geek 男女间会不会多了一种分手的理由：不好意思，你是那 37% 的人……

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏4 回帖

关键词：stata编程 Stata 模拟程序 tata 有意思 编程 程序

在线等高人把下面的程序给模拟出来

感觉利用一个循环语句里面再插一个if语句使得第一个公式取最大值就可以了，但是就是编不出来！急死人啊

ipony 发表于 2012-8-21 21:34
感觉利用一个循环语句里面再插一个if语句使得第一个公式取最大值就可以了，但是就是编不出来！急死人啊

若可直接求解析解，为什么还要模拟？

set more off
set seed 123456

capture program drop sim
program sim, rclass
        args nGG k

        clear
        set obs `nGG'
        gen GG = rnormal()
        
        su GG in 1/`k', meanonly        // 前k个男生中“最优男生”取值
        gen dBetter = (GG > r(max)) // k之后是否有“更优男生”
        gen sBetter = sum(dBetter)        // k之后“更优男生”的排序
        
        su GG, meanonly
        local maxGG = r(max)                // “最优男生”取值
        su GG if sBetter == 1, meanonly        // k之后第一位“更优男生”取值

        return scalar dSUCCESS = (r(max) == `maxGG') // 成功找到“最优男生”
end

local nGG = 30  // 潜在男生数
local pk = 0
local kmax = 0
forval k = 1/`=`nGG'-1' {
        qui simulate dSUCCESS = r(dSUCCESS), reps(10000) nodots: sim `nGG' `k'
        qui su dSUCCESS, meanonly
        if r(mean) > `pk' { // 模拟中成功找到的概率
                local pk = r(mean)
                local kmax = `k'
        }
}

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

不知我模拟程序中存在什么问题（我没有太细看题目，）。
我计算出来的k倒是等于11，但依此方法在10000次随机模拟中找到最优男生的概率只有12.7%。

的Stata实现：
clear
local nGG 30
// P(k)公式计算如下：
set obs `=`nGG'-1'
gen k = _n
gen s = 1/k
gsort -k
gen sums = sum(s)
gen pk = k/`nGG' * sums

// 列出最大pk及其所对应的k
sort pk
list k pk in L

7楼的解析结果确实和5-6楼的模拟结果不同。也不知哪里出现问题，哪位坛友帮忙查查？

voodoo 发表于 2012-8-22 11:57
7楼的解析结果确实和5-6楼的模拟结果不同。也不知哪里出现问题，哪位坛友帮忙查查？

自问自答啦。发现“由于有语句return scalar dSUCCESS = (r(max) == `maxGG') ，应该是数值存储精度的问题”。修改如下：

set more off
set seed 123456

capture program drop sim
program sim, rclass
        args nGG k

        clear
        set obs `nGG'
        gen double GG = rnormal()
        
        su GG in 1/`k', meanonly        // 前k个男生中“最优男生”取值
        gen int dBetter = (GG > r(max)) // k之后是否有“更优男生”
        gen int sBetter = sum(dBetter)  // k之后“更优男生”的排序
        
        su GG, meanonly
        scalar maxGG = r(max)                // “最优男生”取值
        su GG if sBetter == 1, meanonly      // k之后第一位“更优男生”取值

        return scalar dSUCCESS = (r(max) == maxGG) // 成功找到“最优男生”
end

local nGG = 30  // 潜在男生数

local pk = 0
local kmax = 0
forval k = 1/`=`nGG'-1' {
        qui simulate dSUCCESS = r(dSUCCESS), reps(10000) nodots: sim `nGG' `k'
        qui su dSUCCESS, meanonly
        if r(mean) > `pk' { // 模拟中成功找到的概率
                local pk = r(mean)
                local kmax = `k'
        }
}

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

voodoo 发表于 2012-8-22 11:53
的Stata实现：
clear
local nGG 30

好厉害，成功把变量i给用k替代掉了，我之前一直觉得怎么这么多变量，就是写不出来这个求和函数！厉害！