2020-03-08

水際で観察

河川水生昆虫生き物

川に行ったら水際を観察しましょう。

水際は面白い環境で、仔稚魚の生育場になったり、特定の底生生物の生息場になったりします。目を凝らすと案外生き物が見つかります。また、植物や水生昆虫が羽化した抜け殻が漂着します。この漂着物を探すだけでも楽しいものです。

このように、陸からの観察でも水生生物の生活の一側面を知ることができます。水に入る準備がなくても楽しめるのが利点です。ふらっと川に行ってできるよい遊びになります。

陸から撮った写真をいくつか掲載します。

f:id:y_tis:20200308175825j:plain:w400

f:id:y_tis:20200308175921j:plain:w400

f:id:y_tis:20200308175941j:plain:w400

日没前のほんの一時間程度でしたが、良い時間を過ごせました。

2019-11-30

読書感想文 : なぜ、仕事が予定どおりに終わらないのか ? 「時間がない病」の特効薬 ! タスクシュート時間術

読書感想文

本稿は以下の書籍の読書感想文*1です。

佐々木正悟・大橋悦夫 (2014) なぜ、仕事が予定どおりに終わらないのか ? 「時間がない病」の特効薬 ! タスクシュート時間術. 技術評論社 *2.

思ったように時間を作れないことがきっかけで読み始めた本です。作業ごとに時間を見積もり計画を立てるというのは、私も行ってきました。しかし、なかなか予定通りに進まないものです。何を改善すればよいのかを考える材料になればと思い手に取りました。

本書の指摘の内、特に重要なのは、1 分以上かかる作業は必ず時間を見積もることです。実際に作業を行うにあたり、ほとんどの場合小さな準備時間が入ります。この準備時間も見積もりに入れておかなければなりません。小さな休憩時間でも 1 分以上使うなら見積もりに入れておく必要があります。仕事はやった順にしか終わらないため、小さな時間も考慮しなければ、いつの間にか仕事が想定よりも遅れてしまうことになります。一日のあらゆる行動を記録し、それぞれどれだけ時間がかかっているかを図ることで、一日に収まる仕事量であるかを判断することが重要です。

読了後、私はこれまでの時間の使い方を振り返りました。その結果、今までは準備時間の見積もりが出来ていなかったと考えました。例えば、帰宅する前の食器の片づけなど、小さな作業の時間が積み重なり、見積もりとの大きな差になっているのではないかと考えました。実際に、異なる作業間に 5--10 分程度の時間を使うことを想定すると、見積もりより予定が遅れることが少なくなりました。さらに、これまでの仕事量の想定は、一日には収めるには無理のあるものであったことがわかりました。

注意しなければならないのは休憩時間と準備時間を区別することです。作業後に 30 分休む場合、作業後に 5 分の片づけがあったとしても、休憩を 25 分にしてはなりません。そのようにしてしまうと、休憩が減り、その後の仕事の能率が下がってしまいます。本書でも、休憩時間を予備時間とすることを禁じています。同様に、睡眠時間も考慮しなければ、一日が十分な長さであるかを図れません。

以上のように、一日の中の全ての 1 分以上の行動と休憩を記録して、無理のない仕事量はどれだけであるかを把握しなければ、正確に仕事時間を見積もることができません。正確に見積もれないため、仕事が時間通りに終わらないという事態になります。私は仕事以外にかかっている時間の見積もりに失敗していたために、一日に収まる仕事量を見誤っていたということです。

上記の事柄以外にも「先送りを繰り返して時間を失わない技術」など、時間の使い方についても本書は言及しています。ただし、いわゆる隙間時間の活用のような、小さい時間の活用は唱えられません。隙間時間にできることは限られている上、隙間時間はあまり存在しないためです。一日は活動と休息で構成され、一日を活動でいっぱいにすることはむしろ避けるべきことであるとされます。

