BMP に正規分布乱数を描く

正規分布乱数

Fortran の組み込み関数 random_number の乱数は区間 [0, 1) の一様分布ですが、これを用いてボックス=ミュラー法で正規分布に変換してみます。うまくゆけば平均 0、分散 1 の正規分布が出来上がります。そして、それを昨日作った BMP 用ルーチンでプロットしてみます（横軸がないのでヘロヘロですがｗ）。

ボックス=ミュラー法　wikipedia　

以下のように、二つの n 要素の一様乱数の入った配列から、2n 要素の正規分布乱数を求めています。

        integer, parameter :: n = 10**6
        real :: x(n), y(n), r(n), t(n), g(2 * n)
        call random_seed()
        call random_number(x)
        call random_number(y)
        r = sqrt(-2.0 * log(x))
        t = 2.0 * Pi * y
        g(    1:    n) = r * cos(t)
        g(n + 1:2 * n) = r * sin(t)

実行結果

二百万点について計算しています。平均、分散、グラフの見た感じは、予想通りになっている気がします。

 average= -2.2025492E-04  variance=  0.9992810
続行するには何かキーを押してください . . .

プログラム

Intel Fortran v.18, v.19beta で block 内で構造体を parameter としたときにうまく初期化が出来ないようです。文法的には、いいような気がするのですが？？

      module m_bmp
          implicit none

          type :: t_bmp_file_header
              sequence  
              integer(2) :: bfType = transfer('BM', 0_2, 1) ! BitMap
              integer(4) :: bfSize          ! file size in bytes
              integer(2) :: bfReserved1 = 0 ! always 0
              integer(2) :: bfReserved2 = 0 ! always 0
              integer(4) :: bfOffBits
          end type t_bmp_file_header


          type :: t_bmp_info_header
              sequence
              integer(4) :: biSize     = Z'28' ! size of bmp_info_header ; 40bytes 
              integer(4) :: biWidth
              integer(4) :: biHeight
              integer(2) :: biPlanes   = 1 ! always 1
              integer(2) :: biBitCount
              integer(4) :: biCompression = 0 !0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
              integer(4) :: biSizeImage
              integer(4) :: biXPelsPerMeter = 3780 ! 96dpi
              integer(4) :: biYPelsPerMeter = 3780 ! 96dpi 
              integer(4) :: biClrUsed      = 0
              integer(4) :: biClrImportant = 0 
          end type t_bmp_info_header


          type :: t_rgb
              sequence
              character :: b, g, r
          end type t_rgb 


          type :: t_bmp
              integer :: nx, ny
              type (t_rgb), allocatable :: rgb(:, :)
          contains 
              procedure :: wr => wr_bmp
              procedure :: rd => rd_bmp
          end type t_bmp 


          interface t_bmp
              module procedure :: init_bmp 
          end interface t_bmp

      contains   

          type (t_bmp) function init_bmp(nx, ny) result(bmp)
              integer, intent(in) :: nx, ny
              allocate(bmp%rgb(nx, ny))
          end function init_bmp


          subroutine wr_bmp(bmp, fn)
              class (t_bmp), intent(in) :: bmp
              character (len = *), intent(in) :: fn
              type (t_bmp_file_header) :: bmp_file_header
              type (t_bmp_info_header) :: bmp_info_header
              integer :: i
              associate(nx => size(bmp%rgb, 1), ny => size(bmp%rgb, 2))
                  bmp_file_header%bfSize      = 14 + 40 + 0 + (3 * nx + mod(nx, 4)) * ny
                  bmp_file_header%bfOffBits   = 14 + 40
                  bmp_info_header%biWidth     = nx
                  bmp_info_header%biHeight    = ny
                  bmp_info_header%biBitCount  = 24 
                  bmp_info_header%biSizeImage = (3 * nx + mod(nx, 4)) * ny
                  open(9, file = fn//'.bmp', access = 'stream', status = 'unknown')
                  write(9) bmp_file_header
                  write(9) bmp_info_header
                  write(9) (bmp%rgb(:, i), repeat(achar(0), mod(nx, 4)), i = 1, ny)
                  close(9)
              end associate
          end subroutine wr_bmp


          subroutine rd_bmp(bmp, fn)
              class (t_bmp), intent(out) :: bmp
              character (len = *), intent(in) :: fn
              type (t_bmp_file_header) :: bmp_file_header
              type (t_bmp_info_header) :: bmp_info_header
              integer :: i
              character :: dummy(4)
              associate(nx => bmp_info_header%biWidth, ny => bmp_info_header%biHeight)
                  open(10, file = fn//'.bmp', access = 'stream', status = 'old')
                  read(10) bmp_file_header
                  read(10) bmp_info_header
                  allocate(bmp%rgb(nx, ny))
                  read(10) (bmp%rgb(:, i), dummy(:mod(nx, 4)), i = 1, ny)
                  close(10)
              end associate  
          end subroutine rd_bmp
      end module m_bmp

      module m_pic
          use m_bmp
          implicit none

          type, extends(t_bmp) :: t_pic
          contains
              procedure :: point, line
          end type t_pic


          interface t_pic
              module procedure :: init_pic
          end interface t_pic

      contains    

         type (t_pic) function init_pic(nx, ny) result(bmp)
              integer, intent(in) :: nx, ny
              allocate(bmp%rgb(nx, ny), source = t_rgb(achar(0), achar(0), achar(0)))
         end function init_pic


         subroutine point(pic, ix, iy, rgb)
              class (t_pic), intent(in out) :: pic
              integer      , intent(in) :: ix, iy
              type (t_rgb) , intent(in) :: rgb
              pic%rgb(ix, iy) = rgb
         end subroutine point


         subroutine line(pic, ix0, iy0, ix1, iy1, rgb)
              class (t_pic), intent(in out) :: pic
              integer      , intent(in) :: ix0, iy0, ix1, iy1
              type (t_rgb) , intent(in) :: rgb
              integer :: ix, iy, mx, my
              real :: dx, dy, x, y
              mx = ix1 - ix0
              my = iy1 - iy0
              if (mx == 0 .and. my == 0) return
              if (abs(mx) > abs(my)) then
                  dy = my / real(mx)
                  y  = iy0 
                  do ix = 0, mx, sign(1, mx)
                      y = dy * ix
                      iy = nint(y)
                      call pic%point(ix0 + ix, iy0 + iy, rgb)
                  end do    
              else   
                  dx = mx / real(my)
                  x  = ix0 
                  do iy = 0, my, sign(1, my)
                      x = dx * iy
                      ix = nint(x)
                      call pic%point(ix0 + ix, iy0 + iy, rgb)
                  end do    
              end if    
         end subroutine line
    end module m_pic


    program Gaussian
        implicit none
        real, parameter :: Pi = 4.0 * atan(1.0)
        integer, parameter :: n = 10**6
        real :: x(n), y(n), r(n), t(n), g(2 * n)
        call random_seed()
        call random_number(x)
        call random_number(y)
        r = sqrt(-2.0 * log(x))
        t = 2.0 * Pi * y
        g(    1:    n) = r * cos(t)
        g(n + 1:2 * n) = r * sin(t)

        block 
            real :: av, sig2
            av   = sum(g) / (2 * n)
            sig2 = sum((g - av)**2) / (2 * n)
            print *, 'average=', av, ' variance=', sig2
        end block

        block
            use m_pic
            integer, parameter :: nx = 640, ny = 480
        !     type (t_pic), parameter :: White = t_rgb(achar(255), achar(255), achar(255)) ! ?syntax error?
            type (t_rgb) :: white 
            integer :: i, k, ih(nx) = 0
            type (t_pic) :: pic
            white = t_rgb(achar(255), achar(255), achar(255))
            do i = 1, 2 * n                                    ! histgram
                k = g(i) * 50 + nx / 2 
                ih(k) = ih(k) + 1      
            end do   
            pic = t_pic(nx, ny)
            call pic%line(     1, 1, nx    ,  1, white)        ! x-axis 
            call pic%line(nx / 2, 1, nx / 2, ny, white)        ! y-axis
            do i = 1, nx                                       ! plot Gaussian
                k = ny / 50 * ih(i) / int(sqrt(real(n))) + 1   
                call pic%point(i, k, white)
            end do   
            call pic%wr('gaussian')
        end block    
    end program Gaussian

Gaussian 二次元プロット　H30.6.26 追加

生成した二次元正規分布を点として書いて見ます。--水素原子 1s 軌道の確率分布相当になると思います。-- ごめんなさい、間違ってました。exp(-2x) の形でした。量子化学計算プログラム・ガウシアンでなら本当ですけどｗ

    program Gaussian
        implicit none
        real, parameter :: Pi = 4.0 * atan(1.0)
        integer, parameter :: n = 10**6
        real :: x(n), y(n), r(n), t(n), g(2 * n)
        call random_seed()
        call random_number(x)
        call random_number(y)
        r = sqrt(-2.0 * log(x))
        t = 2.0 * Pi * y
        g(    1:    n) = r * cos(t)
        g(n + 1:2 * n) = r * sin(t)

        block
            use m_pic
            integer, parameter :: nx = 1000, ny = 1000
        !     type (t_pic), parameter :: White = t_rgb(achar(255), achar(255), achar(255)) ! ?syntax error?
            integer :: i, ix, iy, ih(nx) = 0
            type (t_pic) :: pic
            type (t_rgb) :: white 
            white = t_rgb(achar(255), achar(255), achar(255))
            pic = t_pic(nx, ny)
            do i = 1, n 
                ix = 100000 * g(i    ) / sqrt(real(n)) + nx / 2
                iy = 100000 * g(i + n) / sqrt(real(n)) + ny / 2
                call pic%point(ix, iy, white)
            end do   
            call pic%wr('S')
        end block    
    end program Gaussian

（縮尺がまずいと升目状の線が見えますが、見かけのものです。）