裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

Julia の小ネタ--001 ベクトルの中から特定の値を持つ要素を除いたベクトルを得る

2021年02月22日 | ブログラミング

ベクトルの中から特定の値を持つ要素を除いたベクトルを得る

R でもやる,インプレースではない方法
元のベクトルを変えるなら a に代入。
別の変数に代入すれば a は元のまま。

a = [2, 1, 4, 3, 4, 5]
a = a[a .!= 4]
a # [2, 1, 3, 5]

以下は全てインプレース。
deleteat!(), filter!() は結果も返すので,別の変数に代入することもできる。
findfirst(), findlast(), findall() の第 1 引数は関数。

該当する要素を全て除いたベクトルを得る。

a = [2, 1, 4, 3, 4, 5]
deleteat!(a, findall(isequal(4), a))
a # [2, 1, 3, 5]

a = [2, 1, 4, 3, 4, 5]
filter!(x -> x != 4, a) # λ関数(無名関数)
a # [2, 1, 3, 5]

a = [2, 1, 4, 3, 4, 5]
filter!(!=(4), a) # != も関数
a # [2, 1, 3, 5]

a = [2, 1, 4, 3, 4, 5]
filter!(!isequal(4), a) # 関数はいくつか用意されている。
a # [2, 1, 3, 5]

該当する最初の要素だけを除く場合。

a = [2, 1, 4, 3, 4, 5]
deleteat!(a, findfirst(isequal(4), a))
a # [2, 1, 3, 4, 5]

該当する最後の要素だけを除く場合。

a = [2, 1, 4, 3, 4, 5]
deleteat!(a, findlast(isequal(4), a))
a # [2, 1, 4, 3, 5]

問題

整数ベクトル x の中に 5 が 3 個以上あるのは確かである。
2番目に出てくる 5 だけを削除したベクトルを求めなさい。

x = [3, 5, 2, 4, 6, 5, 8, 4, 5] のとき,[3, 5, 2, 4, 6, 8, 4, 5] を求める。

例解

x = [3, 5, 2, 4, 6, 5, 8, 4, 5]
deleteat!(x, findall(isequal(5), x)[2])

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia に翻訳--001 日本地図を描く(白地図のみ)

2021年02月22日 | ブログラミング

#=
Julia の修行をするときに,いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

日本地図を描く(白地図のみ)
http://aoki2.si.gunma-u.ac.jp/R/map.html

draw.map <- function(fn) # 境界線データのあるファイル名
{
        data <-read.table(fn, header=FALSE) # x, y 座標が組みになっている
        data[data[,1]==0 & data[,2]==0,] <- NA # x, y 座標が共に 0 であるのは,一連の境界線の終わりを意味する
        plot(data, type = "l", axes=FALSE, bty="n", xlab="", ylab="", asp=1) # (NA, NA) は,結線されない
}

ファイル名: drawmap.jl  関数名: drawmap

翻訳するときに書いたメモ

  • 列名のないファイルをデータフレームに読むときは,header=["x", "y"] のように名前を付けて読む。
  • 欠損値は missing で表される。
  • 元のファイルに欠損値があるときは,欠損値を表す文字列を missingstring="missing" で指定しないと,含まれる列が String 型になってしまう。
  • 元のファイルに欠損値がなくて,読み込み後のデータ処理で missing を割り当てる場合は,allowmissing!() で missing を許容するように変換する必要がある。
  • 行の選択のときなどで,複数の比較演算結果を論理演算するときは,比較演算は (data.x .== 0) .& (data.y .== 0) のように,( ) でくくること。
  • plot される x, y データ中に missing があると,その点は描画されない。これを利用することができる。
  • plot 引数で,よく使うもの
    線の色        linecolor=:black
    軸の描画      showaxis=false
    ティックマーク ticks=false
    グリッド      grid=false
    線の凡例      label=""
    アスペクト比   aspect_ratio=1

=#

using Plots, CSV, DataFrames
function drawmap(fn)
        data = CSV.read(fn, DataFrame, header=["x", "y"])
        allowmissing!(data)
        data[(data.x .== 0) .& (data.y .== 0), :] .= missing
        plot(data.x, data.y, linecolor=:black, showaxis=false,
             ticks=false, grid=false, label="", aspect_ratio=1)
end

drawmap("jpn.dat")
drawmap("map.data/toukyou")

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

逆引き辞書 Python, R, Julia (その 4)

2021年02月21日 | ブログラミング

メッシュグリッド

# Python
a, b = np.meshgrid([1,2,3], [4, 5, 6, 7])
a
# array([[1, 2, 3],
#        [1, 2, 3],
#        [1, 2, 3],
#        [1, 2, 3]])
b
# array([[4, 4, 4],
#        [5, 5, 5],
#        [6, 6, 6],
#        [7, 7, 7]])

# R
x = matrix(0, 4, 3)

col(x) # 列番号を要素とする行列
#      [,1] [,2] [,3]
# [1,]    1    2    3
# [2,]    1    2    3
# [3,]    1    2    3
# [4,]    1    2    3
#
row(x) # 行番号を要素とする行列
#  [,1] [,2] [,3]
# [1,]    1    1    1
# [2,]    2    2    2
# [3,]    3    3    3
# [4,]    4    4    4
#
row(x)+3 # 加工すればよい
#      [,1] [,2] [,3]
# [1,]    4    4    4
# [2,]    5    5    5
# [3,]    6    6    6
# [4,]    7    7    7

# Julia
function meshgrid(x, y)
    nr = length(y)
    nc = length(x)
    transpose(reshape(repeat(x, nr), nc, :)), reshape(repeat(y, nr), nr, :)
end

a, b = meshgrid([1,2,3], [4, 5, 6, 7])
# a : [1 2 3; 1 2 3; 1 2 3; 1 2 3]
# b : [4 4 4 4; 5 5 5 5; 6 6 6 6; 7 7 7 7]

配列の連結

# Python
a = np.array([1.0, 2.0, 3.0])           # [ 1.  2.  3.]
b = np.array([5.0, 6.0, 7.0])           # [ 5.  6.  7.]
np.vstack((a, b))                       # [[ 1.,  2.,  3.], # np でも同じ
                                        #  [ 5.,  6.,  7.]]
np.hstack((a, b))                       # [ 1.,  2.,  3.,  5.,  6.,  7.]
#
x = np.array([[1, 2, 3], [4, 5, 6]])
y = np.array([[10, 11, 12], [20, 21, 22]])

np.vstack((x, y))                       # [[ 1,  2,  3], # rbind
                                        #  [ 4,  5,  6],
                                        #  [10, 11, 12],
                                        #  [20, 21, 22]]

np.hstack((x, y))                       # [[ 1,  2,  3, 10, 11, 12], # cbind
                                        #  [ 4,  5,  6, 20, 21, 22]]
#
np.concatenate((a, b))                  # [ 1.,  2.,  3.,  5.,  6.,  7.]

np.concatenate(([a], [b]))              # [[ 1.,  2.,  3.],
                                        #  [ 5.,  6.,  7.]]
#
np.concatenate((x, y), axis=0)  # [[ 1,  2,  3],
                                        #  [ 4,  5,  6],
                                        #  [10, 11, 12],
                                        #  [20, 21, 22]]

np.concatenate((x, y), axis=1)  # [[ 1,  2,  3, 10, 11, 12],
                                        #  [ 4,  5,  6, 20, 21, 22]]
#
np.c_[a, b]                             # [[ 1.,  5.],              R: cbind()
                                        #  [ 2.,  6.],
                                        #  [ 3.,  7.]]

np.r_[a, b]                             # [ 1.,  2.,  3.,  5.,  6.,  7.]     R: rbind()
#
np.c_[x, y]                             # [[ 1,  2,  3, 10, 11, 12],
                                        #  [ 4,  5,  6, 20, 21, 22]]

np.r_[x, y]                             # [[ 1,  2,  3],
                                        #  [ 4,  5,  6],
                                        #  [10, 11, 12],
                                        #  [20, 21, 22]]

