******************************************** ;
*          ＳＡＳプログラム［１］            ;
******************************************** ;

* テストデータ生成（実際は観測データが存在） ;

data x ;
    do n = 1 to 100 ;
        r = ranuni( 0 ) ;
        select ;
            when( 1 <= n <= 35 ) x = 1 ;
            when(36 <= n <= 65 ) x = 2 ;
            when(66 <= n <= 93 ) x = 3 ;
            when(94 <= n <=100 ) x = 4 ;
        end ;
        output ;
    end ;
run ;
proc sort data = x out = x( keep = x ) ;
    by r ;
run ;


* ブートストラップ標本 2000 個作成する ;

data bootstrp ;
    do boot = 1 to 2000 ;
        do n = 1 to size ;
            obs = ceil( ranuni( 0 ) * 100 ) ;
            set x nobs = size point = obs ;
            output ;
        end ;
    end ;
    stop ;
run ;


* 標本比率の計算 ;

proc freq data = bootstrp  noprint ;
    by boot ;
    tables x / out = freq ;
run ;

data dif ;
    retain p1 p2 ;
    set freq ;
    by boot ;
    select( x ) ;
        when( 1 ) p1 = percent ;
        when( 2 ) p2 = percent ;
        otherwise ;
    end ;
    if last.boot then
        do ;
            dif = ( p1 - p2 ) / 100 ;
            output ;
        end ;
run ;


* 推定量 dif の平均と分散を計算する ;

title 'シミュレーション結果［１］' ;
proc univariate data= dif plot ;
    var dif ;
run ;
title ;


******************************************** ;
*          ＳＡＳプログラム［２］            ;
******************************************** ;

data bootstrp ;
    size = 100 ;
    do boot = 1 to 2000 ;
        p1 = 0 ;
        p2 = 0 ;
        do n = 1 to size ;
            x = ceil( ranuni( 0 ) * size ) ;
            select ;
                when( 1 <= x <= 35 ) p1 + 1 ;
                when(36 <= x <= 65 ) p2 + 1 ;
                otherwise ;
            end ;
        end ;
        dif = ( p1 - p2 ) / size ;
        output ;
    end ;
    stop ;
run ;

title 'シミュレーション結果［２］' ;
proc univariate data= bootstrp  plot ;
    var dif ;
run ;