なぜ仕事が予定通りに終わらないのか、という問いに対し、それは時間があると思い込んで仕事を詰めすぎているからだ、というのが本書の回答です。本書で明らかにされるのは、時間捻出のための技術というよりも、実際に使える時間と思い込みの差がどれだけあるのかをあぶりだす技術です。現在やりたいことができているかに関係なく、もっと時間を捻出できるはずだと思っている人こそ、一読の価値があると思います。

*1:書評よりも読書感想文の方がしっくりきたので、読書感想文タグを作りました

*2:

なぜ、仕事が予定どおりに終わらないのか? ~「時間ない病」の特効薬!タスクシュート時間術

作者: 佐々木正悟,大橋悦夫
出版社/メーカー: 技術評論社
発売日: 2014/04/09
メディア: 単行本（ソフトカバー）
この商品を含むブログ (6件) を見る

2019-09-22

mt19937ar の一部を Modern Fortran で書いてみた

Fortran 統計学

統計計算を行うのに何かと乱数が必要になります。Fortran には組み込み関数として random_number() がありますが、
コンパイラによって実装が異なります*1 *2。環境によって品質が変わると困るので、自分で実装したものを使うのが安全だと思います。私はメルセンヌツイスターを用いています。これまでは mt19937 をモジュールとして実装していましたが、この度 mt19937ar をクラスにしてみました。記事の表題に Modern Fortran とあるのはモジュールやクラスを使ったから、くらいの意味です*3。

以下に mt19937ar.c を移植したものを示します*4。都合により、mt19937ar.c にある init_by_array() および genrand_int31() は実装していません。また、関数の名称を変えています。

コメントはできるだけ付けました。しかしながら、私がメルセンヌツイスターの数学的な部分を理解していないので、あまり気の利いたコメントが書けませんでした。コメントの書式は Javadoc 風ですが、ドキュメント生成ツールを使うためではなく、変数の説明を書くのにちょうどよい書式だと思ったからです。

実装ミスをはじめ、間違い等ございましたら、ご指摘くださいませ。

