按照特定要求生成随机分布

有如下变量（该变量为计数数据）：
var
2
0
3
12
5
8
.
.
.
20

问题一：
生成1000个新变量（随机分布），要求新变量的数值总和与原变量数值总和相等，并且新变量依旧为计数数据（非负整数）。

我根据要求写了如下代码：
sum var
return list
forvalues i=1(1)1000 {
gen var`i'=rpoisson(r(mean))
}

以上代码基本上可以解决问题，但是不完美。因为是计数数据，我因此使用了泊松分布，让新生成的泊松分布的均值与原变量保持一致。但问题是根据以上代码新生成的新变量并不能严格保证数值总和和原变量数值总和完全一致（尽管比较接近）。因此，希望能找到完美的解决方案。

问题二：
生成1000个新变量，新变量要求保持原变量的所有数值不变，只是改变观察值出现的顺序。

可用附件数据进行试验。

谢谢！


对于问题一，根据夏日贵志的思路，我改写代码如下（具体思路及出现的问题，详见下面的跟帖）：
*Simulation Distribution*
forvalues i=1/1000 {
use "test.dta", clear
gen order=_n
sum var
local sum=r(sum)
gen var`i'=rpoisson(r(mean))
sum var`i'
local diff=`sum'-r(sum)
gen random=rnormal() if var`i'!=0
replace random=-1 if var`i'==0
sort random
if `diff'>0 {
replace var`i'=var`i'+1 in f/`diff'
}
if `diff'<0 {
replace var`i'=var`i'-1 in `diff'/l
}
save "var`i'.dta", replace
}
*Merge*
use test.dta, clear
gen order=_n
forvalues i=1/1000 {
merge 1:1 order using "var`i'.dta", nogen
}
keep var*
save "simulation.dta", replace
*Erase*
forvalues i=1/1000 {
erase "var`i'.dta
}

对于问题二，根据夏日贵志的思路，我改写代码如下（夏日贵志的代码生成了一个新变量，我的问题要求生成1000个新变量，因此我在夏日贵志思路的基础上，简单嵌套了循环程序）：
forvalues i=1/1000 {
use "test.dta", clear
gen random`i'=rnormal()
sort random`i'
gen order=_n
drop random
rename var var`i'
save "var`i'.dta", replace
}
use test.dta, clear
gen order=_n
forvalues i=1/1000 {
merge 1:1 order using "var`i'.dta", nogen
}
keep var*
save "simulation.dta", replace
forvalues i=1/1000 {
erase "var`i'.dta"
}

以上代码的运算速度还是很快的，只是在merge的步骤相对费时一些（merge 1000次）。

第一个问题略麻烦一点，不过还算有意思。

1. * get sum of var
   
2. su var
   
3. local sum=r(mean)*r(N)
   
4. * gen var2 randomly - distribution does not matter
   
5. * uniform makes sure there is no negative number
   
6. gen var2=runiform()
   
7. * get the two sums to be the same
   
8. su var2
   
9. replace var2=var2/(r(mean)*r(N))*`sum'
   
10. * round var2
    
11. replace var2=round(var2)
    
12. * check by how much the two sums differ
    
13. su var2
    
14. local diff=`sum'-r(mean)*r(N)
    
15. * randomly choose a few observations of var2 to apply correction
    
16. gen var3=runiform() if var2>=1
    
17. sort var3
    
18. replace var2=var2+1 in 1/`diff'
    
19. drop var3
    
20. * verify the two means(i.e., sums) are exactly equal
    
21. su var var2

复制代码

第二个问题就很简单了
gen id=_n
preserve
gen i=rnormal()
sort i
replace id=_n
rename var var2
drop i
save tmp, replace
restore
merge 1:1 id using tmp, nogenerate
drop id
erase tmp.dta

非常感谢夏日贵志的参与，思路很巧妙。

夏目贵志 发表于 2017-3-11 23:32
第一个问题略麻烦一点，不过还算有意思。

非常感谢你的帖子，第一个问题目前确实能保证总和相等，但是两个变量的标准差之间的差距太大了。我一开始之所以选择泊松分布，因为计数数据大体服从泊松分布。但严格来讲这也有些问题。泊松分布满足均值等于方差，但实际上样本未必满足这个条件。另外，randomly choose a new observation of var2 to apply correction，其实并不是很随机的。生成的local diff确实可以巧妙地解决总和相等的问题，但是replace var2这步总改变从1到diff位置的计数值，而且就是在前几个位置上加1。这毫无疑问可以满足均值相等的问题，但是分配数值的方法并不随机。再有个小问题，我要生成1000个新变量，也许需要扩展到10000个，目前测试文件的观察值只有130，实际问题的观察值会比这个大得多，因此这个也要同时考虑代码运算速度和内存容量的问题。

我根据你的生成local diff的思路，依然使用泊松分布，改写代码如下：

*Simulation Distribution*
forvalues i=1/1000 {
use "test.dta", clear
gen order=_n
sum var
local sum=r(sum)
gen var`i'=rpoisson(r(mean))
sum var`i'
local diff=`sum'-r(sum)
sort var`i'
if `diff'>0 {
replace var`i'=var`i'+1 in f/`diff'
}
if `diff'<0 {
replace var`i'=var`i'-1 in `diff'/l
}
save "var`i'.dta", replace
}
*Merge*
use test.dta, clear
gen order=_n
forvalues i=1/1000 {
merge 1:1 order using "var`i'.dta", nogen
}
keep var*
save "simulation.dta", replace
*Erase*
forvalues i=1/1000 {
erase "var`i'.dta"
}