# R
a = c(1, 2, 3)
b = c(5, 6, 7)
x = matrix(1:6, ncol=3, byrow=TRUE)
y = matrix(c(10, 11, 12, 20, 21, 22), ncol=3, byrow=TRUE)

rbind(a, b)
cbind(a, b) # np.c_ に相当
c(a, b)
append(a, b)

rbind(x, y)
cbind(x, y)

# Julia
a = [1, 2, 3]
b = [5, 6, 7]
x = [1 2 3; 4 5 6]
y = [10 11 12; 20 21 22]

vcat(a, b) # [1, 2, 3, 4, 5, 6]
hcat(a, b) # [1 5; 2 6; 3 7]

vcat(x, y)
# 1   2   3
# 4   5   6
# 10  11  12
# 20  21  22

hcat(x, y)
# 1  2  3  10  11  12
# 4  5  6  20  21  22

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

逆引き辞書 Python, R, Julia (その 3)

2021年02月21日 | ブログラミング

次元数の変更

# Python
a = np.reshape(np.arange(9), (3, 3)) # array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
b = np.reshape(np.arange(9), (3, -1))
c = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
d = np.reshape(c, (9,))              # array([0, 1, 2, 3, 4, 5, 6, 7, 8])

# Julia
vec1 = collect(1:6);   # [1, 2, 3, 4, 5, 6]
reshape(vec1, 2, 3);   # [1 3 5; 2 4 6] # 2x3 行列
reshape(vec1, (2, 3)); # [1 3 5; 2 4 6] # 2x3 行列
reshape(vec1, 2, :);   # [1 3 5; 2 4 6]  # 2x3 行列
reshape(vec1, 6, :);   # 6x1 行列

リサイズ(インプレースで変更)

# Python
d = np.arange(4); d.resize((6,)); d # [0, 1, 2, 3, 0, 0]
a = np.array([[1, 2, 3], [4, 5, 6]])
a.resize((6,))                      # [1, 2, 3, 4, 5, 6]; 行列からベクトル
#
e = np.arange(4)
e.resize((4, 2)) # [[0, 1], [2, 3], [0, 0], [0, 0]] # ベクトルから行列

行列をベクトルに変更,スカラーさえも

# Python
x = np.array([[1,2,3,4],[5,6,7,8]]); np.ravel(x)
np.ravel(1) # [1]

# R
x = matrix(1:8, byrow=TRUE, ncol=4)
c(t(x)) # 列優先
c(x)    # 行優先なら c(x)