! mt19937ar.c を Fortran へ移植したプログラムです。
! ただし、init_by_array() および genrand_int31() は実装していません。
! 
! [0.0, 1.0]、[0.0, 1.0)、(0.0, 1.0) からの一様乱数を生成します。
!
! 元のプログラムは以下のリンク先にあります。
! http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
!
!
! ----- 移植元ライセンス ここから----- 
!/* 
!   A C-program for MT19937, with initialization improved 2002/1/26.
!   Coded by Takuji Nishimura and Makoto Matsumoto.
!
!   Before using, initialize the state by using init_genrand(seed)  
!   or init_by_array(init_key, key_length).
!
!   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
!   All rights reserved.                          
!
!   Redistribution and use in source and binary forms, with or without
!   modification, are permitted provided that the following conditions
!   are met:
!
!     1. Redistributions of source code must retain the above copyright
!        notice, this list of conditions and the following disclaimer.
!
!     2. Redistributions in binary form must reproduce the above copyright
!        notice, this list of conditions and the following disclaimer in the
!        documentation and/or other materials provided with the distribution.
!
!     3. The names of its contributors may not be used to endorse or promote 
!        products derived from this software without specific prior written 
!        permission.
!
!   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
!   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
!   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
!   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
!   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
!   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
!   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
!   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
!   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!
!   Any feedback is very welcome.
!   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
!   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
!*/
! ----- 移植元ライセンス ここまで----- 
!
module Class_mt19937ar
    use, intrinsic :: iso_fortran_env
    implicit none

    private

    !  パラメータ
    integer(int32), parameter :: n = 624
    integer(int32), parameter :: m = 397
    ! 計算に使うために定義したパラメータです。
    integer(int32), parameter :: n_1 = n - 1
    integer(int32), parameter :: n_2 = n - 2
    integer(int32), parameter :: m_1 = m - 1
    integer(int32), parameter :: n_m_1 = n - m - 1
    integer(int32), parameter :: n_m = n - m
    integer(int32), parameter :: m_n = m - n

    ! 2147483647 は 0x7fffffff の符号付き 10 進表記です。
    integer(int32), parameter :: lowerMask = 2147483647

    ! 0x80000000 の符号付き 10 進表記です。
    integer(int32), parameter :: upperMask = -lowerMask - 1


    ! mt19937 による一様乱数生成器です。
    type, public :: mt19937ar
        integer(int32), private :: mt(0:n-1)
        integer(int32), private :: mti

        contains

        ! [0, 1] から 32 bit 精度の一様乱数を生成します。
        procedure, pass :: generateOnCC

        ! [0, 1) から 32 bit 精度の一様乱数を生成します。
        procedure, pass :: generateOnCO

        ! (0, 1) から 32 bit 精度の一様乱数を生成します。
        procedure, pass :: generateOnOO

        ! [0, 1) から 53 bit 精度の一様乱数を生成します。
        procedure, pass :: generateOnCO_53bit
    end type

 
    ! コンストラクタです。
    interface new_mt19937ar
      module procedure :: new_mt19937ar
    end interface


    ! public にする関数
    public new_mt19937ar
     

    contains


    ! コンストラクタです。シード値を指定して乱数生成器を初期化します。
    !
    ! @param inutSeed 乱数シード
    ! @return 乱数生成器
    function new_mt19937ar(inputSeed) result(newGenerator)
        ! 入力
        integer(int32), optional, intent(in) :: inputSeed

        ! 変数
        integer(int32) seed

        ! 戻り値
        type(mt19937ar) newGenerator


        ! シード値を決定します。
        if(present(inputSeed)) then
            ! 指定されている場合は、指定の値をシード値とします。
            seed = inputSeed
        else
            ! 指定されていない場合はデフォルト値 5489 を用います。
            seed = 5489
        end if

        ! -1 は 0xffffffffUL の符号あり 10 進表記です。
        newGenerator%mt(0) = iand(seed, -1)

        block 
            ! ループ カウンター
            integer(int32) i

            do i = 1, n_1
                newGenerator%mt(i) = 1812433253 &
                  & * (ieor(newGenerator%mt(i-1), ishft(newGenerator%mt(i-1), -30))) &
                  & + i

                ! -1 は 0xffffffffUL の符号あり 10 進表記です。
                newGenerator%mt(i) = iand(newGenerator%mt(i), -1);
            end do
            newGenerator%mti = n
        end block
    end function


    ! [0x0, 0xffffffff] の範囲から符号あり整数の乱数を返します。
    ! 
    ! @param this 乱数生成器
    ! @return [0x0, 0xffffffff] の範囲からの符号あり整数の乱数
    function generate_int32(this) result(y)
        ! 引数
        class(mt19937ar), intent(inout) :: this

        ! 戻り値
        integer(int32) y


        if (this%mti >= n) then
            call generateNWordsAtOneTime(this, y)
        end if

        y = this%mt(this%mti)
        this%mti = this%mti + 1

        ! Tempering
        ! ishft(i, j) で j =  n とすると左に n bit シフトします。
        ! ishft(i, j) で j = -n とすると右に n bit シフトします。
        y = ieor(y, ishft(y, -11))
        ! -1658038656 は 9d2c5680 の符号あり 10 進表記です。
        y = ieor(y, iand(ishft(y, 7), -1658038656))
        ! -272236544 は efc60000 の符号あり 10 進表記です。
        y = ieor(y, iand(ishft(y, 15), -272236544))
        y = ieor(y, ishft(y, -18))
    end function


    ! mti >= n の時に、内部状態を更新します。
    ! 
    ! @param this 乱数生成器
    ! @param y 符号あり整数
    subroutine generateNWordsAtOneTime(this, y)
        ! 引数
        class(mt19937ar), intent(inout) :: this
        integer(int32), intent(inout) :: y

        ! パラメータ
        ! -1727483681 は 9908b0df の符号なし 10 進表記です。
        integer(int32), parameter :: matrixA = -1727483681
        integer(int32), parameter :: mag01(0:1) = (/0, matrixA/)


        ! ループ カウンター
        integer(int32) kk

        do kk = 0, n_m_1
            y = ior(iand(this%mt(kk)  , upperMask), &
                  & iand(this%mt(kk+1), lowerMask))

            this%mt(kk) = ieor(ieor(this%mt(kk+m), ishft(y, -1)), &
                             & mag01(iand(y, 1)))
        end do

        do kk = n_m, n_2
            y = ior(iand(this%mt(kk)  , upperMask), &
                  & iand(this%mt(kk+1), lowerMask))

            this%mt(kk) = ieor(ieor(this%mt(kk+m_n), ishft(y, -1)), &
                             & mag01(iand(y, 1)))
        end do

        y = ior(iand(this%mt(n_1), upperMask), iand(this%mt(0), lowerMask))
        this%mt(n_1) = ieor(ieor(this%mt(m_1), ishft(y, -1)), &
                          & mag01(iand(y, 1)))

        this%mti = 0
    end subroutine 


    ! 以下、[0, 1]、[0, 1)、(0, 1) から一様乱数を発生させる関数です。
    ! 関数名の on の後につく O は 端点を含まないことを (open の意味) 、
    ! C は端点を含むことを (close の意味) 表します。
    !
    ! 例えば、
    ! generateOnCO --- [0, 1) からの一様乱数
    ! generateOnOO --- (0, 1) からの一様乱数
    ! を表します。


    ! [0, 1] から 32 bit 精度の一様乱数を発生させます。
    !
    ! @param this 乱数生成器
    ! @return [0, 1] からの一様乱数
    function generateOnCC(this) result(outputRandomNumber)
        ! 引数
        class(mt19937ar), intent(inout) :: this

        ! 変数
        integer(int32) intRandomNumber
        real(real64) convertedIntRandomNumber

        ! 戻り値
        real(real64) outputRandomNumber


        ! 符号ありの整数乱数を取得します。
        intRandomNumber = generate_int32(this)

        ! 符号なし整数としての値を、実数として取得します。
        convertedIntRandomNumber & 
          & = convertIntegerToUnsignedIntegerAsReal(intRandomNumber)

        ! [0, 1] の値に変換します。
        ! 4294967295 は 2^32-1 です。
        outputRandomNumber = convertedIntRandomNumber / 4294967295.0d0
    end function


    ! [0, 1) から 32 bit 精度の一様乱数を発生させます。
    !
    ! @param this 乱数生成器
    ! @return [0, 1) からの一様乱数
    function generateOnCO(this) result(outputRandomNumber)
        ! 引数
        class(mt19937ar), intent(inout) :: this

        ! 変数
        integer(int32) intRandomNumber
        real(real64) convertedIntRandomNumber

        ! 戻り値
        real(real64) outputRandomNumber


        ! 符号ありの整数乱数を取得します。
        intRandomNumber = generate_int32(this)

        ! 符号なし整数としての値を、実数として取得します。
        convertedIntRandomNumber & 
          & = convertIntegerToUnsignedIntegerAsReal(intRandomNumber)

        ! [0, 1) の値に変換します。
        ! 4294967295 は 2^32 です。
        outputRandomNumber = convertedIntRandomNumber / 4294967296.0d0
    end function


    ! [0, 1) から 53 bit 精度の一様乱数を発生させます。
    !
    ! @param this 乱数生成器
    ! @return [0, 1) からの一様乱数
    function generateOnCO_53bit(this) result(outputRandomNumber)
        ! 引数
        class(mt19937ar), intent(inout) :: this

        ! 変数
        integer(int32) a_int
        integer(int32) b_int
        real(real64) a_double
        real(real64) b_double

        ! 戻り値
        real(real64) outputRandomNumber


        ! 符号ありの整数乱数を取得します。
        a_int = ishft(generate_int32(this), -5)
        b_int = ishft(generate_int32(this), -6)
       
        ! 符号なし整数としての値を、実数として取得します。
        a_double = convertIntegerToUnsignedIntegerAsReal(a_int)
        b_double = convertIntegerToUnsignedIntegerAsReal(b_int)

        ! [0, 1) の値に変換します。 
        ! 4294967295 は 2^32 です。
        outputRandomNumber = (a_double * 67108864.0 + b_double) &
                             & / 9007199254740992.0d0
    end function


    ! (0, 1) から 32 bit 精度の一様乱数を発生させます。
    !
    ! @param this 乱数生成器
    ! @return (0, 1) からの一様乱数
    function generateOnOO(this) result(outputRandomNumber)
        ! 引数
        class(mt19937ar), intent(inout) :: this

        ! 変数
        integer(int32) intRandomNumber
        real(real64) convertedIntRandomNumber

        ! 戻り値
        real(real64) outputRandomNumber


        ! 符号ありの整数乱数を取得します。
        intRandomNumber = generate_int32(this)

        ! 符号なし整数としての値を、実数として取得します。
        convertedIntRandomNumber = &
          & convertIntegerToUnsignedIntegerAsReal(intRandomNumber)

        ! (0, 1) の値に変換します。 
        ! 4294967295 は 2^32 です。
        outputRandomNumber = (convertedIntRandomNumber + 0.5d0) &
                             & / 4294967296.0d0
    end function


    ! 符号あり整数を、符号なし整数を倍精度実数に変換した値にします。
    ! Fortran には unsigned long がないため、この関数で unsigned long
    ! としての値に変換します。
    ! 
    ! @param inputInteger 整数
    ! @return 符号あり整数を倍精度実数に変換した値
    function convertIntegerToUnsignedIntegerAsReal(inputInteger) result(outputReal)
        ! 引数
        integer(int32), intent(in) :: inputInteger

        ! 戻り値
        real(real64) outputReal


        if (inputInteger < 0) then
            ! 4294967296 は 2^32 です。
            outputReal = dble(inputInteger) + 4294967296.d0
        else
            outputReal = dble(inputInteger)
        end if 
    end function
