裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

敵の敵は味方っていうじゃないですか...

2021年01月07日 | ブログラミング

まあ,はっきりいって,

R の敵は Python か?

じゃあ,Python の敵が Julia だったら,

Julia は R の味方じゃないですか...ニッコリ

Python なら できるんだけど,R には苦手,R じゃできないことって,どんなことがありますか?

Julia だったらできることなんじゃないですか?しかも,Python より上手に。

そんなことがあったら,知らせて下さい。

コメント

Julia で フィッシャーの正確確率検定

2021年01月07日 | ブログラミング

以前,Python でフィッシャーの正確確率検定を行うという記事を書いた。

今回は,もう Python は捨てようと,Julia でフィッシャーの正確確率検定を行う方法について記事を書く。

(1) http://netlib.org/toms/ から,643.gz をダウンロードする。

(2) 643.gz を展開してできる 643 というファイルを開き,1965 行目あたりからあるデータ行を削除する。

(3) 5行目(process the data file と書かれている行)の次に,以下のプログラムを挿入する。

      function fisher(nrow, ncol, table)
      implicit none

      integer    i, j, ncol, nrow
      double precision emin, expect, percnt, pre, prt, table(*)
      intent(in) nrow, ncol, table
      double precision fisher
      external   fexact
      expect = 0.0
      percnt = 0.0
      emin = 0.0
      call fexact (nrow, ncol, table, nrow, expect, percnt, emin, prt,
     &             pre)
      fisher = pre
      return
      end

(4) コマンドラインに以下のコマンドを入力する。

$ gfortran fisher.f -o libfisher.so -shared -fPIC

(5) エラーなくコンパイルできたら,コマンドラインに以下のコマンドを入力し,_fisher_ を含む行があるか確認する。

$ nm libfisher.so

     :
000000000000665f T _fisher_
     :

(6) 以下のプログラムを,ファイル名 fisher.jl でセーブする。

function fisher(tbl)
    nrow, ncol = size(tbl)
    ccall(("fisher_", "libfisher.so"),
        Float64,
        (Ref{Int32}, Ref{Int32}, Ref{Float64}),
        Int32(nrow), Int32(ncol), Array{Float64, 2}(tbl))
end

(7) Julia プログラムで,以下のようにして利用する。

include("fisher.jl") # fisher.jl のパス名に注意
p = fisher([1 3 6; 4 2 3])
println("p value = $p")

(8) julia の REPL だと

$ julia

  : 略

julia> include("fisher.jl")
fisher (generic function with 1 method)

julia> p = fisher([1 3 6; 4 2 3])
0.36348481240121944

julia> println("p value = $p")
p value = 0.36348481240121944

R で確認するなら

julia> using RCall

julia> R"fisher.test(matrix(c(1, 3, 6, 4, 2, 3), byrow=TRUE, ncol=3))"
RObject{VecSxp}

 Fisher's Exact Test for Count Data

data:  matrix(c(1, 3, 6, 4, 2, 3), byrow = TRUE, ncol = 3)
p-value = 0.3635
alternative hypothesis: two.sided

julia> R"options(digits=15)"
RObject{VecSxp}
$digits
[1] 7

julia> R"fisher.test(matrix(c(1, 3, 6, 4, 2, 3), byrow=TRUE, ncol=3))$p.value"
RObject{RealSxp}
[1] 0.363484812401219

コメント

Julia から FORTRAN の関数を呼ぶ

2021年01月07日 | ブログラミング

ネット上の情報は,Julia のバージョンが低かったり,新しい書き方の FORTRAN プログラムの場合(?)だったりして,そのまま使えなかった。

試行錯誤の結果,やっとまともに動く例を示すことができた。

まずは,古典的な書式で書かれた FORTRAN プログラム(関数副プログラム,サブルーチン副プログラム)

ファイル名 test.f

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     関数副プログラム

      function test1(nrow, ncol, table) ! 関数名が戻り値を表す

      implicit none
      integer i, j
      integer nrow, ncol           ! 配列のサイズ情報は引数で渡す
      double precision table(*)    ! 配列の受けとり方法は使い方次第(一次元配列の場合)
      intent(in) nrow, ncol, table ! 入力引数
      double precision test1       ! 戻り値(関数副プログラム名)の型の宣言

      do 20 i=1, nrow
          write (*, 99997) (table(i+(j-1)*nrow), j=1,ncol) ! 一次元配列として使用
   20 continue
