i am BEST

IBM i のいろいろ、についてです。(主にデータベースとプログラミングについて)  http://www.iforum.ne.jp/index.php?topic=blogger13 から引越してきました。

2011年06月

ILE RPG で「乱数によるシミュレーション」(続き)

ILE RPG で「乱数によるシミュレーション」」の続きです。

前回は乱数の発生させ方のいろいろを紹介しましたが、今回は、前回収まりきらなかったポアソン乱数と、そうした乱数を利用した「大数の法則」と「中心極限定理」のシミュレーション例です。

Program 9-4 (ポアソン乱数)

      * Program 9-4 (ポアソン乱数)
     H DFTACTGRP(*NO) BNDDIR('QC2LE')                                                               
      *                                                                                             
     D pornd           PR            10i 0                                                          
     D  a                             8f                                                            
      *                                                                                             
     D exprnd          PR             8f                                                            
     D  a                             8f                                                            
      *                                                                                             
     D srand           PR                  extproc('srand')                                         
     D  var                          10u 0 value                                                    
      *                                                                                             
     D rand            PR            10i 0 extproc('rand')                                          
      *                                                                                             
     D log             PR             8f   extproc('log')                                           
     D  var                           8f   value                                                    
      *                                                                                             
     D time            PR            20i 0 extproc('time')                                          
     D  var                          20I 0                                                          
      *                                                                                             
     D a               S             10i 0                                                          
     D i               S             10i 0                                                          
     D numrnd          S             10i 0                                                          
     D seed1           S              8f   inz(1)                                                   
     D seed3           S              8f   inz(3)                                                   
     D NULL            S             20i 0 inz(x'00')                                               
      *                                                                                             
      /free                                                                                       
          srand(%uns(time(NULL))) ;                                                                 
                                                                                                    
          numrnd = 200 ;                                                                            
                                                                                                    
          for i = 1 to numrnd ;                                                                     
            a = pornd(seed3) ;                                                                      
            dsply %char(%dech(a:15:5)) ;                                                            
          endfor ;                                                                                  
                                                                                                    
          return ;                                                                               
      /end-free                                                                                     
      *                                                                                             
     P pornd           B                                                                            
     D pornd           PI            10i 0                                                          
     D  t                             8f                                                            
      *                                                                                             
     D t1              S              8f                                                            
     D k               S             10i 0 inz(0)                                                   
      *                                                                                             
      /free                                                                                       
          t1 = exprnd(seed1) ;                                                                      
                                                                                                    
          dow (t1<t) ;                                                                              
            k = k+1 ;                                                                               
            t1 = t1 + exprnd(seed1) ;                                                               
          enddo ;                                                                                   
                                                                                                    
          return k ;                                                                            
      /end-free                                                                                     
     P pornd           E                                                                            
      *                                                                                             
     P exprnd          B                                                                            
     D exprnd          PI             8f                                                            
     D  a                             8f                                                            
      *                                                                                             
     D b               S              8f                                                            
     D c               S              8f                                                            
      *                                                                                             
      /free                                                                                       
          b = rand()/32768.0 ;                                                                      
                                                                                                    
          if b = 0 ;                                                                                
            b = 1/32768.0 ;                                                                         
          endif ;                                                                                   
                                                                                                    
          c = -(1/a)*log(b) ;                                                                       
                                                                                                    
          return c ;                                                                            
      /end-free                                                                                     
     P exprnd          E                                                                            

Program 9-5 (大数の法則)

      * Program 9-5 (大数の法則)
     H DFTACTGRP(*NO) BNDDIR('QC2LE')                                                               
      *                                                                                             
     D mean1           PR             8f                                                            
     D  n                            10i 0                                                          
      *                                                                                             
     D rand            PR            10i 0 extproc('rand')                                          
      *                                                                                             
     D i               S             10i 0                                                          
     D n               S             10i 0                                                          
     D numrnd          S             10i 0                                                          
     D a               S              8f                                                            
      *                                                                                             
      /free                                                                                         
          numrnd = 200 ;                                                                            
          n = 5 ;                                                                                   
                                                                                                    
          for i = 1 to numrnd ;                                                                     
            a = mean1(n) ;                                                                          
            dsply %char(%dech(a:15:5)) ;                                                            
          endfor ;                                                                                  
                                                                                                    
          return ;                                                                                 
      /end-free                                                                                     
      *                                                                                             
     P mean1           B                                                                            
     D mean1           PI             8f                                                            
     D  n                            10i 0                                                          
      *                                                                                             
     D a               S              8f                                                            
     D b               S              8f                                                            
     D c               S              8f                                                            
     D i               S             10i 0                                                          
      *                                                                                             
      /free                                                                                         
          a = 0 ;                                                                                   
                                                                                                    
          for i = 1 to n ;                                                                          
            b = rand()/32768.0 ;                                                                    
            a = a+b ;                                                                               
          endfor ;                                                                                  
                                                                                                    
          c = a/n ;                                                                                 
                                                                                                    
          return c ;                                                                              
      /end-free                                                                                     
     P mean1           E                                                                            

