Fortran2008 による「ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門」の実装

README.md

README

Original

Gaussian_Metropolis.c

Used Environment

• GNU Fortran (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4
• gnuplot Version 5.4 patchlevel 3

main.f90

! [ORIGINAL]
! https://github.com/masanorihanada/MCMC-Sample-Codes/blob/master/Gaussian_Metropolis.c

! [REFERENCE]
! https://qiita.com/mt_west/items/ba9f674c4ee104f6b73d

program Gaussian_Metropolis

    ! required MODULE(s)
    use , intrinsic :: iso_fortran_env

    ! require all variables to be explicitly declared
    implicit none

    ! definition: user-defined TYPE
    type :: type_sample

        ! field(s) of this TYPE
        real    (REAL64) :: value
        logical          :: acceptance

    end type



    call execute_sampling( 10_INT32 ** 7 , 0.5_REAL64 , '10e7.dat' )



    contains



    pure elemental function action(x)

        ! dummy argument(s) for this FUNCTION
        real(REAL64) , intent(in) :: x

        ! return value of this FUNCTION
        real(REAL64) :: action

        action = 0.5_REAL64 * x * x

        ! eqn.(4.15) @ p.065 
        ! action = - log( exp( - 0.5_REAL64 * (x - 3.0_REAL64) ** 2 ) + exp( - 0.5_REAL64 * (x + 3.0_REAL64) ** 2 ) )

    end function



    function validity_Metropolis( sample_new, sample_old ) result(validity)

        ! dummy argument(s) for this FUNCTION
        real(REAL64) , intent(in) :: sample_new
        real(REAL64) , intent(in) :: sample_old

        ! declaration: variable(s) for this SUBROUTINE
        real(REAL64) :: criterion
            
        ! return value of this FUNCTION
        logical :: validity

        call random_number(criterion)
        validity = ( exp( action(sample_old) - action(sample_new) ) .gt. criterion )

    end function



    subroutine execute_Metropolis_method ( n_trials, step_size, sample )

        ! dummy arguments for this SUBROUTINE
        integer ( INT32       ) , intent( in    ) :: n_trials
        real    ( REAL64      ) , intent( in    ) :: step_size
        type    ( type_sample ) , intent( inout ) :: sample    (n_trials)

        ! declaration: variable(s) for this SUBROUTINE
        real(REAL64) :: x
        
        ! declaration: support variable(s) for this SUBROUTINE
        integer(INT32) :: iter

        ! STEP.01
        ! set the initial configuration
        x = 0.0_REAL64

        ! STEP.02
        ! generate samples using Metropolis method
        do iter = 1_INT32 , n_trials
        block

            ! declaration: variable(s) for this SUBROUTINE
            real(REAL64) :: x_backup
            real(REAL64) :: x_variant

            ! STEP.01
            ! get the variant of sample from a PRNG
            call random_number(x_variant)
            x_variant = 2.0_REAL64 * step_size * (x_variant - 0.5_REAL64)

            ! STEP.02
            ! update the sample
            x_backup = x
            x        = x + x_variant

            ! STEP.03
            ! save the generated samples
            sample(iter)%acceptance = validity_Metropolis( x, x_backup )
            sample(iter)%value      = x

            ! STEP.04
            ! update the configuration
            if ( .not. sample(iter)%acceptance ) then
                x = x_backup
            end if

        end block
        end do

    end subroutine



    subroutine execute_sampling ( n_trials, step_size, file_path )

        ! dummy arguments for this SUBROUTINE
        integer   ( INT32  ) , intent(in) :: n_trials
        real      ( REAL64 ) , intent(in) :: step_size
        character ( len= * ) , intent(in) :: file_path

        ! declaration: variable(s) for this SUBROUTINE
        integer :: save_unit

        ! declaration: variable(s) for this SUBROUTINE
        type(type_sample) , allocatable :: sample(:)

        ! STEP.01
        ! check the given arguments
        if ( n_trials .lt. 1_INT32 ) then
            stop 'The number of trials must be greater than ZERO.'
        end if

        if ( step_size .le. 0.0_REAL64 ) then
            stop 'The step size must be greater than ZERO.'
        end if

        ! STEP.02
        ! allocation
        allocate( sample(n_trials) )

        ! STEP.03
        ! open a file to save the result
        open( &!
            newunit = save_unit     , &!
            file    = file_path     , &!
            access  = 'stream'      , &!
            action  = 'write'       , &!
            form    = 'unformatted' , &!
            status  = 'unknown'       &! 
        )

        ! STEP.04
        ! execute the Metropolis method
        call execute_Metropolis_method( n_trials, step_size, sample(:) )

        ! STEP.05
        ! save the simulation result
        call export_samples( save_unit, sample(:) )

        ! STEP.06
        ! close the used file
        close( unit= save_unit, status= 'keep' )

    end subroutine



    subroutine export_samples ( save_unit, samples )

        ! dummy argument(s) for this SUBROUTINE
        integer               , intent(in) :: save_unit 
        type    (type_sample) , intent(in) :: samples   (:)
        
        ! declaration: variable(s) for this SUBROUTINE
        real(REAL64) :: denominator
        real(REAL64) :: mean_raw
        real(REAL64) :: mean_square

        ! declaration: support variable(s) for this SUBROUTINE
        integer(INT32) :: iter

        mean_raw    = 0
        mean_square = 0

        do iter = 1_INT32 , size( samples(:) )

            associate( sample => samples(iter) )

                if ( sample%acceptance ) then

                    denominator = 1.0_REAL64 / iter

                    mean_raw    = ( 1.0_REAL64 - denominator ) * mean_raw    + denominator * sample%value
                    mean_square = ( 1.0_REAL64 - denominator ) * mean_square + denominator * sample%value * sample%value

                    write( unit= save_unit ) &!
                        sample%value , &!
                        mean_raw     , &!
                        mean_square

                end if

            end associate

        end do

    end subroutine

end program

make_graph.plt

# [ORIGINAL]
# https://github.com/masanorihanada/MCMC-Sample-Codes/blob/master/Gaussian_Metropolis.c

# [REFERENCE]
# https://qiita.com/tell/items/5209b92d5f525ca7b028

reset session

# DATA_FILE_PATH = '10e5.dat'
DATA_FILE_PATH = '10e6.dat'

stats DATA_FILE_PATH binary endian=little format='%3double' using ($1) nooutput

VAL_XRANGE = ( abs(STATS_max) > abs(STATS_min) ) ? abs(STATS_max) : abs(STATS_min)
VAL_XRANGE = VAL_XRANGE * 1.2

WIDTH_BIN = (STATS_max - STATS_min) / ( ceil( log10(STATS_records) / log10(2.0) ) + 1 )
BIN(x)    = WIDTH_BIN * floor( 0.5 + x / WIDTH_BIN )

set boxwidth WIDTH_BIN
set format   y "{10}^{%L}"
set key      outside
set key      right center
set key      Left
set key      reverse
set xlabel   "accepted sample"
set ylabel   "relative frequency"
set logscale y 10
set xrange   [ - VAL_XRANGE : VAL_XRANGE ]
set yrange   [ 10 ** ( - ceil( log10(STATS_records) ) ) : 1.0 ]
set title    sprintf("Number of accpeted samples=%d", STATS_records)

TARGET_DIST(x) = exp(-0.5*x*x)/sqrt(2*pi)

# [REFERENCE]
# eqn.(4.15) @ p.065 
# TARGET_DIST(x) = ( exp( -0.5 * (x-3)**2 ) + exp( -0.5 * (x+3)**2 ) ) / ( 2 * sqrt(2*pi) )


plot \
    DATA_FILE_PATH \
    binary endian=little format='%3double' \
    using ( BIN($1) ):( 1.0 / (WIDTH_BIN * STATS_records) ) \
    smooth freq \
    with boxes \
    title "ACCEPTED SAMPLES" \
    , \
    TARGET_DIST(x)

pause -1

VAL_YRANGE = 10 ** ( ceil( log10(STATS_records) ) )

set   xlabel "number of accepted samples"
unset ylabel
set   format x '{10}^{%L}'
set   format y '%2.1f'
set   logscale x 10
unset logscale y
set   xrange [ 1.0 : VAL_YRANGE ]
unset yrange
set   xzeroaxis linetype -1 dashtype 2
set   arrow 1 from 1.0, 1.0 to VAL_YRANGE, 1.0 nohead linetype -1 dashtype 2

plot \
    DATA_FILE_PATH \
    binary endian=little format='%3double' \
    using 2 \
    with lines \
    title '<x>' \
    ,\
    DATA_FILE_PATH \
    binary endian=little format='%3double' \
    using 3 \
    with lines \
    title '<x^2>'

Makefile

# compiler selection
FC = gfortran

# compiler option: Common Options
FFLAGS_COMMON = -ffree-line-length-none -fimplicit-none -pedantic -std=f2008 -Wall -Werror -Wextra

# compiler option: Release mode
FFLAGS_RELEASE = ${FFLAGS_COMMON} -O3

# compiler option: Debug mode
FFLAGS_DEBUG = ${FFLAGS_COMMON} -O0 -s -fbacktrace -fbounds-check -g

# target of compilation
TARGET = ./Gaussian_Metropolis.exe

# object codes
OBJS = ./main.o

# suffix rule
.SUFFIXES:
	.o .f90

all: $(TARGET)

$(TARGET): $(OBJS)
	$(FC) -o $@ $(OBJS)

%.o: %.f90
	$(FC) $(FFLAGS) -c $<

clean:
	rm ./*.exe ./*.mod ./*.o ./*.smod

debug_mode:
	make clean; \
	make FFLAGS="$(FFLAGS_DEBUG)"

release_mode:
	make clean; \
	make FFLAGS="$(FFLAGS_RELEASE)"

run_exe:
	time $(TARGET)