# Julia
x = [1 2 3 4; 5 6 7 8]
vec(x') # [1, 2, 3, 4, 5, 6, 7, 8]
vec(x)  # [1, 5, 2, 6, 3, 7, 4, 8]

a = 1
[a]

単位行列

# Python
np.identity(3)
np.eye(3)
np.eye(3, 5) # 3x5 行列

# R
diag(3)

# Julia
using LinearAlgebra
Array{Int}(I, 3, 3)    # 3x3 行列
Array{Float64}(I, 3, 3)
Array{Float64}(I, 2, 3) # 2x3 行列

対角行列

# Python
np.diag([1, 3, 9])

# R
diag(c(1, 3, 9)) # 対角要素が c(1, 3, 9)
diag(3, 5)       # 対角要素が 3 の 5x5 行列

# Julia
diagm(0 => 1:3)
diagm(0 => [1, 3, 9])
diagm(1 => [1, 3, 9])                   # 対角要素の 1 つ上。従って 4x4 行列
diagm(-2 => [1, 3, 9])                  # 対角要素の 2 つ下。従って 5x5 行列
diagm(1 => [1, 3, 9], -2 => [2, 6, 13]) # 複数指定も可
Diagonal([1, 2, 3])
Diagonal(ones(3))
Diagonal(ones(Int, 3))

対角成分へ代入

# Python
r = np.reshape(np.arange(1, 10), (3, 3))
diag2 = np.eye(3) == 1
r[diag2] = 999 # [[999, 2, 3], [4, 999, 6], [7, 8, 999]]

# R
r = matrix(1:9, byrow=TRUE, 3)
diag(r) = 999

Julia には用意されていないようなので,独自関数を定義する。

function diag!(a, b)
    [a[i, i] = b[i] for i=1:length(b)]
end

a = [1 2 3; 4 5 6; 7 8 9]
b = [10 20 30]

diag!(a, b)
a

対角成分の取り出し

# Python
r = np.reshape(np.arange(1, 10), (3, 3))
np.diag(r)      # [1, 5, 9]
np.diag(r, 1)   # [2, 6] 対角線の 1 つ上
np.diag(r, -2)  # [7] 対角線の 2 つ下

# R
r = matrix(1:9, byrow=TRUE, 3)
diag(r)         # c(1, 5, 9)

# Julia
r = [1 2 3; 4 5 6; 7 8 9]
diag(r)         # [1, 5, 9]
diag(r, 1)      # [2, 6]
diag(r, -2)     # [7]

三角行列

要素が 1 の下三角行列の作成

# Python
np.tri(3)
np.tri(2, 5, 2)
np.tri(3, 5, -1)

# R
a = matrix(0, 3, 3)
a[lower.tri(a, diag=TRUE)] = 1 # np.tri(n) に相応


下三角行列(上三角行列)の抽出

# Python
a = np.reshape(np.arange(1, 10), (3,3))
# array([[1, 2, 3],
#        [4, 5, 6],
#        [7, 8, 9]])
np.tril(a)
# array([[1, 0, 0],
#        [4, 5, 0],
#        [7, 8, 9]])
np.tril(a, -1)         # 対角成分より一つ下(対角を含まない)
# array([[0, 0, 0],
#        [4, 0, 0],
#        [7, 8, 0]])
np.triu(a)
# array([[1, 2, 3],
#        [0, 5, 6],
#        [0, 0, 9]])
np.triu(a, -1)         # 対角成分より一つ下(対角を含まない)
# array([[1, 2, 3],
#        [4, 5, 6],
#        [0, 8, 9]])

# R
ifelse(lower.tri(a, diag=TRUE), a, 0)
ifelse(lower.tri(a), a, 0)            # 対角成分より一つ下
ifelse(upper.tri(a, diag=TRUE), a, 0)
ifelse(upper.tri(a), a, 0)            # 対角成分より一つ上

# Julia
a = [1 2 3; 4 5 6; 7 8 9]
[i <= j ? a[j, i] : 0 for j = 1:n, i = 1:n] # 不等式部分でどこから取り出すかを指定
[i < j ? a[j, i] : 0 for j = 1:n, i = 1:n]
[i >= j ? a[j, i] : 0 for j = 1:n, i = 1:n]
[i >= j-1 ? a[j, i] : 0 for j = 1:n, i = 1:n]

下三角行列の要素をベクトルとして抽出

# Python
a = np.reshape(np.arange(1, 10), (3,3))
# array([[1, 2, 3],
#        [4, 5, 6],
#        [7, 8, 9]])
#
b = np.tril(np.ones_like(a)) > 0
# array([[ True, False, False],
#        [ True,  True, False],
#        [ True,  True,  True]])
a[b]    # [1, 4, 5, 7, 8, 9]
#
b2 = np.triu(np.ones_like(a)) > 0
# array([[ True,  True,  True],
#        [False,  True,  True],
#        [False, False,  True]])
a[b2]    # [1, 2, 3, 5, 6, 9]
#
c = np.tril(np.ones_like(a), -1) > 0
# array([[False, False, False],
#        [ True, False, False],
#        [ True,  True, False]])
a[c]    # [4, 7, 8]
#
c2 = np.triu(np.ones_like(a), -1) > 0
# array([[ True,  True,  True],
#        [ True,  True,  True],
#        [False,  True,  True]])
a[c2]    # array([1, 2, 3, 4, 5, 6, 8, 9])

# R
a = matrix(1:9, 3, byrow=TRUE)
#
a[lower.tri(a, diag=TRUE)] # c(1, 4, 7, 5, 8, 9); 列優先になる
a[upper.tri(a, diag=TRUE)] # c(1, 2, 5, 3, 6, 9)
#
a[lower.tri(a)] # c(4, 7, 8)
a[upper.tri(a)] # c(2, 3, 6)

# Julia
a = [1 2 3; 4 5 6; 7 8 9]
using LinearAlgebra
a[tril(trues(size(a)))] # [1, 4, 7, 5, 8, 9]; 列優先になる
a[triu(trues(size(a)))] # [1, 2, 5, 3, 6, 9]

下三角行列(上三角行列)への代入

# Python
r = np.reshape(np.arange(1, 10), (3, 3))
# array([[1, 2, 3],
#        [4, 5, 6],
#        [7, 8, 9]])
p = 100 * r
# array([[100, 200, 300],
#        [400, 500, 600],
#        [700, 800, 900]])
lower = np.tri(3) == 1
r[lower] = p[lower]
print(r)
# array([[100,   2,   3],
#        [400, 500,   6],
#        [700, 800, 900]])

# R
r = matrix(1:9, 3, byrow=TRUE)
#      [,1] [,2] [,3]
# [1,]    1    2    3
# [2,]    4    5    6
# [3,]    7    8    9
p = 100 * r
#      [,1] [,2] [,3]
# [1,]  100  200  300
# [2,]  400  500  600
# [3,]  700  800  900
ifelse(lower.tri(p, diag=TRUE), p, r)
#      [,1] [,2] [,3]
# [1,]  100    2    3
# [2,]  400  500    6
# [3,]  700  800  900
ifelse(upper.tri(p, diag=TRUE), p, r)
#      [,1] [,2] [,3]
# [1,]  100  200  300
# [2,]    4  500  600
# [3,]    7    8  900

# Julia
r = [1 2 3; 4 5 6; 7 8 9]
# 1  2  3
# 4  5  6
# 7  8  9
p = 100 * r
# 100  200  300
# 400  500  600
# 700  800  900
n = size(p, 1)
[i > j ? r[i, j] : p[i, j] for j = 1:n, i = 1:n]
# 100    4    7
# 200  500    8
# 300  600  900
[i < j ? r[i, j] : p[i, j] for j = 1:n, i = 1:n]
# 100  400  700
#   2  500  800
#   3    6  900

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

逆引き辞書 Python, R, Julia (その 2)

2021年02月21日 | ブログラミング

初期化した配列(ベクトル,行列)

ゼロベクトル

# Python
np.zeros(3)

# R
integer(3)
numeric(3)
logical(3)

# Julia
zeros(3)       # [0.0, 0.0, 0.0]
zeros(Int, 3)  # [0, 0, 0]
zeros(Bool, 3) # [0 0 0]
falses(3)      # [0 0 0]

ゼロマトリックス

# Python
np.zeros((3, 4))

# R
matrix(0, 3, 4)

# Julia
zeros(3, 4)         # [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
zeros(Int, 3, 4)    # [0 0 0; 0 0 0; 0 0 0]
zeros(Bool, 3, 4)   # [0 0 0; 0 0 0; 0 0 0]

1 ベクトル

# Python
np.ones(3)

# R
rep(1, 3)

# Julia
ones(3)         # [1.0, 1.0, 1.0]
ones(Int, 3)    # [1, 1, 1]
ones(Bool, 3)   # [1, 1, 1]
trues(3)        # [1, 1, 1]

1 マトリックス

# Python
np.ones((3, 4))

# R
matrix(1, 3, 4)

# Julia
ones(3, 4)      # [1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0]

ones(Int, 3, 4) # [1 1 1; 1 1 1; 1 1 1]

一定値で初期化するベクトル

# Python
np.repeat(123, 3)
np.full(3, 123)

# R
rep(123, 3)

# Julia
fill(123, 3) # [123, 123, 123]

一定値で初期化するマトリックス

# Python
np.array([[123]*4 for i in range(3)])
np.full((3, 4), 123)

# R
matrix(123, 3, 4)

# Julia
fill(123, 3, 4) # [123 123 123 123; 123 123 123 123; 123 123 123 123]

初期化しないベクトル,行列

# Python
np.empty(5)
np.empty(5, 7)

# Julia
V = Array{Float64, 1}(undef, 5)     # 型の指定は Int64, Int, Bool など
A = Array{Float64, 2}(undef, 5, 7)

既存のベクトル,行列と同じサイズのベクトル,行列

# Python
a = np.ones(5)
np.ones_like(a)
np.zeros_like(a)
np.empty_like(a)
#
b = np.ones((2, 3))
np.ones_like(b)
np.zeros_like(b)
np.empty_like(b)

# Julia
a = ones(5)
a2 = similar(a) # サイズだけが同じで初期化されていないベクトル
b = zeros(2, 3)
b2 = similar(b) # サイズだけが同じで初期化されていない行列

ソート

# Python
a = [4, 1, 3, 5, 2]
sorted(a)
sorted(a, reverse=True)
#
a.sort()
a.sort(reverse=True)   # インプレースでソート
#
b = np.array([[4, 3, 5], [3, 2, 1]])
b.sort(axis=1)         # インプレース・次元指定でソート
#
c = np.sort(b, axis=1) # np.sort はインプレースではない

# R
a = c(4, 1, 3, 5, 2)
sort(a)
sort(a, decreasing=TRUE)
#
x = matrix(c(4, 3, 5, 3, 2, 1), byrow=TRUE, ncol=3)
t(apply(x, 1, sort))  # 次元指定でソート

# Julia
a = [4, 1, 3, 5, 2];
sort(a)
sort(a, rev=true)
#
sort!(a) # インプレースでソート
sort!(a, rev=true) # インプレースで逆順ソート

逆転

# Python
a = [4, 1, 3, 2, 5]
a.reverse()
a[::-1]
#
reversed(a) # 逆順のイテレータ

# R
a = c(4, 1, 3, 2, 5)
rev(a)
a[seq(length(a), 1, by=-1)]

# Julia
a = [4, 1, 3, 2, 5]
reverse(a)      # [5, 2, 3, 1, 4]
a[end:-1:begin] # [5, 2, 3, 1, 4]

小さい順の要素のインデックス

x[o] とすると x がソートされるような o

# Python
x = np.array([2, 1, 2, 3, 4, 5, 4, 3, 4, 6])
o = x.argsort() # array([1, 0, 2, 3, 7, 4, 6, 8, 5, 9]) # 0 から始まる
x[o]            # array([1, 2, 2, 3, 3, 4, 4, 4, 5, 6])
(-x).argsort()  # array([9, 5, 4, 6, 8, 3, 7, 0, 2, 1])

# R
x = c(2, 1, 2, 3, 4, 5, 4, 3, 4, 6)
o = order(x)                           # c(2, 1, 3, 4, 8, 5, 7, 9, 6, 10) # 1 から始まる
x[o] = c(1, 2, 2, 3, 3, 4, 4, 4, 5, 6)
order(x, decreasing=TRUE)              # c(10, 6, 5, 7, 9, 4, 8, 1, 3, 2)

# Julia
x = [2, 1, 2, 3, 4, 5, 4, 3, 4, 6]
o = sortperm(x) # [2, 1, 3, 4, 8, 5, 7, 9, 6, 10]
x[o]            # [1, 2, 2, 3, 3, 4, 4, 4, 5, 6]

逆順ソートに対応するベクトル ro を返す。
ro = sortperm(x, rev=true) # [10, 6, 5, 7, 9, 4, 8, 1, 3, 2]
x[ro] # [6, 5, 4, 4, 4, 3, 3, 2, 2, 1]

等差数列

公差=1

# Python
list(range(5))
np.array(range(5))
np.arange(5)

# R
0:4
seq(0, 4)

# Julia
Vector(0:4)                     # [0, 1, 2, 3, 4]
collect(0:4)                    # [0, 1, 2, 3, 4]
range(0, 4, step=1)             # 0:1:4 # イテレータ
collect(range(0, 4, step=1))    # [0, 1, 2, 3, 4]

公差≠1

# Python
np.arange(10, 13, 0.5) # array([10., 10.5, 11., 11.5, 12., 12.5])

# R
seq(10, 12.5, 0.5)     # c(10, 10.5, 11, 11.5, 12, 12.5)

# Julia
Vector(10:0.5:12.5)                 # [10.0, 10.5, 11.0, 11.5, 12.0, 12.5]
collect(10:0.5:12.5)                # [10.0, 10.5, 11.0, 11.5, 12.0, 12.5]
range(10, 12.5, step=0.5)           # 10.0:0.5:12.5 # イテレータ
collect(range(10, 12.5, step=0.5))  # [10.0, 10.5, 11.0, 11.5, 12.0, 12.5]

長さ指定

# Python
np.linspace(1, 6, 5) # [1., 2.25, 3.5, 4.75, 6.]

# R
seq(1, 6, length=5)  # c(1, 2.25, 3.5, 4.75, 6)

# Julia
collect(1:(6-1)/4:6)            # [1.0, 2.25, 3.5, 4.75, 6.0]
range(1, 6, length=5)           # 1.0:1.25:6.0 # イテレータ
collect(range(1, 6, length=5))  # [1.0, 2.25, 3.5, 4.75, 6.0]

等比数列

# Python
np.logspace(0, 2, num = 3, dtype=int) # array([1, 10, 100]) i.e. 10^0, 10^1, 10^2
np.logspace(0, 5, num = 6, base=3)    # array([1., 3., 9., 27., 81., 243.])

# R
10^(0:2) # c(1, 10, 100)
3^(0:5)  # c(1, 3, 9, 27, 81, 243)

10 .^ (0:2) # [1, 10, 100]
3 .^ (0:5)  # [1, 3, 9, 27, 81, 243]

乱数

乱数の種

# Python
import random
random.seed(123)
#
np.random.seed(123)

# R
set.seed(123)

# Julia
using Random
Random.seed!(123)

一様乱数

# Python
random.random()                  # [0, 1) の一様乱数を 1 個
random.uniform(min, max)         # [min, max) の一様乱数
random.randint(min, max)         # [min, max] の整数一様乱数
np.random.rand(n)                # [0, 1) の一様乱数を n 個
np.random.randint(min, max, n)   # [min, max) の整数一様乱数を n 個

# R
runif(n, min, max)               # [min, max) の一様乱数を n 個
sample(min:max, n, replace=TRUE) # [min, max] の整数一様乱数 を n 個
sample(max, n, replace=TRUE)     # [1, max] の整数一様乱数 を n 個

# Julia
rand(n)         # Float64 の [0, 1) の一様乱数
rand(Int8, n)   # Int8 の一様乱数を n 個
rand(Int32, n)  # Int32 の一様乱数を n 個

無作為抽出

# Python
random.randint(min, max)            # [min, max] の整数一様乱数 1 個
np.random.randint(min, max, n)      # [min, max) の整数一様乱数 n 個
np.random.sample(["A", "B", "C"], n, replace=True)

# R
sample(min:max, n, replace=TRUE)    # [min, max] の整数一様乱数 n 個
sample(c("A", "B", "C"), n, replace=TRUE)

# Julia
n = 3; min0=1; max0=6
rand(min0:max0, n)                  # min0:max0 から n 個を無作為抽出
rand(["A", "B", "C"], n)            # ベクトル要素を無作為抽出
rand(('x', 'y', :z), n)             # ベクトル要素を無作為抽出
rand(2, 3)                          # 2x3 の一様乱数行列

ベクトルをシャッフル(重複なしの無作為抽出)

# Python
np.shuffle(x)

# R
sample(x)

# Julia
x = collect(1:10)
shuffle(x)

正規乱数

# Python
random.normalvariate(mean, std)     # 平均値 = mean, 標準偏差 = std の正規乱数 1 個
random.gauss(mean, std)             # 平均値 = mean, 標準偏差 = std の正規乱数 1 個
np.random.normal(mean, std, n)      # 平均値 = mean, 標準偏差 = std の正規乱数 n 個
np.random.normal(mean, std, (n, m)) # nxm の正規乱数行列

# R
rnorm(n, mean, sd)                  # 平均値 = mean, 標準偏差 = sd の正規乱数 n 個
matrix(rnorm(n*m, mean, sd), n, m)  # nxm の正規乱数行列

# Julia
Random.seed!(123)
randn(3)                            # 長さ 3 の正規乱数ベクトル(平均値 = 0,標準偏差 = 1 に固定)
randn(3, 4)                         # 3x4 の正規乱数行列

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

逆引き辞書 Python, R, Julia (その 1)

2021年02月21日 | ブログラミング

in 演算子

# Python
a = range(1, 11)
3 in range(1, 11) # True
3 in a            # True
not(3 in a)       # False
20 not in a       # True

# R
a = 1:10
3 %in% a        # TRUE
!(3 %in% a)     # FALSE
20 %in% a       # FALSE

# Julia
a = 1:10        # 1:10
b = collect(a)  # [1, 2, ..., 9, 10]
typeof(a)       # UnitRange{Int64}
typeof(b)       # Array{Int64,1}
3 in a          # true
3 in b          # true
!(3 in a)       # false
20 in a         # false

どこにあるか

# Python
a = [3, 1, 0, 2, 4, 0, 3, 5, 4, 2]
a.index(4)                 # 4
#
b = np.array([3, 2, 1, 2, 3, 4, 1])
c = np.where(b == 1)       # (array([2, 6]),)
c[0]                       # array([2, 6])
c[0][0]                    # 2
c[0][1]                    # 6

# R
a = c(3, 1, 0, 2, 4, 0, 3, 5, 4, 2)
which(a == 4)               # 5, 9

# Julia
a = [3, 1, 0, 2, 4, 0, 3, 5, 4, 2]
findfirst(==(4), a)         # 5
findfirst(isequal(4), a)    # 5
#
findfirst(>(3), a)          # 5
findall(iseven, a)          # [3, 4, 5, 6, 9, 10]
findall(isodd, a)           # [1, 2, 7, 8]
#
findfirst(==(3), a)         # 1
findlast(==(3), a)          # 7
findall(==(3), a)           # [1, 7]

いくつあるか

# Python
a = [3, 1, 0, 2, 4, 0, 3, 5, 4, 2]
a.count(3)                  # 2

# R
a = c(3, 1, 0, 2, 4, 0, 3, 5, 4, 2)
length(which(a == 3))       # 2

# Julia
a = [3, 1, 0, 2, 4, 0, 3, 5, 4, 2]
count(==(3), a)             # 2
count(isequal(3), a)        # 2

三項演算子

# Python
i = 6
a = "even" if i % 2 == 0 else "odd"     # 'even'

# R
i = 6
a = if (i %% 2 == 0) "even" else "odd"  # "even"

# Julia
i = 6
a = i % 2 == 0 ? "even" : "odd"         # "even"

条件により実行される演算を選択する

# Python
import numpy as np
a = np.array([2, 3, 5, 6, 8])
b = np.where(a % 2 == 0, a*10, 0)     # array([20,  0,  0, 60, 80])

# R
a = c(2, 3, 5, 6, 8)
b = ifelse(a %% 2 == 0, a*10, 0)      # c(20, 0, 0, 60, 80)

# Julia
a = [2, 3, 5, 6, 8];
b = map(x -> x % 2 == 0 ? 10x : 0, a) # [20, 0, 0, 60, 80]

リスト(配列)同士の連結

# Python
a = [1, 2, 3]
b = [11, 12]
x = a + b           # [1, 2, 3, 11, 12]; インプレースではない
#
import operator as op
y = op.concat(a, b) # [1, 2, 3, 11, 12]
#
a.extend(b)         # [1, 2, 3, 11, 12]; a にインプレースで連結
#
a2 = np.array([1, 2, 3])
b2 = np.array([11, 12])
x = np.hstack((a2, b2)) # array([1, 2, 3, 11, 12]); インプレースではない

# R
a = c(1, 2, 3)
b = c(11, 12)
x = c(a, b)         # c(1, 2, 3, 11, 12); インプレースではない
y = append(a, b)    # c(1, 2, 3, 11, 12); インプレースではない

# Julia
a = [1, 2, 3];
b = [11, 12];
x = vcat(a, b)      # [1, 2, 3, 11, 12]; インプレースではない
append!(a, b)       # [1, 2, 3, 11, 12]; a にインプレースで連結

リストに変数を連結

# Python
a = [1, 2, 3]
b = 5
x = a + [b]         # [1, 2, 3, 5]; インプレースではない
a.append(b)         # [1, 2, 3, 5]; a にインプレースで連結
#
a2 = np.array([1, 2, 3])
b = 5
x = np.hstack((a2, b)) # [1, 2, 3, 5]; インプレースではない

# R
a = c(1, 2, 3)
b = 5
x = c(a, b)         # c(1, 2, 3, 5); インプレースではない
y = append(a, b)    # c(1, 2, 3, 5); インプレースではない

# Julia
a = [1, 2, 3]
b = 5
x = vcat(a, b)      # [1, 2, 3, 5]; インプレースではない
append!(a, b)       # [1, 2, 3, 5]; a にインプレースで追加

リストに変数を挿入

# Python
a = [1, 2, 3]
b = 5
a.insert(0, b)      # [5, 1, 2, 3]; a にインプレースで挿入

# R
a = c(1, 2, 3)
b = 5
x = append(a, b, 0) # c(5, 1, 2, 3); インプレースではない

# Julia
a = [1, 2, 3]
b = 5
insert!(a, 1, b)    # [5, 1, 2, 3]; a にインプレースで挿入

要素の全削除

# Python
a = [1, 2, 3]
a.clear()
a                 # []; 要素がなくなるのであって,a 自体がなくなるわけではない

# R
a = c(1, 2, 3)
a = NULL
a                 # NULL; 要素がなくなるのであって,a 自体がなくなるわけではない

# Julia
a = [1, 2, 3]
empty!(a)
a                 # Int64[]; 要素がなくなるのであって,a 自体がなくなるわけではない

要素の削除

# Python
a = [1, 2, 3]
del a[1]          # インデックスを指定して要素を削除
a                 # [1, 3]
#
a = [1, 2, 3]
a.remove(3)       # 指定する値を持つ要素をインプレースで削除
a                 # [1, 2]

# R
a = c(1, 2, 3)
a = a[-2]         # インプレースではないため,結果を反映させるためには代入が必要

# Julia
a = [1, 2, 3];
deleteat!(a, 2)   # a からインプレースで削除
a # [1, 3]

要素のポップ

# Python
a = [3, 2, 1, 4, 3]
a.pop()             # 3; 取り出された要素
a                   # [3, 2, 1, 4]
#
a = [3, 2, 1, 4, 3 ]
a.pop(1)            # 2; 取り出された要素
a                   # [3, 1, 4, 3]

# R 以下の関数では,取り出した要素は把握できない
pop = function(a, i=length(a)) a[-i]
#
a = c(3, 2, 1, 4, 3)
a = pop(a)
a                   # c(3, 2, 1, 4)
#
a = c(3, 2, 1, 4, 3)
a = pop(a, 3)
a                   # c(3, 2, 4, 3)

# Julia
a = [3, 2, 1, 4, 3];
pop!(a)             # 3; 取り出された要素
a                   # [3, 2, 1, 4] # a からインプレースで最後の要素を取り出す
#
a = [3, 2, 1, 4, 3];
popfirst!(a)        # 3; 取り出された要素
a                   # [2, 1, 4, 3] # a からインプレースで最初の要素を取り出す

要素のプッシュ

# Python
a = [3, 2, 1, 4, 3]
a.append(88)        # 最後の要素としてインプレースで付け加える(insert を使ってもよい)
a.insert(3, 99)     # [3, 2, 1, 99, 4, 3, 88]; 3 番目の要素としてインプレースで挿入

# R
a = c(3, 2, 1, 4, 3)
a = append(a, 88)        # インプレースではないので,代入しないと変化なし
a = append(a, 99, 3)     # インプレースではないので,代入しないと変化なし
a                        # c(3, 2, 1, 99, 4, 3, 88)

# Julia
a = [3, 2, 1, 4, 3];
pushfirst!(a, 10)        # a にインプレースで先頭に要素を押し込む
push!(a, 20)             # a にインプレースで最後尾に要素を押し込む
a                        # [10, 3, 2, 1, 4, 3, 20]

リストの繰り返し

# Python
a = [1, 2, 3]
b = a*3                  # [1, 2, 3, 1, 2, 3, 1, 2, 3]
#
np.tile(a, 3)            # array([1, 2, 3, 1, 2, 3, 1, 2, 3])
#
np.repeat(a, 3)          # array([1, 1, 1, 2, 2, 2, 3, 3, 3])
#
np.repeat(np.array([1,2]), [2,3])    # array([1, 1, 2, 2, 2])

# R
rep(c(1, 2, 3), 3)       # c(1, 2, 3, 1, 2, 3, 1, 2, 3)
#
rep(c(1, 2, 3), each=3)  # c(1, 1, 1, 2, 2, 2, 3, 3, 3)
#
rep(c(1, 2), c(2, 3))    # c(1, 1, 2, 2, 2)

# Julia
repeat([1, 2, 3], 3)     # [1, 2, 3, 1, 2, 3, 1, 2, 3]
#
repeat([1 2 3], inner=(1, 2))               # [1 1 2 2 3 3]; 1x6 行列
reshape(repeat([1 2 3], inner=(1,2)), 6, :) # [1, 1, 2, 2, 3, 3]; 6x1 行列
#
# 関数定義
repeateach(x, n) = vec([j for i = 1:n, j in x])
a = repeateach([1, 2, 3], 3)        # [1, 1, 1, 2, 2, 2, 3, 3, 3]
a = repeateach(["a", "b", "c"], 2) # ["a", "a", "b", "b", "c", "c"]

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

数学問題を解きましたということなんだが

2021年02月19日 | ブログラミング

曲線 y = x^3 上に,点 A(-1, -1),B(1, 1),P(t, t^3) を置く。ただし,0 < t < 1。
直線 AP と x 軸のなす角を α,直線 PB と x 軸のなす角を β とする。
tan α,tan β を t を用いて表せ。ただし,0 < α < π/2, 0 < β < π/2 とする。
tan∠APB を t を用いて表せ。
tan∠APB を最小にする t の値を求めよ。

using Plots
x = -1.1:0.01:1.1
y = x .^ 3;
p = plot(x, y, xlabel="\$x\$", ylabel="\$y\$", label="",
   xlims=(-1.5, 1.5), ylims=(-1.1, 1.1),
   size=(600, 400), aspect_ratio=1);
x2 = [-1, 0.5, 1];
y2 = x2 .^ 3;
plot!(x2, y2, label="", linecolor=:black);
scatter!(x2, y2, label="", markercolor=:black);
annotate!(-0.65, -1, text("\$A(-1, -1)\$", 10));
annotate!(1.2, 1, text("\$B(1, 1)\$", 10));
annotate!(0.7, 0.4^3, text("\$P(t, t^3)\$", 10));
annotate!(-0.25, 0.15, text("\$y = x^3\$", :blue, 15));
display(p)

例のごとく,詰めが甘い。

問題 1 は,図を描いて,tangent の定義をそのまま書けばよいと。

α = atan((1-t^3)/(1-t))
β =
atan((1+t^3)/(1+t))

simplify() してもよいけど...

問題2 は,求める θ を 問題 1 の α,β を用いて表すだけ(ちょっと,混乱するけど)

その後は,手計算でやるんだろうけど,SymPy を使ってずるをする。

using SymPy

@syms θ α β t

# 問題 1
# α = atan((1-t^3)/(1-t))  後に示すように,simplify してもよいけどそのまま
# β =
atan((1+t^3)/(1+t))

# 問題 2
θ = π + atan((1-t^3)/(1-t)) - atan((1+t^3)/(1+t))
# 導関数
diff(θ)
# 導関数 = 0 つまり,接線が 0 になるところを求める
solve(diff(θ))
# 以下の内 0 < t < 1 を満たすのは
sqrt(6)/3
# 上のでもよいが,単に数値計算して
0.8164965809277259

# α,βを簡約して,以下のようにも

simplify((t^3+1)/(t+1))
# α = t^2-t+1
simplify((1-t^3)/(1-t))

# β = t^2-t+1
θ = π + atan(t^2-t+1) - atan(t^2+t+1)
diff(θ)

a = solve(diff(θ))
a[1].evalf()
# 0.816496580927726

a[2].evalf()
# 0.816496580927726

数値計算で解を示すなら,二分法かニュートンラフソン法で求解するのが常套手段かな?

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia でコンソールから入力

2021年02月18日 | ブログラミング

コンソールから入力

基本は,

  • readline() で文字列を読む。
  • 複数のフィールドがあるなら,split(文字列)でフィールドを分ける。
  • 分割された文字列を,整数に変換するには parse(Int, 文字列),実数に変換するには parse(Float64, 文字列)。
  • 複数の整数,実数変換は parse.() のようにドット記法を使う。
  • 変換された要素を配列(ベクトル)に代入,または複数の変数に代入する。

文字列として 1 行入力

入力:"abcdefg 12345"

s = readline() # "abcdefg 12345"

区切り文字で区切って複数の文字列として配列(ベクトル)へ入力

デフォルトでは空白区切り

入力:"abcdefg 12345"

s = split(readline()) # 要素数 2 のベクトル: "abcdefg", "123"

カンマで区切るとき

入力:"123,456,789"

s = split(readline(), ",") # 要素数 3 のベクトル: "123", "456", "789"

複数の変数へ入力

入力:"123 456 789"

a, b, c = split(readline()) # 3 変数へ入力: "123", "456", "789"
a # "123"
b # "456"
c # "789"

1 行に整数 1 つ

入力:"12345"

n = parse(Int, readline()) # 12345

1 行に実数 1 つ

入力:"12345"

f = parse(Float64, readline()) # 12345.0

入力:"12.345"

f = parse(Float64, readline()) # 12.345

複数の整数を複数の変数へ入力

入力:"12345 67890 11 22"

i, j, k, l = parse.(Int, split(readline())) # i = 12345, j = 67890, k = 11, l = 22

複数の実数を複数の変数へ入力

入力:"0.12 4e2 -3.5e-3"

a, b, c = parse.(Float64, split(readline())) # a = 0.12, b = 400.0, c = -0.0035

複数の整数を配列(ベクトル)へ入力

入力:"12345 67890 11 22"

n = parse.(Int, split(readline())) # 要素数 4 のベクトル: 12345, 67890, 11, 22

複数の整数を配列(ベクトル)へ入力する(リスト内包表記; 非推奨)

入力:"23 45 67"

m = [parse(Int, x) for x in split(readline())] # 23, 45, 67

複数の実数を配列(ベクトル)へ入力

入力:"0.12 4e2 -3.5e-3"

f = parse.(Float64, split(readline())) # 要素数 3 のベクトル: 0.12, 400.0, -0.0035

複数の実数を配列(ベクトル)へ入力する(リスト内包表記; 非推奨)

入力:"23 45 67"

g = [parse(Float64, x) for x in split(readline())] # 23.0, 45.0, 67.0

読み込む変数の型が複数あるとき(文字列,整数,実数が混ざっているとき)

入力:"abc 123 45.67"

a, b, c = split(readline())
a # "abc"
b = parse(Int, b) # 123
c = parse(Float64, c) # 45.67

文字列を入力する(split() しない)と文字に分解されて集合に読み込まれる

入力:"qwretudhddhdhshkfh"

s = Set(readline()) # 'f', 'w', 'd', 'e', 'h', 's', 'u', 'r', 't', 'k', 'q'

文字列として集合に読み込ませるには,split() する

入力:"abc 123 *&* @@"

s2 = Set(split(readline())) # "abc", "*&*", "@@", "123"

複数行からの入力

1 行に 1 個ずつの整数を n 個(n 行から)入力する

n = 4 # これも入力するなら n = parse(Int, readline())
入力:
"123"
"456"
"789"
"111"

m = [parse(Int, readline()) for s in 1:n] # 123, 456, 789, 111

実数の入力なら parse(Float64, readline()),
文字列そのままを入力するなら readline() のまま

1 行に複数個の整数,別々の配列(ベクトル)の要素として入力

入力:
"1 2 3"
"4 5 6"
"7 8 9"
"0 1 2"

n = 4
a = zeros(Int, n)
b = zeros(Int, n)
c = zeros(Int, n)
for i = 1:n
    a[i], b[i], c[i] = parse.(Int, split(readline()))
end
a # 1, 4, 7, 0
b # 2, 5, 8, 1
c # 3, 6, 9, 2

1 行に複数,入力の都度,処理

n = 5
for i = 1:n
    a, b, c = parse.(Int, split(readline()))
    # a, b, c を使用する処理
end

2 次元配列(行列)へ入力

n x m 行列に読み込む。入力時には m の情報は使わない。1 行当たりの入力個数に注意。

入力:
"1 2 3 4 5"
"6 7 8 9 10"
"11 12 13 14 15"

n, m = 3, 5
A = zeros(Int, n, m)
for i =1:n
    A[i, :] = parse.(Int, split(readline()))
end
A
# 3×5 Array{Int64,2}:
#   1   2   3   4   5
#   6   7   8   9  10
#  11  12  13  14  15

特定の文字列を読むまで繰り返し入力

入力:
"1"
"2"
"3"
"finish"

s = 0
while true
    e = readline()
    e == "finish" && break
    s += parse(Int, e)
end
s # 6

空行が入力されるまで繰り返し入力

入力:
"1"
"2"
"3"
"finish"

s = 0
while true
    try
        e = readline()
        s += parse(Int, e)
    catch
        break
    end
end
s # 6

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia で GLM

2021年02月18日 | ブログラミング

############# R

> library(COUNT)
> example(affairs)

> glmaffp <- glm(naffairs ~ kids + yrsmarr2 + yrsmarr3 + yrsmarr4 + yrsmarr5,
+                family = poisson, data = affairs)
> summary(glmaffp)

Call:
glm(formula = naffairs ~ kids + yrsmarr2 + yrsmarr3 + yrsmarr4 + 
    yrsmarr5, family = poisson, data = affairs)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9668  -1.9364  -1.5412  -0.9274   7.0799  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.34038    0.09182   3.707  0.00021
kids         0.28809    0.09371   3.074  0.00211
yrsmarr2    -1.18431    0.17058  -6.943 3.84e-12
yrsmarr3    -0.45650    0.10536  -4.333 1.47e-05
yrsmarr4    -0.11823    0.09896  -1.195  0.23220
yrsmarr5     0.03119    0.09912   0.315  0.75303

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2925.5  on 600  degrees of freedom
Residual deviance: 2797.0  on 595  degrees of freedom
AIC: 3303

Number of Fisher Scoring iterations: 7

> exp(coef(glmaffp))
(Intercept)        kids    yrsmarr2    yrsmarr3    yrsmarr4    yrsmarr5 
  1.4054793   1.3338755   0.3059569   0.6334955   0.8884915   1.0316802 

############# Julia

using GLM, RDatasets
affairs = dataset("COUNT", "affairs");

glmaffp = glm(@formula(NAffairs ~ Kids + YrsMarr2 + YrsMarr3 + YrsMarr4 + YrsMarr5),
  affairs, Poisson())

StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Array{Float64,1},Poisson{Float64},LogLink},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

