编写stationary bootstrap 程序

假设数据集中有变量X，含有10个观测，所谓stationary bootstra的方法是首先从这10个任意取一个观测，比如说是第6个观测，那么下一个观测的取法是这样设计的，以概率p（事先是给定的，比如说是0.4）取除第6个观测以外的其他9个观测，而以概率1-p取第7个观测，也就是在第一个观测的观测号上递增一个，如此类推，第三个观测是在第二个观测号上按照上述概率来选取，直到取10次为止，请问如何使用SAS来完成上述要求的程序。

用R很容易写bootstrap的程序，SAS里用MACRO也能实现

data test;

array x{10};
do i=1 to 10;
x(i)=i*2;
end;

k=6;
y=x(k);
put  y= ;

do j=1 to 10;
p=0.4;
q=ranuni(1);
if q>p and k<10 then do;
  k+1;
  y=x(k) ;
  put  y= ;
  end;
else if q>p and k=10 then do;
  y=x(k) ;
  put  y= ;
  end;
else if q<=p then do;
  k=floor(ranuni(1)*9)+1;
  y=x(k);
  put y= ;
  end;
end;

run;

这里简单举出一个例子，更普遍的情况需要高手写成MACRO即可

以概率p（事先是给定的，比如说是0.4）取除第6个观测以外的其他9个观测

应该还是从10个观测里取还是从剩余的9个里取？

by my understanding:

1. data have;
   
2.         drop x1 j;
   
3.         N =10;
   
4.         p =.4;
   
5.         do sample =1 to 100;
   
6.                 x1 =.;
   
7.                 j =ranuni(1);
   
8.                 x =ceil(N*j);
   
9.                 output;
   
10.                 do _n_ =1 by 1 while(_n_ <=N-1);
    
11.                         x1 =x;  
    
12.                         j =ranuni(1);
    
13.                         if j <=p then do;
    
14.                                 x =ceil((N-1)*ranuni(1));
    
15.                                 if x =x1 then x =N;  
    
16.                         end;
    
17.                                 else do;
    
18.                                         x =x +1;  if x >N then x =1; *if x overflows then set as the first one;
    
19.                                 end;
    
20.                         output;
    
21.                 end;
    
22.         end;
    
23. run;
    
24. 
    
25. proc sgplot;
    
26. histogram x;
    
27. run;

复制代码

jingju11 发表于 2011-8-18 01:27
by my understanding:

我的理解是应该从其他9个中取，你的程序我有的地方看不懂，能否给出注释一下。样本容量假设为10，按照上述抽取方法抽取10个观测就完成一次抽样，如果把上述抽样进行比如100次，则如何编写呢？do _n_ =1 by 1 while(_n_ <=N-1);

                        x1 =x;  

                        j =ranuni(1);

                        if j <=p then do;

                                x =ceil((N-1)*ranuni(1));

                                if x =x1 then x =N;  

                        end;

这部分我是没有看懂，_n_没有终止值吗，  if x =x1 then x =N; 的作用是什么？能否在IML中实现这个程序，拜托。

学习了，呵呵~~~

我的做法是这样的：
假设N个记录，编号1-N.。预设概率p0.
1.从N个数里任意抽出一个x。输出
2.生成一个概率p ～uniform。
2a 如果p <=P0,那么从N-1（除掉上一次抽出的数）里任取一个；
2b 如果p >P0，那么取出x+1 的那个数。如果x+1>N,那么取数1.
输出
3.重复2，直到取出N个数。

jingju11 发表于 2011-8-18 09:57
我的做法是这样的：
假设N个记录，编号1-N.。预设概率p0.
1.从N个数里任意抽出一个x。输出

那这样的话取x+1的概率是1-p0+p0/(N-1)啦？不知道楼主的意思到底是什么。。如果要控制取x+1的概率为1-p0的话，那么应该是p<=p0的时候从除去x和x+1以外的数中等概率取一个。

cathy3212 发表于 2011-8-18 16:09
那这样的话取x+1的概率是1-p0+p0/(N-1)啦？不知道楼主的意思到底是什么。。如果要控制取x+1的概率为1-p0的 ...

我本身没有读过stationary boots的原文。
我的程序里，从第二个数开始，取x+1(X =LAGx) 的概率是1-p0.。取任一的非x的概率是P0/(N-1).
你的说法是对的：综合两个过程，取x+1(X =LAGx) 的概率是1-p0 +P0/(N-1).

京剧