但是这个还是不能解决新变量标准差与原变量标准差之间差距较大的问题。我也想到了使用正态分布规定均值和标准差的方法，利用round取整，但是这会产生许多负值，可以replace负值为0，但是改动后标准差的差距又大了。另外，我觉得sort变量后首末位几个数据加减1的方法不够随机。这实际上是把最大值变小，把最小值变大。

非常感谢你的帖子，第一个问题目前确实能保证总和相等，但是两个变量的标准差之间的差距太大了。我一开始之所以选择泊松分布，因为计数数据大体服从泊松分布。但严格来讲这也有些问题。泊松分布满足均值等于方差，但实际上样本未必满足这个条件。

至于我为什么用了uniform，其实只是随手挑了一个。你要是用rnormal()^2也是一样的。并不是我觉得uniform一定对，别的一定不对。只是需要这些随机数不为负而已。要求样本的矩和总体完全一样本来就是不合理的。你就算从现在的样本里随机抽样，也没有任何办法保证的。所以建议你别纠结这一点了。

另外，randomly choose a new observation of var2 to apply correction，其实并不是很随机的。生成的local diff确实可以巧妙地解决总和相等的问题，但是replace var2这步总改变从1到diff位置的计数值，而且就是在前几个位置上加1。这毫无疑问可以满足均值相等的问题，但是分配数值的方法并不随机。

这个并没有不随机。
gen var3=runiform() if var2>=1
sort var3
其实就是在随机排序。

再有个小问题，我要生成1000个新变量，也许需要扩展到10000个，目前测试文件的观察值只有130，实际问题的观察值会比这个大得多，因此这个也要同时考虑代码运算速度和内存容量的问题。

我的这个方法并没有什么效率格外不高的地方。在我电脑上一千万观测值的情况下泡一下我这个程序也不到十秒钟。
r; t=8.08 10:15:44

但是这个还是不能解决新变量标准差与原变量标准差之间差距较大的问题。另外，我觉得sort变量后首末位几个数据加减1的方法不够随机。这实际上是把最大值变小，把最小值变大。

刚才解释过了。标准差之间有差距是正常的。你换一个和原分布相近的就能缓解这个问题。想要完全一致就算不是不可能，我个人觉得也没什么道理。sort的问题已经说过是随机sort的了。

另外，如果你回想一下自由度的概念，就能发现：你要保证均值和方差，就要牺牲自由度。换句话说就是不可能让你所有的观测值都随机。总有需要特别调整的。

夏目贵志 发表于 2017-3-12 22:18
至于我为什么用了uniform，其实只是随手挑了一个。你要是用rnormal()^2也是一样的。并不是我觉得uniform一 ...

对于标准差的问题，我同意你的观点，但是，目前的问题是标准差的差距过大，我希望差距可以减小一些。

对于随机分布的问题，你通过gen var3，再sort的方法，我理解了，确实可以保证随机。但是，有一点我没完全明白，我改写的代码中，加入了判断diff是正还是负的判断，如果diff大于0，则加1，如果diff小于0，则减1。你是怎么保证diff都是正值的呢？

另外，我是要生成新变量1000次，用你的代码，多测试几次，会出现如下错误：
. replace var2=var2+1 in 1/`diff'
'0' invalid observation number
r(198);
我改写的代码中由于加入了判断diff是否大于0的情况避免了这个问题。
所以你能不能显示一下嵌套了1000次循环的完整的代码。代码的执行效率受写循环的方法影响很大。我在上面代码中的循环速度还是可以的。但是，如果是10000个观察值，模拟1000次，再把所有文件合并成一个包含10001个变量的文件，我的代码是没办法在10秒内完成的，因为merge那步比较费时间。所以你说的10秒内完成，如果是10000个观察值，模拟1000次，再合并，那速度真的太快了。但是如果不包括所有以上步骤，那就不算快。

通过生成随机数，再按随机数排序保证随机加减1的方法，我改写我上面帖子中的代码如下：
*Simulation Distribution*
forvalues i=1/1000 {
use "test.dta", clear
gen order=_n
sum var
local sum=r(sum)
gen var`i'=rpoisson(r(mean))
sum var`i'
local diff=`sum'-r(sum)
gen random=rnormal() if var`i'!=0
replace random=-1 if var`i'==0
sort random
if `diff'>0 {
replace var`i'=var`i'+1 in f/`diff'
}
if `diff'<0 {
replace var`i'=var`i'-1 in `diff'/l
}
save "var`i'.dta", replace
}
*Merge*
use test.dta, clear
gen order=_n
forvalues i=1/1000 {
merge 1:1 order using "var`i'.dta", nogen
}
keep var*
save "simulation.dta", replace
*Erase*
forvalues i=1/1000 {
erase "var`i'.dta
}
这个程序如果用在10000个观察值的变量下，我用timer测试了时间，用时大概140秒。当然，机器配置不同，有的会快些，有的会慢些。

xingxf 发表于 2017-3-12 23:10
对于标准差的问题，我同意你的观点，但是，目前的问题是标准差的差距过大，我希望差距可以减小一些。

...

我只是提供了一个思路而已。这个和最后能做production use的代码不是一个概念。你指出的那个问题确实是对的。diff的那个地方检查一下符号也是很trivial的事情。

生成1000次什么的可以通过在1000 x _N的数据里生成一次来实现。不用一个一个的merge 1000次。

祝好运~

夏目贵志 发表于 2017-3-13 02:46
我只是提供了一个思路而已。这个和最后能做production use的代码不是一个概念。你指出的那个问题确实是对 ...

非常感谢你的思路。

至于你说在一个文件里一次生成1000次的问题能具体一点么？10000个观察值，1000个新变量，如果只是按我开头问题中呈现的方法，那还可以。但是如果考虑到你提供的思路，那中间的处理过程实际上要增加许多变量，不单个生成文件再merge的话，Stata处理起来恐怕更费时间。

高手过招！学习了！