NAffairs ~ 1 + Kids + YrsMarr2 + YrsMarr3 + YrsMarr4 + YrsMarr5

Coefficients:
───────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      z  Pr(>|z|)  Lower 95%   Upper 95%
───────────────────────────────────────────────────────────────────────────
(Intercept)   0.340379    0.0918004   3.71    0.0002   0.160453   0.520304
Kids          0.288088    0.0936936   3.07    0.0021   0.104452   0.471725
YrsMarr2     -1.18431     0.170531   -6.94    <1e-11  -1.51854   -0.850077
YrsMarr3     -0.456502    0.10535    -4.33    <1e-4   -0.662984  -0.250021
YrsMarr4     -0.11823     0.0989574  -1.19    0.2322  -0.312183   0.0757227
YrsMarr5      0.0311887   0.0991201   0.31    0.7530  -0.163083   0.225461
───────────────────────────────────────────────────────────────────────────

exp.(coef(glmaffp))

6-element Array{Float64,1}:
 1.4054795343529018
 1.333875284127698
 0.3059569619008904
 0.6334954669694786
 0.8884915158052847
 1.0316802145278257

############# R

> require(MASS)
> glmaffnb <- glm.nb(naffairs ~ kids + yrsmarr2 + yrsmarr3 + yrsmarr4 + yrsmarr5,
+                    data=affairs)
> summary(glmaffnb)