99997 format (1x, 10f7.1)

      test1 = 0.0                  ! 戻り値の定義
      do 10 i = 1, nrow * ncol
          test1 = test1 + table(i)
   10 continue
      test1 = test1 / (nrow * ncol)
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     サブルーチン副プログラム

      subroutine test2(nrow, ncol, table, mean, cells) ! 戻り値は引数に指定

      implicit none
      integer i, j
      integer nrow, ncol              ! 配列のサイズ情報は引数で渡す
      double precision table(nrow, *) ! 配列のサイズは最終次元のみ * でよい
      double precision mean           ! 戻り値の宣言
      integer cells                   ! 戻り値の宣言
      intent(in) nrow, ncol, table    ! 入力引数
      intent(out) mean, cells         ! 戻り値の引数

      mean  = 0.0                     ! 戻り値の設定
      cells = nrow * ncol             ! 戻り値の設定

      do 20 i = 1, nrow               ! 配列要素の並び順は Julia と同じ(列優先)
          write (*, 99997) (table(i, j), j = 1, ncol) ! 二次元配列として使用
99997     format (1x, 10f7.1)
          do 10 j = 1, ncol
              mean = mean + table(i, j)
   10     continue
   20 continue
      mean = mean / cells
      end

次に,test.f をコンパイルし,動的ライブラリ(libtest.so)を生成する。

$ gfortran test.f -o libtest.so -shared -fPIC

できあがった動的ライブラリの中身を確認する。

$ nm libtest.so

出力される中に

0000000000000af4 T _test1_
0000000000000cd5 T _test2_

のような行がある。

_test1__test2_ を Julia プログラムの ccall 関数の第 1 引数として,("test1_", "libtest.so") ("test2_", "libtest.so")  として使う。

できあがった動的ライブラリを以下のような Julia プログラム test.jl から呼び出す。

# FORTRAN 関数副プログラムの例
function FuncProg(tbl)
    n, m = size(tbl)   # 配列サイズは引数で渡さなければならない
    nrow = Int32(n)    # 入力引数 Int64 を Int32 に変換
    ncol = Int32(m)    # 入力引数 Int64 を Int32 に変換
    table = Array{Float64}(tbl) # 入力引数 Array{Int64,2}, を Array{Float64,2} に変換
    # 入力引数 Int32, Int32, Array{Float64,2}
    # 戻り値  Float64 FORTRAN での double precision あるいは real*8
    mean = ccall(("test1_", "libtest.so"), # nm コマンドでは '_test1_' となっている
        Float64, # FORTRAN 関数プログラムからの戻り値の型
        (Ref{Int32}, Ref{Int32}, Ref{Float64}), # 引数の型
        nrow, ncol, table) # 実引数
    println("mean = $mean")
end

# FORTRAN サブルーチン副プログラムの例
function SubProg(tbl)
    n, m = size(tbl)   # 配列サイズは引数で渡さなければならない
    nrow = Int32(n)    # 入力引数 Int64 を Int32 に変換
    ncol = Int32(m)    # 入力引数 Int64 を Int32 に変換
    table = Array{Float64}(tbl) # 入力引数 Array{Int64,2}, を Array{Float64,2} に変換
    mean = [0.0]       # 戻り値,スカラーのときでも必ずベクトル(値は何でもよい)
    cells = [Int32(0)] # 戻り値は Int32 なので,受け取る変数も Int32 で定義
    # 入力引数: Int32, Int32, Array{Float64,2} と 戻り値 Float64, Int32
    # FORTRAN からの戻り値はないので Nothing を指定
    ccall(("test2_", "libtest.so"),
        Nothing,
        (Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Int32}),
        nrow, ncol, table, mean, cells)
    println("cells = $(cells[1])   mean = $(mean[1])") # スカラーとして取り出すときは [1] などとする
end

# FORTRAN プログラムへの引数

tbl = [1 2 3; 4 5 6]

# FORTRAN 関数副プログラム test1 を呼び出す FuncProg() の使用例
FuncProg(tbl)

# FORTRAN サブルーチン副プログラム test2 を呼び出す SubProg() の使用例
SubProg(tbl)

例えば以下のようにコマンドラインから実行する。

$ julia test.jl

結果は次のようになる。

     1.0    2.0    3.0
     4.0    5.0    6.0
mean = 3.5
     1.0    2.0    3.0
     4.0    5.0    6.0
cells = 6   mean = 3.5

 

コメント