Julia Array操作の基本をメモ。

Pythonでやっていたlistの操作にすっかり慣れてしまったけど、JuliaのArrayではどうやるか確かめる。
Julia0.3で試して動いたものを掲載。

配列の生成

a = [i for i in 1:10]     #1から10を要素に持つ配列生成
a = [1:10]                #同上
b = reshape([1:9], (3,3)) #3x3の2次元配列。
b = [1 2 3; 4 5 6; 7 8 9] #3x3の2次元配列。上とは行、列が逆(転置の関係)

要素数

length(a)

要素アクセス

a[1]
a[1:end]   #要素1から最後まで
a[2:3:8] #要素2から8を3個おきに

追加,挿入

push!(a, 13)         #要素と同じ型追加
append!(a, [1,2,3])  #末尾に[1,2,3]を追加
prepend!(a, [1,2,3]) #先頭に[1,2,3]を追加
insert!(a, 3, 15)    #3番目に15を挿入する。(要素数は増える)
vcat(a, [11,12,13])  #aに末尾に[11,12,13]を追加(aは非破壊になる)
hcat(b, [10,11,12])  #bに列を追加
vcat(b, [10,11,12]') #bに行を追加('は転置)

要素数があらかじめわかっているなら、始めに配列を確保して代入したほうが速いっぽい。

置換

a[1] = 15
a[1:4] = [4,4,4,4]
#a[2:5]を削除して[0,0,0]を挿入。戻り値は削除した元のaの[2:5]の値。
splice!(a, 2:5, [0,0,0]) 

削除

shift!(a)  #最初の要素を削除
pop!(a)    #最後の要素を削除

検索

3 in a
findfirst(a, 4) #はじめに4が出る要素番号、なければ0を返す
findin(a, 4)    #4の存在する要素数をArrayに入れて返す)

Juliaで固定配列

JuliaでArrayを使えばベクトルや行列計算できるけど、2とか3次元の点を表すのに可変長の配列を使うのはどうも無駄が多い気がする。自前で

immutable Vec
    x::Float64
    y::Float64
    z::Float64
end

とかやったけど、演算子とかを定義するのが面倒。

そんな問題を解決してくれるのがImuutableArrays。
Julia0.3にはすでに採用されているので、

Using ImuutableArrays
v3 = Vector3{Float64}(1, 2, 3)
v2 = Vector2{Float64}(4, 5)

と書ける。デフォルトでは読み込んだ時に4次元まで作ってくれて、演算子の定義とかもやってくれている。(これがメタプログラミングの力か…)
しかも上で自分で定義したのより速い…

行列は

typealias Mat3d Matrix3x3{Float64}

とすれば使える。

Juliaで簡単な回帰

Juliaの練習がてら簡単な線形回帰を行った。

以下のことが分かる練習にはちょうど良いシンプルなコードだ。

• 乱数をどう発生させるか
• 計画行列をどう作るか
• Gadflyを使って2つのプロットを重ねるにはどうするか
• Gadflyを使って軸の設定をどうするか

これらがわかれば自作する当面の計算にはこまらない気がする…

linreg()はバイアス項を勝手に計算してくれるみたいなので、入れなくてよいみたいだ。

# #Regression Example

using Gadfly
using Distributions

N = 15
xdata = (rand(N) - 0.5) * 2 * pi
σ² = 0.5
function noisy(x, r=rand(Normal(0, σ²))
    return sin(x) + r
end

#X results in Nx9 array
t = noisy(xdata);
X = Float64[x^i for x=xdata, i=[1:9]] #bias is not needed for linreg(). it will interpolate automatically

#regression by linreg
coeffs = linreg(X, t)

#plotting
f = x -> coeffs[1] + sum([coeffs[i]*x^(i-1) for i = [2:10]])

x_reg = linspace(-pi, pi, 100)
y_reg = [f(xi) for xi=x_reg]

plot(layer(x=x_reg, y=y_reg, Geom.line), layer(x=xdata, y=t, Geom.point),
    Coord.Cartesian(xmin=-pi, xmax=pi, ymin=-1.5, ymax=1.5))