end module


! メイン関数です。
! 使用例として、乱数生成器から 4 種類の乱数を生成します。
program MainProgram
    use, intrinsic :: iso_fortran_env
    use Class_mt19937ar
    implicit none

    ! 変数
    ! 乱数生成器
    type(mt19937ar) generator
    integer(int32) file_randomNumber

    ! ループ カウンター
    integer(int32) i


    ! 乱数生成器を生成します。
    generator = new_mt19937ar(5489)
    open(newunit=file_randomNumber, file='fortran_mt19937ar_CC.csv', status='replace')
    ! 乱数をファイルに出力します。
    do i = 1, 10000
        write(file_randomNumber, *) generator%generateOnCC() 
    end do
    close(file_randomNumber)

    ! 乱数生成器のシードを再設定します。
    generator = new_mt19937ar(5489)
    open(newunit=file_randomNumber, file='fortran_mt19937ar_CO.csv', status='replace')
    ! 乱数をファイルに出力します。
    do i = 1, 10000
        write(file_randomNumber, *) generator%generateOnCO() 
    end do
    close(file_randomNumber)

    ! 乱数生成器のシードを再設定します。
    generator = new_mt19937ar(5489)
    open(newunit=file_randomNumber, file='fortran_mt19937ar_OO.csv', status='replace')
    ! 乱数をファイルに出力します。
    do i = 1, 10000
        write(file_randomNumber, *) generator%generateOnOO() 
    end do
    close(file_randomNumber)

    ! 乱数生成器のシードを再設定します。
    generator= new_mt19937ar(5489)
    open(newunit=file_randomNumber, file='fortran_mt19937ar_CO_53bit.csv', status='replace')
    ! 乱数をファイルに出力します。
    do i = 1, 10000
        write(file_randomNumber, *) generator%generateOnCO_53bit() 
    end do
    close(file_randomNumber)