Call:
glm.nb(formula = naffairs ~ kids + yrsmarr2 + yrsmarr3 + yrsmarr4 + 
    yrsmarr5, data = affairs, init.theta = 0.1188516427, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8253  -0.8155  -0.7611  -0.6002   1.9331  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.33180    0.31247   1.062  0.28830
kids         0.27309    0.31104   0.878  0.37995
yrsmarr2    -1.19445    0.42962  -2.780  0.00543
yrsmarr3    -0.38936    0.35449  -1.098  0.27205
yrsmarr4    -0.08137    0.38166  -0.213  0.83116
yrsmarr5     0.07220    0.40464   0.178  0.85838

(Dispersion parameter for Negative Binomial(0.1189) family taken to be 1)

    Null deviance: 344.44  on 600  degrees of freedom
Residual deviance: 332.40  on 595  degrees of freedom
AIC: 1504.6

Number of Fisher Scoring iterations: 1


              Theta:  0.1189 
          Std. Err.:  0.0127 

 2 x log-likelihood:  -1490.6260 
> exp(coef(glmaffnb))
(Intercept)        kids    yrsmarr2    yrsmarr3    yrsmarr4    yrsmarr5 
  1.3934701   1.3140193   0.3028717   0.6774934   0.9218498   1.0748732 

############# Julia

glmaffnb = glm(@formula(NAffairs ~ Kids + YrsMarr2 + YrsMarr3 + YrsMarr4 + YrsMarr5),
  affairs, NegativeBinomial(0.11), LogLink())

StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Array{Float64,1},NegativeBinomial{Float64},LogLink},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

