書いてあること
自作の回帰モデルを,julia の機械学習フレームワークであるMLJ.jl
から"最低限"使えるようにする方法が書いてあります.
ただし,公式パッケージとして使えるようにするのではなく,あくまでも自分の計算機の中だけで使えるようにします.
MLJ.jl
そのものについてはほとんど書いていないので,知りたい人は公式ドキュメント(https://alan-turing-institute.github.io/MLJ.jl/dev/ )を眺めたり,ドキュメントの Getting Started で遊んでみてください.
環境
julia> versioninfo()
Julia Version 1.7.0
Commit 3bf9d17731 (2021-11-30 12:12 UTC)
Platform Info:
OS: macOS (x86_64-apple-darwin19.5.0)
CPU: Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, icelake-client)
Environment:
JULIA_NUM_THREADS = 4
準備
まず,追加する自作モデルを用意します.
適当な作業用のディレクトリで REPL に入って,
(@v1.7) pkg> generate MyMLJModel
Generating project MyMLJModel:
MyMLJModel/Project.toml
MyMLJModel/src/MyMLJModel.jl
shell> cd MyMLJModel
として説明用のモジュールを作って,そのディレクトリに移動しておきます.
./src/MyMLJModel.jl
の内容を以下のように編集します.
これはテキトーに作った説明用の自作モデルなので,モデルとしては非常にアレです.
このモデルは一応線形回帰になってますが,最小二乗法的なことはせずに,a
回乱数を振って一番良かったものを重みとします.
パラメータを指定してデータを突っ込むと何か出てくるブラックボックスとでも思ってください.
module MyMLJModel
using Random
export MyModel, myfit!, mypredict
mutable struct MyModel
a::Int
weights::Vector{Float64}
end
MyModel(a) = MyModel(a, Float64[])
function myfit!(model, X, y)
n_vars = size(X)[1]
best_error = Inf
for i = 1:model.a
weights = rand(n_vars)
y_hat = [sum(weights[i] * X[i, j] for i = 1:n_vars) for j = 1:size(X)[2]]
tmp_error = rmse(y_hat, y)
if tmp_error < best_error
model.weights = weights
best_error = tmp_error
end
end
end
mypredict(model, X) = [sum(model.weights[i] * X[i, j] for i = 1:size(X)[1]) for j = 1:size(X)[2]]
rmse(a, b) = sqrt(sum(x -> x^2, a[i] - b[i] for i in eachindex(a)) / length(a))
end # module
とりあえずこのモデルを動かしてみます.
script
というディレクトリを作って,script/example1.jl
でモデルを動かしてみます.
using MyMLJModel
f(x) = 0.2x[1] + 0.4x[2] + 0.6x[3]
X = rand(3, 10)
y = [f(X[:, i]) for i = 1:10]
mymodel = MyModel(100)
myfit!(mymodel, X, y)
yhat = mypredict(mymodel, X)
@show y
@show yhat
テキトーにデータX
, y
を作ってmyfit!
で学習して,mypredict
で予測します.
ここで,X
のサイズが (3, 10)
になっていることに注意してください.X
の行が特徴量に対応しています.つまり,1列(X[:, i]
)が1つのインスタンスに対応しています.
(@v1.7) pkg> activate .
Activating project at `~/hoge/MyMLJModel`
(MyMLJModel) pkg> add Random
Resolving package versions...
Updating `~/hoge/MyMLJModel/Project.toml`
[9a3f8284] + Random
Updating `~/hoge/MyMJLModel/Manifest.toml`
[9a3f8284] + Random
[ea8e919c] + SHA
[9e88b42a] + Serialization
julia> using Revise
julia> include("script/example1.jl")
[ Info: Precompiling MyMLJModel [68c9983a-15f9-41e9-9249-e301e2c6b7d7]
┌ Warning: Package MyMLJModel does not have Random in its dependencies:
│ - If you have MyMLJModel checked out for development and have
│ added Random as a dependency but haven't updated your primary
│ environment's manifest file, try `Pkg.resolve()`.
│ - Otherwise you may need to report an issue with MyMLJModel
└ Loading Random into MyMLJModel from project dependency, future warnings for MyMLJModel are suppressed.
y = [0.5601453907761953, 0.6879910740562954, 0.3530135116286392, 0.5567237922697389, 0.5370888314572486, 1.1237237581585515, 0.8041827332557248, 0.2703739272387356, 0.1411133336353244, 0.7799328076341898]
yhat = [0.49971973364491074, 0.7010051972588169, 0.3424426329476714, 0.42568976463547314, 0.5650044702311595, 1.025868358799277, 0.7260945619574919, 0.25485599443926044, 0.12670020159703554, 0.7674023144526892]
10-element Vector{Float64}:
0.49971973364491074
0.7010051972588169
0.3424426329476714
0.42568976463547314
0.5650044702311595
1.025868358799277
0.7260945619574919
0.25485599443926044
0.12670020159703554
0.7674023144526892
warning が出てますが,ここでは問題ないので一旦無視して,@show
で表示されたy
とyhat
を確認してみてください.yhat
がそれっぽい値になってると思います.
ここまでで,説明用の自作のモデルの準備は終了です.
MLJModelInterface.jl
を使って自作モデルを MLJ から使えるようにする
MLJ から自作のモデルを利用できるようにする機能はMLJModelInterface.jl
で提供されています.
以下の3つを実装すると,とりあえず自分のモデルが MLJ から使えるようになります.
- MLJ で使用するモデル型
fit
predict
上から順番に実装していきます.
モデル型を定義する
MLJ におけるモデルとは,学習に使用するハイパーパラメータを格納するものです.
ここで使っている自作モデルのハイパーパラメータは繰り返し回数a
のみなので,MLJのモデルはa
のみを持ちます.学習することで得られるパラメータ(ここでは重み)は含みません.
@mlj_model
というマクロを使ってモデルを定義します.
上の方ではMyModel
という名前の構造体が定義されていましたが,ここではMyRegressor
として新しくMLJ用のモデルを定義します.
import MLJBase
import MLJModelInterface
const MMI = MLJModelInterface
export MyRegressor
MMI.@mlj_model mutable struct MyRegressor <: MMI.Deterministic
a::Int = 100::(_ > 1)
end
パラメータには初期が設定されている必要がります.a::Int = 100
のようにして初期値を定義します.また,オプションとしてパラメータの範囲も定義することができ,::(_ > 1)
と追記することで,a > 1
という制約を設けられます.
@mlj_model
で生成されるコンストラクタの引数はすべてキーワード引数になります.
あと,外部から使うためにexport MyRegressor
も追記しておいてください.
MLJModelInterface
は長いので,本稿のプログラム中ではMMI
と書けるようにしておきます.
fit
を実装する
MMI.fit
を実装します. まず実装例を示します.
function MMI.fit(model::MyRegressor, verbosity, X, y, w = nothing)
Xmat = MLJBase.matrix(X, transpose = true)
mymodel = MyModel(model.a)
myfit!(mymodel, Xmat, y)
fitresult = mymodel
report = nothing
cache = nothing
return (fitresult, cache, report)
end
最初の引数model
には自分のモデルの型を書いておきます.
X
は入力で,行列やDataFrame
などが入ってきます.ここで注意点があるのですが,MLJ ではX
の列が特徴量に対応しています.つまり,1行が1つのインスタンスに対応しています.自作のモデルでそうなっていない場合,X
を変換しておく必要があります.MLJBase.jl
というパッケージのmatrix()
を使うと,DataFrame
を行列に変換したり,転置をとったりできます.
MMI.fit
の内部では,今まで使っていた自作モデルのMyModel
を学習させます.自作モデルで使っていた関数等を使って学習させます.
そして,fitresult
,cache
, report
を定義します.
fitresult
にはpredict
に必要な情報を入れます.predict
では自作モデル側のmypredict
が呼びたくて,そのmypredict
は引数にmymodel
が必要なので,fitresult = mymodel
としておきます.
cache
はここでは使わないのでnothing
とします.不要でもnothing
として返す必要があります.
report
は何かユーザーが知りたい情報を返すのに使うNamedTuple
です.ここでは自作モデルの重みを返すことにします.
predict
を実装する
MMI.predict
を実装します.まず実装例を示します.
function MMI.predict(model::MyRegressor, fitresult, Xnew)
Xmat = MLJBase.matrix(Xnew, transpose = true)
return mypredict(fitresult, Xmat)
end
最初の引数model
には自分のモデルの型を書いておきます.
ここでもfit
と同様に,X
をmypredict
に渡すためにMLJBase.matrix
で変換しておきます.
あとはmypredict
を呼んでその結果を返すだけです.
自作モデルをMLJから使ってみる
これで自作モデルをMLJから使うための最低限の用意が整いました.
試しに,以下に示すscript/example2.jl
で動作を確認してみます.
using MyMLJModel
using MLJ
f(x) = 0.2x[1] + 0.4x[2] + 0.6x[3]
X = rand(10, 3) # n × p
y = [f(X[i, :]) for i = 1:10]
model = MyRegressor(a = 100)
mach = machine(model, X, y)
evaluate!(mach, measure = rms)
実行すると,以下のような結果が出てくると思います.
julia> include("script/example2.jl")
Evaluating over 6 folds: 100%[=========================] Time: 0:00:00
PerformanceEvaluation object with these fields:
measure, measurement, operation, per_fold,
per_observation, fitted_params_per_fold,
report_per_fold, train_test_pairs
Extract:
┌────────────────────────┬─────────────┬───────────┬─────────────────────────────────────────────────┐
│ measure │ measurement │ operation │ per_fold │
├────────────────────────┼─────────────┼───────────┼─────────────────────────────────────────────────┤
│ RootMeanSquaredError() │ 0.0691 │ predict │ [0.14, 0.0514, 0.0208, 0.0332, 0.00209, 0.0701] │
└────────────────────────┴─────────────┴───────────┴─────────────────────────────────────────────────┘
チューニングもしてみます.
model = MyRegressor(a = 100)
r = range(model, :a, lower = 10, upper = 1000)
tuned_model = TunedModel(model = model, resampling = CV(nfolds = 5), tuning = Grid(resolution = 10), range = r, measure = rms)
fit!(mach)
julia> mach.report.plotting.measurements
10-element Vector{Float64}:
0.02442430654777976
0.06935442889813258
0.01598692079500873
0.04161092341230241
0.0268491382201212
0.056923740156152744
0.03366051708222645
0.039598241639119326
0.03212408357200946
0.03651355084042113
julia> fitted_params(mach).best_model
MyRegressor(
a = 890)
julia> mach.report.best_history_entry
(model = MyRegressor,
measure = [RootMeanSquaredError()],
measurement = [0.01598692079500873],
per_fold = [[0.016621197285246875, 0.008638543600675262, 0.007427216316183303, 0.025226911542651647, 0.01534467137495609]],)
a
回試行できるのでa
が大きい方が良さそうですが,10から1000までしか動かしてませんし運ゲーなので,それっぽい結果は得られてると思います.
これで自作モデルをMLJから最低限使えるようになりました.