Program 9-6 (中心極限定理)

      * Program 9-6 (中心極限定理)
     H DFTACTGRP(*NO) BNDDIR('QC2LE')                                                               
      *                                                                                             
     D mean1           PR             8f                                                            
     D  n                            10i 0                                                          
      *                                                                                             
     D mean2           PR             8f                                                            
     D  n                            10i 0                                                          
      *                                                                                             
     D rand            PR            10i 0 extproc('rand')                                          
      *                                                                                             
     D sqrt            PR             8f   extproc('sqrt')                                          
     D  var                           8f   value                                                    
      *                                                                                             
     D i               S             10i 0                                                          
     D n               S             10i 0                                                          
     D numrnd          S             10i 0                                                          
     D a               S              8f                                                            
      *                                                                                             
      /free                                                                                         
          numrnd = 200 ;                                                                            
          n = 5 ;                                                                                   
                                                                                                    
          for i = 1 to numrnd ;                                                                     
            a = mean2(n) ;                                                                          
            dsply %char(%dech(a:15:5)) ;                                                            
          endfor ;                                                                                  
                                                                                                    
          return ;                                                                                 
      /end-free                                                                                     
      *                                                                                             
     P mean1           B                                                                            
     D mean1           PI             8f                                                            
     D  n                            10i 0                                                          
      *                                                                                             
     D a               S              8f                                                            
     D b               S              8f                                                            
     D c               S              8f                                                            
     D i               S             10i 0                                                          
      *                                                                                             
      /free                                                                                         
          a = 0 ;                                                                                   
                                                                                                    
          for i = 1 to n ;                                                                          
            b = rand()/32768.0 ;                                                                    
            a = a+b ;                                                                               
          endfor ;                                                                                  
                                                                                                    
          c = a/n ;                                                                                 
                                                                                                    
          return c ;                                                                              
      /end-free                                                                                     
     P mean1           E                                                                            
      *                                                                                             
     P mean2           B                                                                            
     D mean2           PI             8f                                                            
     D  n                            10i 0                                                          
      *                                                                                             
     D a               S              8f                                                            
     D b               S              8f                                                            
      *                                                                                             
      /free                                                                                         
          a = mean1(n) ;                                                                            
          b = sqrt(12.0*n)*(a-0.5) ;                                                                
                                                                                                    
          return b ;                                                                                
      /end-free                                                                                     
     P mean2           E                                                                            

ILE RPG で「乱数によるシミュレーション」

Cによる統計データ解析入門」という本がたまたま手元にあったので、その中の「第09章 乱数によるシミュレーション」で紹介されているプログラム例を ILE RPG に書き換えてみました。(元のプログラムと同じような結果が出てくるのは確認済みです。何せ"乱数によるシミュレーション"なので毎回結果が違うので …)

こうした C の数学関数を使用した数値計算なども難なくできることがわかると思います。(元のプログラムは出力はファイルへの書き出しなのですが、今回は処理をシンプルに見やすくするために画面へのメッセージ出力にしています)

"地の文"なんかはほとんどそのままなことがわかりますね。

ぜひ、本か Web 上の回答例とくらべてみてください。

ILE RPG でできることっていうのは実はけっこう多いんですよ。

Program 9-1 (二項乱数)

      * Program 9-1 (二項乱数)
     H DFTACTGRP(*NO) BNDDIR('QC2LE')                                                               
      *                                                                                             
     D birnd           PR            10i 0                                                          
     D  n                            10i 0                                                          
     D  p                             4f                                                            
      *                                                                                             
     D bern            PR            10i 0                                                          
     D  p                             4f                                                            
      *                                                                                             
     D rand            PR            10i 0 extproc('rand')                                          
      *                                                                                             
     D p               S              4f                                                            
     D i               S             10i 0                                                          
     D k               S             10i 0                                                          
     D n               S             10i 0                                                          
     D numrnd          S             10i 0                                                          
      *                                                                                             
      /free                                                                                                         
          numrnd = 200 ;                                                                            
          n = 5 ;                                                                                   
          p = 0.5 ;                                                                                 
                                                                                                    
          for i = 1 to numrnd ;                                                                     
            k = birnd(n: p) ;                                                                       
            dsply %char(k) ;                                                                        
          endfor ;                                                                                  
                                                                                                    
          return ;                                                                                 
      /end-free                                                                                     
      *                                                                                             
     P birnd           B                                                                            
     D birnd           PI            10i 0                                                          
     D  n                            10i 0                                                          
     D  p                             4f                                                            
      *                                                                                             
     D i               S             10i 0                                                          
     D k               S             10i 0                                                          
      *                                                                                             
      /free                                                                                        
          k = 0 ;                                                                                   
                                                                                                    
          for i = 1 to n ;                                                                          
            k = k + bern(p) ;                                                                       
          endfor ;                                                                                  
                                                                                                    
          return k ;                                                                             
      /end-free                                                                                     
     P birnd           E                                                                            
      *                                                                                             
     P bern            B                                                                            
      *                                                                                             
     D bern            PI            10i 0                                                          
     D p                              4f                                                            
      *                                                                                             
     D j               S             10i 0                                                          
     D a               S              8f                                                            
      *                                                                                             
      /free                                                                                        
          a = rand()/32768.0 ;                                                                      
                                                                                                    
          if a < p ;                                                                                
            j = 1 ;                                                                                 
          else ;                                                                                    
            j = 0 ;                                                                                 
          endif ;                                                                                   
                                                                                                    
          return j ;                                                                              
      /end-free                                                                                     
     P bern            E                                                                            