end program

*1:例えば、GCC 6.5 では KISS algorithm ですが、GCC 7.4 では xorshift1024* になっています。 gcc.gnu.org gcc.gnu.org

*2:様々な乱数を信頼できるアルゴリズムで実装した標準ライブラリを提供していないのは、科学技術計算用言語として致命的だと思います。

*3:Modern Fortran の厳密な定義は知りません。

*4:mt19937ar の説明は以下にあります。同リンク先の mt19937ar.c が移植元のプログラムです。 www.math.sci.hiroshima-u.ac.jp

2018-12-19

Inkscapeからemfを出力しepsに変換するとラスタライズされない

備忘録 Inkscape リンク集

Inkscapeでepsを出力するとラスタライズされてしまい，ベクター画像の意味がなくなります．対処法を探していたら以下のリンクに助けられました．
motchyの備忘録 InkscapeでEPSをエクスポートするときにラスタライズされてしまう時の対処

ラスタライズされる理由はよくわかりませんでした．なお，pdfで出力する場合はラスタライズされませんでした．

2018-12-04

分散の大きさの感覚を養う

Fortran 統計学

何事も実際にやってみるのは大切です．経験することで直観が養われます．

本記事は分散の大きさを感覚的に理解することが目的です．分散は平均とともに統計学で最初に習うものですから，定義を知っている人は多いでしょう，しかし，散布図を見て分散がどのくらいかを感覚的に答えられる人は少ないのではないでしょうか．

