ゼロからできるMCMC / Gaussian_Metropolis.c の Julia への移植

histogram.plt

# [reference]
# https://qiita.com/tell/items/5209b92d5f525ca7b028
# https://ja.wikipedia.org/wiki/ヒストグラム

reset session

set datafile separator comma

# DATA_FILE_PATH = '10e3.csv'
# DATA_FILE_PATH = '10e4.csv'
# DATA_FILE_PATH = '10e5.csv'
DATA_FILE_PATH = '10e6.csv'

stats DATA_FILE_PATH using ($1) nooutput

WIDTH_BIN = 2 * (STATS_up_quartile - STATS_lo_quartile) / ( STATS_records ** (1.0 / 3.0) )
BIN(x)    = WIDTH_BIN * floor( 0.5 + x / WIDTH_BIN )

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

set xrange[ - VAL_XRANGE : VAL_XRANGE ]
set boxwidth WIDTH_BIN

set logscale y 10
set format y "{10}^{%L}"

set key outside
set key right center
set key Left
set key reverse

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

plot \
    DATA_FILE_PATH \
    using ( BIN($1) ):( 1.0 / (WIDTH_BIN * STATS_records) ) \
    smooth freq \
    with boxes \
    title sprintf("N=%d", STATS_records) \
    , \
    TARGET_DIST(x)

main.jl

# [reference]
# ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門
# ISBN 978-4-06-520174-9
# https://github.com/masanorihanada/MCMC-Sample-Codes/blob/master/Gaussian_Metropolis.c

module McmcFromScratch

    using Random

    const FloatSize :: DataType = Float64

    function actionFunction( x::FloatSize )
        # return 0.5 * x * x
        return - log( exp( - 0.5 * (x - 3.0)^2 ) + exp( - 0.5 * (x + 3.0)^2 ) )
    end

    function executeMetropolisMethod( num_trials::Int64, step_size::FloatSize, file_path::String )

        # STEP.01
        # open a file to save result
        save_file = open(file_path, "w")

        # STEP.02
        # initialize the PRNG
        Random.Random.__init__

        # STEP.03
        # initialize
        # * the variable to store a generated sample
        # * the number of acceptance
        sample               = zero(FloatSize)
        num_accepted_samples = zero(num_trials)

        # STEP.04
        # main loop
        for iter in 1:num_trials

            # STEP.01
            # save the current sample
            sample_backup = sample
            action_init   = actionFunction(sample)

            # STEP.02
            # update the sample
            sample_variation  = rand()
            sample_variation  = 2 * step_size * (sample_variation - 0.5)
            sample           += sample_variation
            action_fin        = actionFunction(sample)

            # STEP.03
            # execute the Metropolis test
            if( exp(action_init - action_fin) > rand() )

                # when the new sample was accepted
                num_accepted_samples += 1

                # show the result
                println( save_file, sample, ", ", float(num_accepted_samples) / iter )

            else

                # when the new sample was rejected
                sample = sample_backup

            end

        end

        # STEP.05
        close(save_file)

    end

end

using ..McmcFromScratch

McmcFromScratch.executeMetropolisMethod( 10^3, 0.5, "10e3.csv" )
McmcFromScratch.executeMetropolisMethod( 10^4, 0.5, "10e4.csv" )
McmcFromScratch.executeMetropolisMethod( 10^5, 0.5, "10e5.csv" )
McmcFromScratch.executeMetropolisMethod( 10^6, 0.5, "10e6.csv" )