Program 9-2 (正規乱数)

      * Program 9-2 (正規乱数)
     H DFTACTGRP(*NO) BNDDIR('QC2LE')                                                               
      *                                                                                             
     D normrnd         PR             8f                                                            
      *                                                                                             
     D norminv         PR             8f                                                            
     D  a                             8f                                                            
      *                                                                                             
     D rand            PR            10i 0 extproc('rand')                                          
      *                                                                                             
     D sqrt            PR             8f   extproc('sqrt')                                          
     D  var                           8f   value                                                    
      *                                                                                             
     D log             PR             8f   extproc('log')                                           
     D  var                           8f   value                                                    
      *                                                                                             
     D a               S              8f                                                            
     D i               S             10i 0                                                          
     D numrnd          S             10i 0                                                          
      *                                                                                             
      /free                                                                                       
          numrnd = 200 ;                                                                            
                                                                                                    
          for i = 1 to numrnd ;                                                                     
            a = normrnd() ;                                                                         
            dsply %char(%dech(a:15:5)) ;                                                            
          endfor ;                                                                                  
                                                                                                    
          return ;                                                                               
      /end-free                                                                                     
      *                                                                                             
     P normrnd         B                                                                            
     D normrnd         PI             8f                                                            
      *                                                                                             
     D a               S              8f                                                            
     D b               S              8f                                                            
      *                                                                                             
      /free                                                                                       
          a = rand()/32768.0 ;                                                                      
                                                                                                    
          if a = 0 ;                                                                                
            a = 1/32768.0 ;                                                                         
          endif ;                                                                                   
                                                                                                    
          b = norminv(a) ;                                                                          
                                                                                                    
          return b ;                                                                            
      /end-free                                                                                     
     P normrnd         E                                                                            
      *                                                                                             
     P norminv         B                                                                            
      *                                                                                             
     D norminv         PI             8f                                                            
     D  u                             8f                                                            
      *                                                                                             
     D a               S              8f                                                            
     D b               S              8f                                                            
     D w               S              8f                                                            
     D z               S              8f                                                            
      *                                                                                             
      /free                                                                                       
          z = -log(4*u*(1-u)) ;                                                                     
          a = z*(2.0611786-5.7262204/(z+11.640595)) ;                                               
          w = sqrt(a) ;                                                                             
                                                                                                    
          if u < 0.5 ;                                                                              
            b = -w ;                                                                                
          else ;                                                                                    
            b = w ;                                                                                 
          endif ;                                                                                   
                                                                                                    
          return b ;                                                                            
      /end-free                                                                                     
     P norminv         E                                                                            

Program 9-3 (指数乱数)

      * Program 9-3 (指数乱数)
     H DFTACTGRP(*NO) BNDDIR('QC2LE')                                                               
      *                                                                                             
     D exprnd          PR             8f                                                            
     D  a                             8f                                                            
      *                                                                                             
     D rand            PR            10i 0 extproc('rand')                                          
      *                                                                                             
     D log             PR             8f   extproc('log')                                           
     D  var                           8f   value                                                    
      *                                                                                             
     D seed            S              8f   inz(1)                                                   
     D a               S              8f                                                            
     D i               S             10i 0                                                          
     D numrnd          S             10i 0                                                          
      *                                                                                             
      /free                                                                                       
          numrnd = 200 ;                                                                            
                                                                                                    
          for i = 1 to numrnd ;                                                                     
            a = exprnd(seed) ;                                                                      
            dsply %char(%dech(a:15:5)) ;                                                            
          endfor ;                                                                                  
                                                                                                    
          return ;                                                                               
      /end-free                                                                                     
      *                                                                                             
     P exprnd          B                                                                            
     D exprnd          PI             8f                                                            
     D  a                             8f                                                            
      *                                                                                             
     D b               S              8f                                                            
     D c               S              8f                                                            
      *                                                                                             
      /free                                                                                       
          b = rand()/32768.0 ;                                                                      
                                                                                                    
          if b = 0 ;                                                                                
            b = 1/32768.0 ;                                                                         
          endif ;                                                                                   
                                                                                                    
          c = -(1/a)*log(b) ;                                                                       
                                                                                                    
          return c ;                                                                            
      /end-free                                                                                     
     P exprnd          E                                                                            

乱数のプログラミング例には、あとひとつポアソン乱数があるのですが、ここまででずいぶん長くなってしまったので次回に載せたいと思います。

記事検索
プロフィール

i_am_best

タグクラウド
QRコード
QRコード
  • ライブドアブログ