はじめに (対象読者・この記事でわかること)

この記事は、JavaScriptでデータの回帰分析を行いたい開発者、特にPythonのscipy.optimize.curve_fitのような機能をJavaScriptで実現したい方を対象としています。Webブラウザ上やNode.js環境でデータ分析を行いたい場合に役立ちます。

この記事を読むことで、JavaScriptで曲線フィットを実装する方法、既存のライブラリの使い方、シグモイド関数に特化した回帰アルゴリズムの実装方法がわかります。また、実際のコード例を交えて具体的な実装手順を学ぶことができます。

前提知識

この記事を読み進める上で、以下の知識があるとスムーズです。 - JavaScriptの基本的な知識 - 高校数学の関数と微分の基礎知識 - 簡単な統計学の知識（回帰分析の概念）

JavaScriptでの曲線フィットの概要と背景

Pythonの科学計算ライブラリであるSciPyには、curve_fitという強力な関数があり、任意の関数形式に基づいたデータの非線形回帰分析が簡単に行えます。この機能は、実験データの解析や機械学習の前処理など、様々な場面で活用されています。

しかし、JavaScript環境ではこのような高機能な曲線フィットライブラリが標準では提供されていません。Webブラウザ上でデータ分析を行いたい場合や、Node.js環境でデータ処理を行いたい場合に、JavaScriptで同様の機能を実装する必要があります。

JavaScriptで曲線フィットを実装する主な方法として、以下の3つのアプローチが考えられます。

  1. 既存のJavaScriptライブラリを利用する方法
  2. PythonのライブラリをWebAssemblyで利用する方法
  3. 自前でアルゴリズムを実装する方法

それぞれの方法にはメリット・デメリットがあり、用途や要件に応じて適切な方法を選ぶ必要があります。この記事では、これらの方法を具体的に解説します。

JavaScriptで曲線フィットを実装する具体的な方法

方法1: 既存のJavaScriptライブラリを利用する

JavaScriptには、数値計算や統計処理を行うためのライブラリがいくつか存在します。中でも、曲線フィット機能を提供しているものをいくつか紹介します。

ml-regression

ml-regressionは、機械学習の回帰アルゴリズムを提供するライブラリです。線形回帰、多項式回帰、ロジスティック回帰などがサポートされています。

インストール:

Bash

npm install ml-regression

使用例:

Javascript

const LinearRegression = require('ml-regression-linear');
const PolynomialRegression = require('ml-regression-polynomial');

// 線形回帰の例
const x = [1, 2, 3, 4, 5];
const y = [2, 4, 5, 4, 5];

const linearRegression = new LinearRegression(x, y);
console.log(`予測値(x=6): ${linearRegression.predict(6)}`);

// 多項式回帰の例
const polynomialRegression = new PolynomialRegression(x, y, 2);
console.log(`予測値(x=6): ${polynomialRegression.predict(6)}`);

regression

regressionは、様々な回帰アルゴリズムを提供するライブラリです。線形回帰、多項式回帰、指数関数回帰などがサポートされています。

インストール:

Bash

npm install regression

使用例:

Javascript

const regression = require('regression');

const data = [
    [1, 2],
    [2, 4],
    [3, 5],
    [4, 4],
    [5, 5]
];

// 多項式回帰の例
const result = regression.polynomial(data, { order: 2 });
console.log(result.string);
console.log(`予測値(x=6): ${result.predict(6)[1]}`);

方法2: WebAssemblyでPythonライブラリを利用する

JavaScriptの機能だけでは限界がある場合、WebAssemblyを利用してPythonのライブラリをブラウザ上で実行する方法があります。特に、SciPyのcurve_fitのような高度な機能が必要な場合に有効です。

Pyodideの利用

Pyodideは、Pythonの科学計算ライブラリをWebAssemblyで実行するためのプロジェクトです。ブラウザ上でNumPyやSciPyが利用できます。

使用例:

Html

<!DOCTYPE html>
<html>
<head>
    <title>Pyodideでcurve_fit</title>