NAffairs ~ 1 + Kids + YrsMarr2 + YrsMarr3 + YrsMarr4 + YrsMarr5

Coefficients:
──────────────────────────────────────────────────────────────────────────
                  Coef.  Std. Error      z  Pr(>|z|)  Lower 95%  Upper 95%
──────────────────────────────────────────────────────────────────────────
(Intercept)   0.331862     0.246782   1.34    0.1787  -0.151823   0.815546
Kids          0.273032     0.24558    1.11    0.2662  -0.208296   0.75436
YrsMarr2     -1.19469      0.33843   -3.53    0.0004  -1.858     -0.531382
YrsMarr3     -0.389315     0.280041  -1.39    0.1645  -0.938186   0.159557
YrsMarr4     -0.0813931    0.301714  -0.27    0.7873  -0.672741   0.509955
YrsMarr5      0.0721675    0.319973   0.23    0.8216  -0.554969   0.699304
──────────────────────────────────────────────────────────────────────────

exp.(coef(glmaffnb))

6-element Array{Float64,1}:
 1.3935600110230184
 1.3139428682448457
 0.3027969691945509
 0.6775210913867913
 0.9218312041654065
 1.0748353924095455

deviance(glmaffnb) # 314.9101838376581

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

R 4.0.4 が出たけど(問題解消)