実験してみましょう．コンパイラはgfortran4.9.3です．

$x$ は $[0, 2\times3.1415)$ の一様乱数から発生させます．Fortranの組み込み手続きであるrandom_numberを用います．なお，random_numberで用いられている $[0, 1)$ からの一様乱数の発生アルゴリズムは，gfortranではxorshift1024*です*1．今回は $x$ は別に乱数でなくてもよいのですが，ランダムにとったデータぽくみせる演出です．

$y$ は平均 $\sin x$ ，分散 $\sigma^2$ の正規乱数から発生させます．正規乱数のアルゴリズムはZiggurat法です．実装方法は四辻（2010）*2を参照してください．

$\sigma^2= 1.0, 0.1, 0.01, 0.001$ の四つの場合で結果を比較しましょう．

$\sigma^2= 1.0$

$\sigma^2= 0.1$

$\sigma^2= 0.01$

$\sigma^2= 0.001$

正規分布なら各 $x$ において $\sin x \pm 2\sigma$ にその $x$ におけるデータの約95%が入ります．求めるものにもよりますが，分散1.0というのは結構大きく感じます．精密な較正曲線が求められる場面でこんなグラフを出したら怒られそうです．

最後に，今回用いたコードを以下に示します．

program sin_normal
implicit none
double precision x(100)
double precision y_1
double precision y_01
double precision y_001
double precision y_0001
integer n
integer i

n = size(x)

open(10, file="sin_normal.csv" , status="replace")
write(10, *) "x, y_1.0, y_0.1, y_0.01, y_0.001"
call random_number(x) 
x = 2.0d0 * 3.1415d0 * x
do i = 1, n
	y_1 = sin(x(i)) + random_normal(0.0d0, 1.0d0)
	y_01 = sin(x(i)) + random_normal(0.0d0, 1.0d-1)
	y_001 = sin(x(i)) + random_normal(0.0d0, 1.0d-2)
	y_0001 = sin(x(i)) + random_normal(0.0d0, 1.0d-3)
	write(10, *) x(i), ",", y_1, ",", y_01, ",", y_001, ",", y_0001
end do
close(10)

contains