</head>
<body>
    <script src="https://cdn.jsdelivr.net/pyodide/v0.21.3/full/pyodide.js"></script>
    <script>
        async function main() {
            // Pyodideの初期化
            await loadPyodide();

            // Pythonコードの実行
            pyodide.runPython(`
                import numpy as np
                from scipy.optimize import curve_fit
                import js

                # データの準備
                x = np.array([1, 2, 3, 4, 5])
                y = np.array([2, 4, 5, 4, 5])

                # シグモイド関数の定義
                def sigmoid(x, a, b, c):
                    return a / (1 + np.exp(-b * (x - c)))

                # フィッティング
                popt, pcov = curve_fit(sigmoid, x, y)
                print(f"フィッティング結果: a={popt[0]}, b={popt[1]}, c={popt[2]}")

                # JavaScriptの配列に結果を渡す
                result = [float(popt[0]), float(popt[1]), float(popt[2])]
                js.results = result
            `);

            // 結果の表示
            console.log("フィッティング結果:", pyodide.globals.get('results').toJs());
        }

        main();
    </script>
</body>
</html>

方法3: 自前でアルゴリズムを実装する

既存のライブラリやWebAssemblyが利用できない場合、または特定の要件を満たすカスタムアルゴリズムが必要な場合は、自前で実装する必要があります。ここでは、シグモイド関数に特化した回帰アルゴリズムの実装例を紹介します。

シグモイド関数のフィッティング

シグモイド関数は、以下の式で表されます。

f(x) = L / (1 + e^(-k(x - x0)))

ここで、Lは関数の最大値、kは傾き、x0はxの中心位置を表します。

この関数のフィッティングを実装するには、非線形最小二乗法を用いる必要があります。ここでは、レーベンベルグ・マルカート法（Levenberg-Marquardt algorithm）を実装します。

Javascript

class SigmoidRegression {
    constructor() {
        this.params = [1, 1, 1]; // 初期パラメータ [L, k, x0]
    }

    // シグモイド関数
    sigmoid(x, L, k, x0) {
        return L / (1 + Math.exp(-k * (x - x0)));
    }

    // 残差の計算
    residuals(params, x, y) {
        return y.map((yi, i) => yi - this.sigmoid(x[i], params[0], params[1], params[2]));
    }

    // ヤコビアンの計算
    jacobian(params, x) {
        return x.map(xi => {
            const exp_term = Math.exp(-params[1] * (xi - params[2]));
            const denom = 1 + exp_term;

            return [
                1 / denom, // dL/dparam
                params[0] * exp_term * (xi - params[2]) / (denom * denom), // dk/dparam
                params[0] * exp_term * params[1] / (denom * denom) // dx0/dparam
            ];
        });
    }

    // レーベンベルグ・マルカート法の実装
    fit(x, y, maxIterations = 100, tolerance = 1e-6) {
        let params = [...this.params];
        let lambda = 0.001; // ダンピングパラメータ

        for (let iter = 0; iter < maxIterations; iter++) {
            const residuals = this.residuals(params, x, y);
            const jacobian = this.jacobian(params, x);

            // J^T * J + lambda * I の計算
            const jtj = this.multiplyJTJ(jacobian);
            const jtjLambda = jtj.map((row, i) => 
                row.map((val, j) => val + (i === j ? lambda : 0))
            );

            // J^T * r の計算
            const jtr = this.multiplyJTr(jacobian, residuals);

            // 連立方程式の解を求める (JTJ * delta = JTr)
            const delta = this.solveLinearSystem(jtjLambda, jtr);

            // パラメータの更新
            const newParams = params.map((val, i) => val + delta[i]);

            // 残差二乗和の評価
            const newResiduals = this.residuals(newParams, x, y);
            const currentError = this.sumOfSquares(residuals);
            const newError = this.sumOfSquares(newResiduals);

            if (newError < currentError) {
                // 改善した場合
                params = newParams;
                lambda *= 0.1; // lambdaを減少
            } else {
                // 悪化した場合
                lambda *= 10; // lambdaを増加
                continue;
            }

            // 収束判定
            const maxDelta = Math.max(...delta.map(Math.abs));
            if (maxDelta < tolerance) {
                break;
            }
        }

        this.params = params;
        return this.params;
    }