2021年02月16日 | ブログラミング
2021/02/23 に,

R version 4.0.4 Patched (2021-02-17 r80031) で

解消されたことが確認された

=======================================

日本語がちゃんと表示できないのは,beta 版のまま直っていないなあ

> a = "今日は"
> a
[1] "\u4eca\u65e5は"

> b = "日本語表示カタカナひらがな"
> b
[1] "\u65e5\u672c\u8a9e\u8868\u793aカタカナひらがな"

Mac 版だけの症状かと思ったが,Windows 版でも同じだった。漢字だけ表示できない。

そのあと,色々やっているが,cat による出力では漢字もちゃんと表示できるようだ。

> c = "R は、自由なソフトウェア"
> cat(c)
R は、自由なソフトウェア

sprintf も漢字の出力はだめだが(出力の段階で結局は print を使うから),cat と組み合わせるとオーケーだ。

> sprintf("%s", "R は多くの貢献者による共同プロジェクトです")
[1] "R は\u591aくの\u8ca2\u732e\u8005による\u5171\u540cプロジェクトです"

> cat(sprintf("%s", "R は多くの貢献者による共同プロジェクトです"))
R は多くの貢献者による共同プロジェクトです

> d = sprintf("%s version %d.%d.%d は日本語だめ?", "R", 4, 0, 4)
> cat(d)
R version 4.0.4 は日本語だめ?

当分の間は cat, sprintf でお茶を濁しておこう。

コメント (2)
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

プログラムの書き方(Python)

2021年02月16日 | ブログラミング
if bmi < 18.5:
  msg = "underweight"
elif 18.5 <= bmi < 25:
  msg = "normal weight"
elif 25 <= bmi < 30:
  msg = "over weight"
elif 30 <= bmi < 35:
  msg ="obese"
else:
  msg = "clinically obese. Go to the hospital immediately"

elif の条件判定が冗長

if bmi < 18.5:
  msg = "underweight"
elif bmi < 25:
  msg = "normal weight"
elif bmi < 30:
  msg = "over weight"
elif bmi < 35:
  msg ="obese"
else:
  msg = "clinically obese. Go to the hospital immediately"

 

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Python のリスト内包表記について

2021年02月15日 | ブログラミング

リスト内包表記って,速くもなんともないようです。むしろ遅い。
「1行で簡単に書ける」なんて言葉にだまされて,計算時間が掛かってしまっては本末転倒です。

一重ループ

def a1(n):
  s = 0
  for i in range(1, n+1):
      s += 1/i
  return s

def b1(n):
  return sum([ 1 / i for i in range(1, n+1)])

Python  n=10**6  n=10**7  n=10**8  n=10**9
    a1  0.080    0.840     8.395     83.437
    b1  0.104    1.030    10.314   1103.242

n=10**9 のときの b1 は,ガーベージ・コレクションでもが働いているのか?,計算時間が不定

二重ループ

def a2(n):
  s = 0
  for i in range(1, n+1):
    for j in range(1, n+1):
      s += i/j
  return s

def b2(n):
  return sum([ i / j for i in range(1, n+1) for j in range(1, n+1)])

Python  n=10**3  n=10**4  n=10**5  n=10**6
  a2   0.081     8.217 
  b2   0.099    10.172   

n が 10**5 以上では,これまたガーベージ・コレクションの影響か,終わりそうにないので中止

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

numba が速いってさ(Julia もほぼ同じ)

2021年02月14日 | ブログラミング

処理系ごとの計算所要時間(n は 外側の for ループの 10^n の n)

         n = 0    n = 1  n = 2  n = 3
numba    0.783    1.62   16.9   167
Julia    0.194    1.833  17.531 174.838
Fortran  0.866    8.638  97.865 ---
Java     2.298   25.903 249.377 ---
R       52.671  555.271 ---
Python 229.654 2794.784 ---

  • numba と Julia はいい勝負
  • Julia は Fortran より速い
  • Python おそっ!

