書いてあること

自作の回帰モデルを，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から最低限使えるようになりました．