    // 行列演算のヘルパー関数
    multiplyJTJ(jacobian) {
        const rows = jacobian.length;
        const cols = jacobian[0].length;
        const result = Array(cols).fill().map(() => Array(cols).fill(0));

        for (let i = 0; i < cols; i++) {
            for (let j = 0; j < cols; j++) {
                for (let k = 0; k < rows; k++) {
                    result[i][j] += jacobian[k][i] * jacobian[k][j];
                }
            }
        }

        return result;
    }

    multiplyJTr(jacobian, residuals) {
        const cols = jacobian[0].length;
        const result = Array(cols).fill(0);

        for (let i = 0; i < cols; i++) {
            for (let k = 0; k < jacobian.length; k++) {
                result[i] += jacobian[k][i] * residuals[k];
            }
        }

        return result;
    }

    solveLinearSystem(a, b) {
        const n = a.length;

        // 前進消去
        for (let i = 0; i < n; i++) {
            // ピボット選択
            let maxRow = i;
            for (let k = i + 1; k < n; k++) {
                if (Math.abs(a[k][i]) > Math.abs(a[maxRow][i])) {
                    maxRow = k;
                }
            }
            [a[i], a[maxRow]] = [a[maxRow], a[i]];
            [b[i], b[maxRow]] = [b[maxRow], b[i]];

            // ピボット要素で割る
            const pivot = a[i][i];
            for (let j = i; j < n; j++) {
                a[i][j] /= pivot;
            }
            b[i] /= pivot;

            // その他の行から引く
            for (let k = i + 1; k < n; k++) {
                const factor = a[k][i];
                for (let j = i; j < n; j++) {
                    a[k][j] -= factor * a[i][j];
                }
                b[k] -= factor * b[i];
            }
        }

        // 後退代入
        const x = new Array(n);
        for (let i = n - 1; i >= 0; i--) {
            x[i] = b[i];
            for (let j = i + 1; j < n; j++) {
                x[i] -= a[i][j] * x[j];
            }
        }

        return x;
    }

    sumOfSquares(arr) {
        return arr.reduce((sum, val) => sum + val * val, 0);
    }

    predict(x) {
        return this.sigmoid(x, this.params[0], this.params[1], this.params[2]);
    }
}

// 使用例
const x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
const y = [3, 5, 8, 12, 18, 25, 30, 33, 35, 36];

const regression = new SigmoidRegression();
const params = regression.fit(x, y);

console.log("フィッティング結果:");
console.log(`L = ${params[0]}`);
console.log(`k = ${params[1]}`);
console.log(`x0 = ${params[2]}`);

// 予測値の計算
const newX = [11, 12, 13];
const predictions = newX.map(xi => regression.predict(xi));
console.log("予測値:", predictions);

ハマった点やエラー解決

JavaScriptで曲線フィットを実装する際に遭遇する可能性のある問題とその解決策を以下に示します。

問題1: 収束しない

レーベンベルグ・マルカート法が収束しない場合、初期パラメータの設定が不適切な可能性があります。

解決策: - データの特性に合わせて初期パラメータを設定する - ダンピングパラメータlambdaの初期値を調整する - 最大反復回数を増やす

問題2: 数値的不安定性

ヤコビアンの計算時に数値的不安定が発生することがあります。

解決策: - ヤコビアンの計算に解析微分を使用する（可能であれば） - スケーリングを行う（入力データやパラメータを正規化する） - 数値微分のステップサイズhを調整する

問題3: 計算コストが高い

大量のデータや複雑な関数を扱う場合、計算コストが高くなることがあります。

解決策: - Web Workerを利用して計算を非同期化する - 計算を並列化する（WebAssemblyやSIMDを利用） - アルゴリズムを最適化する（例: より効率的な線形ソルバーを使用）

まとめ

本記事では、JavaScriptで曲線フィットを実装する方法について解説しました。

この記事を通して、JavaScript環境でデータの回帰分析を行うための様々な方法が理解できたことと思います。要件や環境に応じて適切な方法を選択し、効果的なデータ分析を実現してください。

今後は、より高度な回帰アルゴリズムの実装や、パフォーマンス最適化のテクニックについても記事にする予定です。

参考資料