プログラム

numba

from numba import jit

@jit
def hoge(n):
    count = 0
    i = 0
    for k0 in range(10**n):
        for k in range(10**9):
            i += 1
            if i >= 1000:
                count += 1
    return count

julia

function hoge(n)
    count = 0
    i = 0
    for k0 = 1:10^n
    for k = 1:10^9
        i += 1
        if i >= 1000
            count += 1
        end
    end
    end
    count
end

Fortran

program main
    implicit none
    integer(8) n, i, k, k0, count
    n = 0
    count = 0
    i = 0
    do k0 = 1, 10**n
        do k = 1, 10**9
            i = i+1
            if (i >= 1000) then
                count = count + 1
            end if
        end do
    end do
    write(*,*) n, count
end program main

Java

package hoge;

public class hogehoge {

    public static void main(String[] args) {
        long n, i, k, k0, count;
        long startTime = System.currentTimeMillis();
        n = 2;
        count = 0;
        i = 0;
        for (k0 = 1; k0 <= Math.pow(10, n); k0++) {
            for (k = 1; k <= Math.pow(10, 9); k++) {
                i += 1;
                if (i >= 1000) {
                    count += 1;
                }
            }
        }
        System.out.println(n+", "+count);
        long endTime = System.currentTimeMillis();
        System.out.println("処理時間:" + (endTime - startTime) + " ms");
    }
}

R

hoge = function(n) {
    count = 0
    i = 0
    for (k0 in 1:10^n) {
        for (k in 1:10^9) {
            i = i+1
            if (i >= 1000) {
                count = count + 1
            }
        }
    }
    count
}

Python

def hoge(n):
    count = 0
    i = 0
    for k0 in range(10**n):
        for k in range(10**9):
            i += 1
            if i >= 1000:
                count += 1
    return count

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

リスト内包表記は,リストしか使えない Python のひがみ

2021年02月13日 | ブログラミング

別に Python がリストしか使えない訳ではない

以下の例も,numpy ndarray などを使うとか,しかるべき関数を使えばリスト内包表記など使わなくてもできる話である。

しかしなぜか,「リスト内包表記はすごいんだぜ!」的な上から目線の記事が多いので,そうでもないよねというのを書いてみる。

# full_name = "Yamada Taro"
# characters = [char for char in full_name]

fullname = "Yamada Taro"
characters = split(fullname, "")

# Matrix = [[7, 1, 5],
#           [22, 99, 0],
#           [33, 2, 57]]
# row_min = [min(row) for row in Matrix]

M = [7 1 5; 22 90 0; 33 2 57]

rowmin = minimum(M, dims=2)

# People = ["Taro", "Jiro", 'Wan', "Elepahnt", "Tomas", "ken", "sae"]
# target = [name for name in People if len(name) < 4 and name.islower()]

people = ["Taro", "Jiro", "Wan", "Elepahnt", "Tomas", "ken", "sae"]
select(s) = (length(s) < 4) & all("a" .<= split(s, "") .<= "z")
target = people[select.(people)]

# Genius = ["taro", "jiro", "yamada", "ito"]
# target = [name.capitalize() for name in Genius]

genius = ["taro", "jiro", "yamada", "ito"]
target = titlecase.(genius)

# People = ["Taro", "Hiroyuki", "ken", "youu"]
# answer = [name if name.startswith('y') else 'Not Genius' for name in People]

people = ["Taro", "Hiroyuki", "ken", "youu"]
starts(s) = (s[1] == 'y') + 1
answer = ["not Genius", "youu"][starts.(people)]

# Genius = ["taro", "jiro", "yamada"]
# L = []
# for name in Genius:
#     for char in name:
#         L.append(char)

genius = ["taro", "jiro", "yamada"]
split(join(genius), "")

# People = ["Kenny", "Tomas", "tom", "young"]
# L = map(lambda a: a.lower(), People)

people = ["Kenny", "Tomas", "tom", "young"]
select(s) = all("a" .<= split(s, "") .<= "z")
L = lowercase.(people)

# L1 = filter(lambda a: len(a) < 4, People)
# L2 = [a for a in People if len(a) < 4]

people = ["Kenny", "Tomas", "tom", "young"]
L2 = people[length.(people) .< 4]

リスト内包表記とジェネレータのメモリー使用量

# large_list = [x for x in range(1_000_000)] # 8697440
# large_list_g = (x for x in range(1_000_000)) # 96

Base.summarysize(collect(1:1000000))
Base.summarysize(1:1000000)

# list_2d = [
#      [0, 0, 0],
#      [1, 1, 1],
#      [2, 2, 2],
# ]
# flat = [num ** 2 if num % 2 == 0 else num for row in list_2d for num in row]

# list_2d = [
#      [0, 0, 0],
#      [1, 1, 1],
#      [2, 2, 2],
# ]
# flat = []
# for row in list_2d:
#     for num in row:
#         flat.append(num ** 2 if num % 2 == 0 else num)

list2d = [0 0 0; 1 1 1; 2 2 2]
flat = vec(list2d' .^ 2)

いかがでしたでしょうか。リスト内包表記は記述法としてはわかりにくい(そもそも同じ処理を書いても,長くなりがち。長いプログラムはそれだけでよくない)。

対して,ベクトル処理できる関数で対処できる Julia は,関数名を見ただけでもわかる。R も Julia と同じ感じ(だからリスト内包表記などない)。Julia では,内包表記もできる(が,使う局面は少ない)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia と R で整数演算の違い

2021年02月13日 | ブログラミング

using RCall

R"x <- -5:8"      # -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8
x = -5:8          # -5:8
x = collect(-5:8) # -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8

R"x  + 1"  # -4 -3 -2 -1  0  1  2  3  4  5  6  7  8  9
  x .+ 1   # -4 -3 -2 -1  0  1  2  3  4  5  6  7  8  9

R"2*x + 3" # -7 -5 -3 -1  1  3  5  7  9 11 13 15 17 19
  2x .+ 3  # -7 -5 -3 -1  1  3  5  7  9 11 13 15 17 19

R"x %% 2"    #  1  0   1  0   1  0  1  0  1  0  1  0  1  0
  mod.(x, 2) #  1  0   1  0   1  0  1  0  1  0  1  0  1  0
  x .% 2     # -1  0  -1  0  -1  0  1  0  1  0  1  0  1  0
  rem.(x, 2) # -1  0  -1  0  -1  0  1  0  1  0  1  0  1  0

R"x %/% 5"              # -1  -1  -1  -1  -1  0  0  0  0  0  1  1  1  1
  x .÷  5               # -1   0   0   0   0  0  0  0  0  0  1  1  1  1
  div.(x, 5)            # -1   0   0   0   0  0  0  0  0  0  1  1  1  1
  fld.(x, 5)            # -1  -1  -1  -1  -1  0  0  0  0  0  1  1  1  1
  div.(x, 5, RoundDown) # -1  -1  -1  -1  -1  0  0  0  0  0  1  1  1  1
  div.(x, 5)            # -1   0   0   0   0  0  0  0  0  0  1  1  1  1
  div.(x, 5, RoundUp)   # -1   0   0   0   0  0  1  1  1  1  1  2  2  2
  round.(Int, x ./ 5)   # -1  -1  -1   0   0  0  0  0  1  1  1  1  1  2

R"x %% Inf" #  Inf  Inf  Inf  Inf  Inf   0   1   2   3   4   5   6   7   8
  x .% Inf  # -5.0 -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0

R"x  + Inf" # Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
  x .+ Inf  # Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf

R"x  * Inf" # -Inf -Inf -Inf -Inf -Inf NaN Inf Inf Inf Inf Inf Inf Inf Inf
  x .* Inf  # -Inf -Inf -Inf -Inf -Inf NaN Inf Inf Inf Inf Inf Inf Inf Inf

R"x  / Inf" #    0    0    0    0    0   0   0   0   0   0   0   0   0   0
  x ./ Inf  # -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村