function random_normal(mu, sigma2) result(rslt)
	double precision, intent(in)::mu
	double precision, intent(in)::sigma2 
	integer, parameter::k = 8
	integer, parameter::n = 2**k
	double precision, parameter::r = 3.6541528853610088d0
	double precision, parameter::v = 0.00492867323399d0
	double precision x(0:n)
	double precision y
	double precision u1
	double precision u_1
	double precision u__1
	double precision u___1
	double precision ux
	double precision u2
	double precision gu
	double precision gl
	integer sgn
	integer i
	integer j
	integer b
	integer l
	double precision rslt
	
	x(n) = v * exp(r * r * 0.5d0)
	x(n-1) = r
	do l = n-2, 1, -1
		x(l) = sqrt(-2.0d0 * log(exp(-x(l+1) * x(l+1) * 0.5d0) + v / x(l+1)))
	end do
	x(0) = 0.0d0
	
	do 
		!    ステップ 1
		i = 0
		j = 0
		b = 1
		call random_number(u1)
		
		do 
			!    ステップ2
			u_1 = u1 + u1
			if (u_1 < 1.0d0) then
				i = i + b
				u__1 = u_1
			else if(u_1 >= 1.0d0) then
				u__1 = u_1 - 1.0d0
			end if
			!    ステップ3
			j = j + 1
			if (j < k) then
				b = b + b
				u1 = u__1
				cycle ! ステップ2へ戻る
			end if
			exit !ステップ4へ進む
		end do

		!    ステップ4
		u___1 = u__1 + u__1
		if (u___1 < 1.0d0) then
			sgn = 1
			ux = u___1 * x(i+1)
		else
			sgn = -1
			ux = (u___1 - 1.0d0) * x(i+1)
		end if
		
		!    ステップ5
		if (ux < x(i)) then
			y = ux * dble(sgn)
			exit !ステップ7へ進む
		end if
		
		!    ステップ6-a
		if (i == n - 1) then
			y = random_normal_tail(r) * dble(sgn)
			exit !ステップ7へ進む
		!    ステップ6-b1
		else
			gu = exp(-(x(i) * x(i) - ux * ux) * 0.5d0)
			gl = exp(-(x(i+1) * x(i+1) - ux * ux) * 0.5d0)
			!   ステップ6-b2
			call random_number(u2)
			if (u2 * (gu - gl) <= 1.0d0 - gl) then
				y = ux * dble(sgn)
				exit !ステップ7へ進む
			end if
			cycle !ステップ1に戻る
		end if
	end do
	
	!    ステップ7
	rslt = mu + sqrt(sigma2) * y

end function

function random_normal_tail(b) result(rslt)
	double precision, intent(in)::b
	double precision d
	double precision x
	double precision u1
	double precision u2
	double precision rslt
	
	!    ステップ1
	d = b * b
	do
		!    ステップ2
		call random_number(u1)
		call random_number(u2)
		x = sqrt(d - (log(1.0d0 - u1) + log(1.0d0 - u1)))
		!    ステップ3
		if (x * u2 <= b) exit
		cycle
	end do
	rslt = x
end function

end program

*1:詳しいs仕様は以下のリンクを参照してください． gcc.gnu.org

*2:四辻哲章（2010）計算機シミュレーションのための確率分布乱数生成法．プレアデス出版．

計算機シミュレーションのための確率分布乱数生成法

作者: 四辻哲章
出版社/メーカー: プレアデス出版
発売日: 2010/06/01
メディア: 単行本
購入: 3人クリック: 27回
この商品を含むブログ (3件) を見る

2018-11-14

lectureとseminarの違い

備忘録英語語学

lectureとseminarの違いをOALD 8th edition*1で調べました．

lecture

a talk that is given to a group of people to teach them about a particular subject, often as part of a university or college course

seminar

a class at a university or college when a small group of students and a teacher discuss or study a particular topic

a meeting for discussion or training

lectureは誰かから与えられる話で，seminarの方は複数人で話し合うという感じですね．学生同士で勉強会をするときは，誰かが一方的に教えるだけならlecture，あれこれ議論を交えるならseminarという使い分けでしょうか．だれかが最初に一方的に話をしても，のちに議論が交わされるならseminarになるのだとおもいます．

lectureにはもう一つ意味があります．

lecture

a long angry talk that somebody gives to a one person or group of people because they have done something wrong

要はお説教ですね．何かを教える目的で話をするという点で一つ目の意味と共通しています．日本語だと講義と説教では全く字面が違いますが，英単語だとlectureとして一貫しているのが面白いです．単にangry talkではなくlongと付くのもいい味を出しています．

*1:

オックスフォード現代英英辞典第8版

作者: オックスフォード大学出版局
出版社/メーカー: 旺文社
発売日: 2010/09/11
メディア: 単行本
クリック: 5回
この商品を含むブログ (1件) を見る

Oxford Advanced Learner’s Dictionary, 8th edition (Oxford Advanced Learner's Dictionary)

作者: Hornby A. S.
出版社/メーカー: Oxford University Press
発売日: 2013/11/21
メディア: Kindle版
この商品を含むブログ (2件) を見る

2018-10-31

真の分布，統計モデル，事前分布はベイズ統計の基本要素

数学統計学備忘録

本稿ではベイズ統計の基本要素である真の分布，統計モデル，事前分布を説明し，真の分布のベイズ推測とは何かを説明します．参考書にしているのは渡辺（2012）*1，Watanabe（2018）*2および，渡辺澄夫先生の解説記事や講義資料です．

野外観測や実験等で我々はデータ $x^n = (x_1, \dots, x_n)$ をとります． $\forall i, x_i \in \mathbb{R}^N$ とします．これらのデータは真の分布（true distribution） $q(x)$ に従う確率変数 $X^n = (X_1, \dots, X_n)$ の実現値であると考えます．真の分布は一般的に知ることはできません．統計的推測においては，真の分布の予測分布 $p^*(x)$ （predictive distribution）を求めることが目的です．

そこで，我々は散布図やヒストグラムを書いてデータの傾向を見ながら，どんな確率分布ならデータを再現しそうかを考えます．その結果， $p(x|w)$ ならうまくいきそうだなと判断します．確率密度関数 $p(x|w)$ を統計モデルあるいは確率モデルといいます．統計モデルにはパラメータ $w \in W \subset \mathbb{R}^d$ が含まれています．統計モデルの例としては，

$p(x|\mu, \sigma) = \dfrac{1}{\sqrt{2\pi \sigma^2}} \exp \left(-\dfrac{(x - \mu)^2}{2\sigma^2} \right)$

などがあります．この例では $w = (\mu, \sigma)$ です．パラメータ $w$ が与えられた下での $x$ の分布なので，統計モデルは条件付き確率として定義されることに注意しましょう．

パラメータ $w$ と統計モデル $p(x|w)$ について，次の仮説を設定します．

パラメータがある確率分布 $\varphi(w)$ から発生し，データ $x^n$ がモデル $p(x|w)$ から独立に発生した

確率密度関数 $\varphi(w)$ を事前分布（prior distribution）といいます．

パラメータ $w$ とすべてのデータ $x^n = (x_1, \dots, x_n)$ について， $(w, x^n)$ の同時密度関数 $p(x^n, w)$ を考えると，

$\displaystyle p(w, x^n) = p(x^n|w)p(w) = p(x^n|w)\varphi(w) = \varphi(w)\prod_{i=1}^{n} p(x_i|w)$

となります．データ $x^n$ が与えられた時の $w$ の条件付き密度関数は，

$\displaystyle p(w| x^n) = \dfrac{p(w, x^n)}{p(x^n)} =\dfrac{1}{Z} \varphi(w)\prod_{i=1}^{n} p(x_i|w) \\\ Z = \int \varphi(w)\prod_{i=1}^{n} p(x_i|w)dw$