[PHP]log1p関数完全解説 – 高精度なlog(1+x)計算と金融・統計分野での実用的活用法

PHP

PHPで数値計算を行う際、小さな値に対してlog(1+x)を計算する必要がある場面があります。そんなときに威力を発揮するのがlog1p関数です。この記事では、log1p関数の特徴から実際のプログラムでの活用方法まで、サンプルコードを交えながら詳しく解説していきます。

log1p関数とは?

log1p関数は、log(1 + x)を高精度で計算するPHPの数学関数です。単純にlog(1 + x)と書くよりも、特にxが0に近い小さな値の場合に、数値的精度が大幅に向上します。金融計算、統計処理、機械学習などの分野で、精密な計算が要求される場面で重宝される関数です。

基本的な構文

float log1p(float $number)

パラメータ

  • $number: 計算する値(xの部分。-1より大きい値である必要があります)

返り値

  • log(1 + $number)の値(float型)
  • $numberが-1以下の場合は-INFまたはNANを返す

なぜlog1p関数が必要なのか?

精度の問題

<?php
// 小さな値での精度比較
$small_values = [1e-10, 1e-15, 1e-20];

echo "精度比較(小さな値での計算):\n";
foreach ($small_values as $x) {
    $log_method = log(1 + $x);
    $log1p_method = log1p($x);
    
    echo "x = {$x}:\n";
    echo "  log(1 + x)    = " . sprintf("%.20e", $log_method) . "\n";
    echo "  log1p(x)      = " . sprintf("%.20e", $log1p_method) . "\n";
    echo "  差            = " . sprintf("%.20e", abs($log_method - $log1p_method)) . "\n\n";
}

// 理論値との比較(テイラー展開: log(1+x) ≈ x - x²/2 + x³/3 - ...)
echo "理論値との比較:\n";
$x = 1e-10;
$theoretical = $x - ($x * $x) / 2; // 2次まで

echo "x = {$x}:\n";
echo "理論値(2次近似): " . sprintf("%.20e", $theoretical) . "\n";
echo "log1p(x)         : " . sprintf("%.20e", log1p($x)) . "\n";
echo "log(1+x)         : " . sprintf("%.20e", log(1 + $x)) . "\n";
?>

基本的な使用例

基本的な計算

<?php
// 基本的な使用例
$test_values = [0, 0.1, 0.5, 1, 2, -0.5, -0.9];

echo "log1p関数の基本的な使用例:\n";
foreach ($test_values as $x) {
    if ($x > -1) {
        $result = log1p($x);
        $comparison = log(1 + $x);
        
        echo "x = {$x}:\n";
        echo "  log1p({$x}) = {$result}\n";
        echo "  log(1+{$x}) = {$comparison}\n";
        echo "  差 = " . abs($result - $comparison) . "\n\n";
    } else {
        echo "x = {$x}: 定義域外(x > -1が必要)\n\n";
    }
}

// 特別な値
echo "特別な値での動作:\n";
echo "log1p(0) = " . log1p(0) . "\n";      // 0
echo "log1p(e-1) = " . log1p(M_E - 1) . "\n"; // 1
?>

実践的な活用例

金融計算・複利計算

<?php
class FinancialCalculator {
    /**
     * 連続複利の計算
     * log1pを使用してより精密に計算
     */
    public static function continuousCompounding($principal, $rate, $time) {
        if ($rate <= -1) {
            throw new InvalidArgumentException("利率は-100%より大きい必要があります");
        }
        
        // A = P * exp(rt) = P * exp(log(1 + effective_rate))
        // 小さな利率の場合、log1pを使用することで精度が向上
        $log_growth = $rate * $time;
        
        if (abs($log_growth) < 0.1) {
            // 小さな値の場合はlog1pを活用
            $effective_rate = exp(log1p($log_growth - 1)) - 1;
        } else {
            $effective_rate = exp($log_growth) - 1;
        }
        
        return $principal * (1 + $effective_rate);
    }
    
    /**
     * 小さな利率での複利計算
     * 高精度計算が必要な場合
     */
    public static function preciseCompounding($principal, $annual_rate, $years, $compounds_per_year = 1) {
        $period_rate = $annual_rate / $compounds_per_year;
        $total_periods = $compounds_per_year * $years;
        
        if (abs($period_rate) < 1e-6) {
            // 非常に小さな利率の場合、log1pを使用
            $log_factor = $total_periods * log1p($period_rate);
            return $principal * exp($log_factor);
        } else {
            return $principal * pow(1 + $period_rate, $total_periods);
        }
    }
    
    /**
     * リターンの対数変換(対数リターン)
     */
    public static function calculateLogReturns($prices) {
        $log_returns = [];
        
        for ($i = 1; $i < count($prices); $i++) {
            $current_price = $prices[$i];
            $previous_price = $prices[$i - 1];
            
            if ($previous_price > 0 && $current_price > 0) {
                $simple_return = ($current_price - $previous_price) / $previous_price;
                
                // log(1 + return) = log1p(return) で高精度計算
                $log_return = log1p($simple_return);
                $log_returns[] = $log_return;
            }
        }
        
        return $log_returns;
    }
    
    /**
     * ボラティリティの計算
     */
    public static function calculateVolatility($log_returns, $annualization_factor = 252) {
        if (count($log_returns) < 2) {
            return 0;
        }
        
        $mean = array_sum($log_returns) / count($log_returns);
        $variance = 0;
        
        foreach ($log_returns as $return) {
            $variance += pow($return - $mean, 2);
        }
        
        $variance /= count($log_returns) - 1;
        $daily_volatility = sqrt($variance);
        
        // 年率換算
        return $daily_volatility * sqrt($annualization_factor);
    }
}

// 使用例
try {
    echo "金融計算の例:\n\n";
    
    // 小さな利率での複利計算
    $principal = 1000000; // 100万円
    $small_rate = 0.001;  // 0.1%(非常に低い利率)
    $years = 10;
    
    $standard_result = $principal * pow(1 + $small_rate, $years);
    $precise_result = FinancialCalculator::preciseCompounding($principal, $small_rate, $years);
    
    echo "元本: " . number_format($principal) . "円\n";
    echo "利率: " . ($small_rate * 100) . "%\n";
    echo "期間: {$years}年\n";
    echo "標準計算: " . number_format($standard_result, 2) . "円\n";
    echo "高精度計算: " . number_format($precise_result, 2) . "円\n";
    echo "差額: " . number_format(abs($precise_result - $standard_result), 6) . "円\n\n";
    
    // 株価データの対数リターン計算
    $stock_prices = [100, 101.2, 99.8, 102.1, 100.5, 103.2, 101.8];
    $log_returns = FinancialCalculator::calculateLogReturns($stock_prices);
    $volatility = FinancialCalculator::calculateVolatility($log_returns);
    
    echo "株価分析:\n";
    echo "株価データ: [" . implode(', ', $stock_prices) . "]\n";
    echo "対数リターン: [" . implode(', ', array_map(function($x) { return round($x, 6); }, $log_returns)) . "]\n";
    echo "年率ボラティリティ: " . round($volatility * 100, 2) . "%\n";
    
} catch (Exception $e) {
    echo "エラー: " . $e->getMessage() . "\n";
}
?>

統計分析・確率計算

<?php
class StatisticalAnalyzer {
    /**
     * 小さな確率での対数オッズの計算
     * log1pを使用して数値的安定性を確保
     */
    public static function calculateLogOdds($probability) {
        if ($probability <= 0 || $probability >= 1) {
            throw new InvalidArgumentException("確率は0と1の間である必要があります");
        }
        
        // log(p / (1-p)) = log(p) - log(1-p)
        // p が1に近い場合、1-p は非常に小さくなる
        
        if ($probability < 0.5) {
            // pが小さい場合
            return log($probability) - log1p(-$probability);
        } else {
            // pが大きい場合(1-pが小さい)
            return log($probability / (1 - $probability));
        }
    }
    
    /**
     * 相対リスク増加率の計算
     * log1pを使用してより精密に計算
     */
    public static function calculateRelativeRiskIncrease($baseline_risk, $new_risk) {
        if ($baseline_risk <= 0 || $new_risk <= 0) {
            throw new InvalidArgumentException("リスクは正の値である必要があります");
        }
        
        $risk_ratio = $new_risk / $baseline_risk;
        
        if (abs($risk_ratio - 1) < 0.1) {
            // リスク比が1に近い場合、log1pを使用
            return log1p($risk_ratio - 1);
        } else {
            return log($risk_ratio);
        }
    }
    
    /**
     * 情報量(自己情報)の計算
     */
    public static function calculateSelfInformation($probability, $base = M_E) {
        if ($probability <= 0 || $probability > 1) {
            throw new InvalidArgumentException("確率は0より大きく1以下である必要があります");
        }
        
        if ($base === M_E) {
            // 自然対数の場合
            return -log($probability);
        } elseif ($base === 2) {
            // 2を底とする場合(ビット)
            return -log($probability) / log(2);
        } else {
            return -log($probability) / log($base);
        }
    }
    
    /**
     * KLダイバージェンス(相対エントロピー)の計算
     */
    public static function calculateKLDivergence($p_distribution, $q_distribution) {
        if (count($p_distribution) !== count($q_distribution)) {
            throw new InvalidArgumentException("分布の要素数が一致しません");
        }
        
        $kl_divergence = 0;
        
        for ($i = 0; $i < count($p_distribution); $i++) {
            $p = $p_distribution[$i];
            $q = $q_distribution[$i];
            
            if ($p > 0 && $q > 0) {
                $ratio = $p / $q;
                
                if (abs($ratio - 1) < 1e-10) {
                    // p ≈ q の場合、log1p を使用
                    $log_ratio = log1p($ratio - 1);
                } else {
                    $log_ratio = log($ratio);
                }
                
                $kl_divergence += $p * $log_ratio;
            } elseif ($p > 0 && $q <= 0) {
                // qが0でpが正の場合、無限大
                return INF;
            }
            // p = 0の場合は寄与は0
        }
        
        return $kl_divergence;
    }
    
    /**
     * エントロピーの計算
     */
    public static function calculateEntropy($probabilities, $base = M_E) {
        $entropy = 0;
        
        foreach ($probabilities as $p) {
            if ($p > 0) {
                $log_p = ($base === M_E) ? log($p) : log($p) / log($base);
                $entropy -= $p * $log_p;
            }
        }
        
        return $entropy;
    }
}

// 使用例
try {
    echo "\n統計分析の例:\n\n";
    
    // 小さな確率での対数オッズ
    $small_probabilities = [0.001, 0.01, 0.1, 0.5, 0.9, 0.99, 0.999];
    
    echo "対数オッズの計算:\n";
    foreach ($small_probabilities as $p) {
        $log_odds = StatisticalAnalyzer::calculateLogOdds($p);
        echo "P = {$p}: 対数オッズ = " . round($log_odds, 6) . "\n";
    }
    
    // 相対リスク計算
    echo "\n相対リスク増加率:\n";
    $baseline_risk = 0.01; // 1%
    $risk_scenarios = [0.011, 0.012, 0.015, 0.02];
    
    foreach ($risk_scenarios as $new_risk) {
        $risk_increase = StatisticalAnalyzer::calculateRelativeRiskIncrease($baseline_risk, $new_risk);
        $percentage_increase = (exp($risk_increase) - 1) * 100;
        
        echo "ベースライン {$baseline_risk} -> {$new_risk}: ";
        echo "対数増加率 = " . round($risk_increase, 6);
        echo " (+" . round($percentage_increase, 2) . "%)\n";
    }
    
    // KLダイバージェンス
    echo "\nKLダイバージェンス:\n";
    $p_dist = [0.5, 0.3, 0.2];
    $q_dist = [0.4, 0.4, 0.2];
    
    $kl_div = StatisticalAnalyzer::calculateKLDivergence($p_dist, $q_dist);
    echo "P分布: [" . implode(', ', $p_dist) . "]\n";
    echo "Q分布: [" . implode(', ', $q_dist) . "]\n";
    echo "KL(P||Q) = " . round($kl_div, 6) . "\n";
    
} catch (Exception $e) {
    echo "エラー: " . $e->getMessage() . "\n";
}
?>

機械学習・最適化

<?php
class MLOptimizer {
    /**
     * ロジスティック回帰の対数尤度計算
     */
    public static function logisticLogLikelihood($predictions, $actual_labels) {
        $log_likelihood = 0;
        
        for ($i = 0; $i < count($predictions); $i++) {
            $p = $predictions[$i];
            $y = $actual_labels[$i];
            
            // 数値的安定性のためにlog1pを使用
            if ($y == 1) {
                $log_likelihood += log($p);
            } else {
                // log(1-p) を log1p(-p) として計算
                $log_likelihood += log1p(-$p);
            }
        }
        
        return $log_likelihood;
    }
    
    /**
     * 勾配降下法での学習率調整
     */
    public static function adaptiveLearningRate($base_rate, $iteration, $decay_rate = 0.01) {
        // 学習率 = base_rate / (1 + decay_rate * iteration)
        // 小さなdecay_rateの場合、log1pを使用してより精密に計算
        
        $decay_factor = $decay_rate * $iteration;
        
        if ($decay_factor < 0.1) {
            // 小さな減衰の場合
            $denominator = 1 + $decay_factor;
            $log_denominator = log1p($decay_factor);
            
            return $base_rate * exp(-$log_denominator);
        } else {
            return $base_rate / (1 + $decay_factor);
        }
    }
    
    /**
     * 交差エントロピー損失の計算
     */
    public static function crossEntropyLoss($predictions, $true_labels, $epsilon = 1e-15) {
        $total_loss = 0;
        $num_samples = count($predictions);
        
        for ($i = 0; $i < $num_samples; $i++) {
            $p = max(min($predictions[$i], 1 - $epsilon), $epsilon); // クリッピング
            $y = $true_labels[$i];
            
            if ($y == 1) {
                $total_loss -= log($p);
            } else {
                // log(1-p)をlog1p(-p)として計算(pが1に近い場合の精度向上)
                if ($p > 0.9) {
                    $total_loss -= log1p(-$p);
                } else {
                    $total_loss -= log(1 - $p);
                }
            }
        }
        
        return $total_loss / $num_samples;
    }
    
    /**
     * 正則化項の計算(L2正則化)
     */
    public static function l2Regularization($weights, $lambda) {
        $regularization = 0;
        
        foreach ($weights as $weight) {
            $regularization += $weight * $weight;
        }
        
        return 0.5 * $lambda * $regularization;
    }
    
    /**
     * Adam最適化アルゴリズムの更新
     */
    public static function adamUpdate($gradient, &$m, &$v, $beta1 = 0.9, $beta2 = 0.999, $epsilon = 1e-8, $t = 1) {
        // モーメント推定の更新
        $m = $beta1 * $m + (1 - $beta1) * $gradient;
        $v = $beta2 * $v + (1 - $beta2) * $gradient * $gradient;
        
        // バイアス補正
        $m_corrected = $m / (1 - pow($beta1, $t));
        $v_corrected = $v / (1 - pow($beta2, $t));
        
        // 更新量の計算
        $update = $m_corrected / (sqrt($v_corrected) + $epsilon);
        
        return $update;
    }
}

// 使用例
echo "\n機械学習での活用例:\n\n";

// ロジスティック回帰の対数尤度
$predictions = [0.1, 0.8, 0.6, 0.2, 0.9];
$actual_labels = [0, 1, 1, 0, 1];

$log_likelihood = MLOptimizer::logisticLogLikelihood($predictions, $actual_labels);
echo "ロジスティック回帰の対数尤度: " . round($log_likelihood, 6) . "\n";

// 適応的学習率
echo "\n学習率の変化:\n";
$base_learning_rate = 0.01;
for ($iteration = 0; $iteration <= 100; $iteration += 20) {
    $lr = MLOptimizer::adaptiveLearningRate($base_learning_rate, $iteration);
    echo "反復 {$iteration}: 学習率 = " . round($lr, 6) . "\n";
}

// 交差エントロピー損失
$predictions_ce = [0.9, 0.1, 0.8, 0.3, 0.95];
$true_labels = [1, 0, 1, 0, 1];

$ce_loss = MLOptimizer::crossEntropyLoss($predictions_ce, $true_labels);
echo "\n交差エントロピー損失: " . round($ce_loss, 6) . "\n";

// Adam最適化の例
echo "\nAdam最適化の例:\n";
$gradient = 0.1;
$m = 0;
$v = 0;

for ($t = 1; $t <= 5; $t++) {
    $update = MLOptimizer::adamUpdate($gradient, $m, $v, 0.9, 0.999, 1e-8, $t);
    echo "ステップ {$t}: 更新量 = " . round($update, 6) . ", m = " . round($m, 6) . ", v = " . round($v, 6) . "\n";
}
?>

数値解析・科学計算

<?php
class NumericalAnalysis {
    /**
     * 指数関数の高精度近似(テイラー展開)
     * 小さな値に対してlog1pと組み合わせて使用
     */
    public static function expApproximation($x, $terms = 10) {
        if (abs($x) > 1) {
            // 大きな値は標準のexp関数を使用
            return exp($x);
        }
        
        $result = 1;
        $term = 1;
        
        for ($n = 1; $n <= $terms; $n++) {
            $term *= $x / $n;
            $result += $term;
        }
        
        return $result;
    }
    
    /**
     * 数値微分(中心差分)
     * 小さなhに対する高精度計算
     */
    public static function numericalDerivative($func, $x, $h = 1e-8) {
        $f_plus = $func($x + $h);
        $f_minus = $func($x - $h);
        
        return ($f_plus - $f_minus) / (2 * $h);
    }
    
    /**
     * 複合台形公式による数値積分
     */
    public static function trapezoidalIntegration($func, $a, $b, $n = 1000) {
        $h = ($b - $a) / $n;
        $sum = ($func($a) + $func($b)) / 2;
        
        for ($i = 1; $i < $n; $i++) {
            $x = $a + $i * $h;
            $sum += $func($x);
        }
        
        return $h * $sum;
    }
    
    /**
     * log1pを使った高精度計算の例:
     * sinh(x) = (e^x - e^(-x))/2 の小さなxでの計算
     */
    public static function precisesinh($x) {
        if (abs($x) < 1e-8) {
            // 非常に小さな値の場合、テイラー展開を使用
            // sinh(x) = x + x³/6 + x⁵/120 + ...
            $x3 = $x * $x * $x;
            return $x + $x3 / 6;
        } elseif (abs($x) < 0.5) {
            // 小さな値の場合、log1pを活用した計算
            // sinh(x) = x * (1 + x²/6 + x⁴/120 + ...)
            $x2 = $x * $x;
            return $x * (1 + $x2 / 6 + $x2 * $x2 / 120);
        } else {
            return sinh($x);
        }
    }
    
    /**
     * 収束判定を含む反復計算
     */
    public static function iterativeMethod($initial_guess, $iteration_func, $tolerance = 1e-10, $max_iterations = 1000) {
        $x = $initial_guess;
        
        for ($i = 0; $i < $max_iterations; $i++) {
            $x_new = $iteration_func($x);
            
            // 相対誤差の計算にlog1pを使用
            $relative_error = abs($x_new - $x);
            
            if ($x != 0) {
                $relative_error /= abs($x);
            }
            
            if ($relative_error < $tolerance) {
                return [
                    'result' => $x_new,
                    'iterations' => $i + 1,
                    'converged' => true,
                    'final_error' => $relative_error
                ];
            }
            
            $x = $x_new;
        }
        
        return [
            'result' => $x,
            'iterations' => $max_iterations,
            'converged' => false,
            'final_error' => $relative_error ?? INF
        ];
    }
}

// 使用例
echo "\n数値解析の例:\n\n";

// 指数関数の近似
$test_values = [0.1, 0.01, 0.001];
echo "指数関数の近似:\n";
foreach ($test_values as $x) {
    $exact = exp($x);
    $approx = NumericalAnalysis::expApproximation($x);
    $error = abs($exact - $approx);
    
    echo "x = {$x}:\n";
    echo "  正確値: " . sprintf("%.10f", $exact) . "\n";
    echo "  近似値: " . sprintf("%.10f", $approx) . "\n";
    echo "  誤差:   " . sprintf("%.2e", $error) . "\n\n";
}

// 高精度sinh計算
echo "sinh関数の高精度計算:\n";
foreach ([0.001, 0.01, 0.1] as $x) {
    $standard = sinh($x);
    $precise = NumericalAnalysis::precisesinh($x);
    $error = abs($standard - $precise);
    
    echo "x = {$x}: 標準 = " . sprintf("%.10f", $standard);
    echo ", 高精度 = " . sprintf("%.10f", $precise);
    echo ", 差 = " . sprintf("%.2e", $error) . "\n";
}

// 数値微分の例
echo "\n数値微分の例 (f(x) = x²):\n";
$quadratic = function($x) { return $x * $x; };

$test_points = [1, 2, 3];
foreach ($test_points as $x) {
    $analytical = 2 * $x; // 解析的導関数
    $numerical = NumericalAnalysis::numericalDerivative($quadratic, $x);
    $error = abs($analytical - $numerical);
    
    echo "x = {$x}: 解析値 = {$analytical}, 数値解 = " . sprintf("%.10f", $numerical);
    echo ", 誤差 = " . sprintf("%.2e", $error) . "\n";
}

// 反復法の例(平方根の計算:ニュートン法)
echo "\n反復法(平方根計算):\n";
$target = 2;
$newton_iteration = function($x) use ($target) {
    return 0.5 * ($x + $target / $x);
};

$result = NumericalAnalysis::iterativeMethod(1.0, $newton_iteration);

echo "√{$target} の計算:\n";
echo "結果: " . sprintf("%.15f", $result['result']) . "\n";
echo "反復回数: " . $result['iterations'] . "\n";
echo "収束: " . ($result['converged'] ? 'yes' : 'no') . "\n";
echo "最終誤差: " . sprintf("%.2e", $result['final_error']) . "\n";
echo "真値との差: " . sprintf("%.2e", abs($result['result'] - sqrt($target))) . "\n";
?>

エラーハンドリングと特殊ケース

入力値の検証と例外処理

<?php
class SafeLog1p {
    /**
     * 安全なlog1p計算
     */
    public static function calculate($x) {
        // 入力値の検証
        if (!is_numeric($x)) {
            throw new InvalidArgumentException("数値が必要です: " . var_export($x, true));
        }
        
        $x = (float)$x;
        
        // 定義域の確認
        if ($x <= -1) {
            throw new DomainException("log1p(x)はx > -1で定義されます。入力値: {$x}");
        }
        
        // 特殊値の処理
        if (is_infinite($x)) {
            return $x > 0 ? INF : -INF;
        }
        
        if (is_nan($x)) {
            return NAN;
        }
        
        // 計算実行
        $result = log1p($x);
        
        // 結果の検証
        if (!is_finite($result) && is_finite($x) && $x > -1) {
            throw new RuntimeException("計算結果が予期しない値になりました");
        }
        
        return $result;
    }
    
    /**
     * 配列に対する安全な一括計算
     */
    public static function calculateArray($values) {
        $results = [];
        $errors = [];
        
        foreach ($values as $i => $value) {
            try {
                $results[$i] = self::calculate($value);
            } catch (Exception $e) {
                $results[$i] = null;
                $errors[$i] = $e->getMessage();
            }
        }
        
        return [
            'results' => $results,
            'errors' => $errors,
            'success_count' => count($results) - count($errors),
            'error_count' => count($errors)
        ];
    }
}

// テストケース
$test_values = [0, 0.1, -0.5, -0.9, -1, -1.1, 'abc', INF, -INF, NAN, 1e-10, 1e10];

echo "log1p関数の安全な実行テスト:\n";
foreach ($test_values as $value) {
    try {
        $result = SafeLog1p::calculate($value);
        echo "log1p(" . var_export($value, true) . ") = " . 
             (is_finite($result) ? sprintf("%.10f", $result) : $result) . "\n";
    } catch (Exception $e) {
        echo "エラー (" . var_export($value, true) . "): " . $e->getMessage() . "\n";
    }
}

// 配列での一括処理テスト
echo "\n配列での一括処理:\n";
$batch_values = [0.1, -0.5, -1.1, 0.001, 'invalid'];
$batch_results = SafeLog1p::calculateArray($batch_values);

echo "成功: " . $batch_results['success_count'] . "件\n";
echo "エラー: " . $batch_results['error_count'] . "件\n";

if (!empty($batch_results['errors'])) {
    echo "エラー詳細:\n";
    foreach ($batch_results['errors'] as $index => $error) {
        echo "  インデックス {$index}: {$error}\n";
    }
}
?>

精度テストとベンチマーク

<?php
class Log1pPerformanceAnalyzer {
    /**
     * 精度テスト:log1p vs log(1+x)
     */
    public static function precisionTest($test_ranges) {
        $results = [];
        
        foreach ($test_ranges as $range_name => $values) {
            $max_error = 0;
            $total_error = 0;
            $count = 0;
            
            foreach ($values as $x) {
                if ($x > -1) {
                    $log1p_result = log1p($x);
                    $log_result = log(1 + $x);
                    
                    $error = abs($log1p_result - $log_result);
                    $max_error = max($max_error, $error);
                    $total_error += $error;
                    $count++;
                }
            }
            
            $results[$range_name] = [
                'count' => $count,
                'max_error' => $max_error,
                'average_error' => $count > 0 ? $total_error / $count : 0,
                'relative_precision_gain' => $max_error > 0 ? log10(1 / $max_error) : INF
            ];
        }
        
        return $results;
    }
    
    /**
     * パフォーマンステスト
     */
    public static function performanceTest($iterations = 100000) {
        $test_value = 1e-10; // 小さな値でテスト
        
        // log1pのテスト
        $start_time = microtime(true);
        for ($i = 0; $i < $iterations; $i++) {
            log1p($test_value);
        }
        $log1p_time = microtime(true) - $start_time;
        
        // log(1+x)のテスト
        $start_time = microtime(true);
        for ($i = 0; $i < $iterations; $i++) {
            log(1 + $test_value);
        }
        $log_time = microtime(true) - $start_time;
        
        return [
            'iterations' => $iterations,
            'log1p_time' => $log1p_time,
            'log_time' => $log_time,
            'log1p_per_second' => $iterations / $log1p_time,
            'log_per_second' => $iterations / $log_time,
            'performance_ratio' => $log_time / $log1p_time
        ];
    }
    
    /**
     * メモリ使用量の分析
     */
    public static function memoryUsageTest($array_size = 10000) {
        $test_data = [];
        for ($i = 0; $i < $array_size; $i++) {
            $test_data[] = ($i + 1) * 1e-8;
        }
        
        // メモリ使用量測定開始
        $memory_start = memory_get_usage(true);
        
        // log1pでの計算
        $start_time = microtime(true);
        $log1p_results = array_map('log1p', $test_data);
        $log1p_time = microtime(true) - $start_time;
        $memory_after_log1p = memory_get_usage(true);
        
        // log(1+x)での計算
        $start_time = microtime(true);
        $log_results = array_map(function($x) { return log(1 + $x); }, $test_data);
        $log_time = microtime(true) - $start_time;
        $memory_after_log = memory_get_usage(true);
        
        return [
            'array_size' => $array_size,
            'log1p_time' => $log1p_time,
            'log_time' => $log_time,
            'memory_log1p' => $memory_after_log1p - $memory_start,
            'memory_log' => $memory_after_log - $memory_after_log1p,
            'throughput_log1p' => $array_size / $log1p_time,
            'throughput_log' => $array_size / $log_time
        ];
    }
}

// テストの実行
echo "\n精度・パフォーマンス分析:\n\n";

// 精度テスト
$test_ranges = [
    'very_small' => array_map(function($i) { return $i * 1e-15; }, range(1, 100)),
    'small' => array_map(function($i) { return $i * 1e-8; }, range(1, 100)),
    'medium' => array_map(function($i) { return $i * 0.001; }, range(1, 100)),
    'large' => array_map(function($i) { return $i * 0.1; }, range(1, 100))
];

echo "精度テスト結果:\n";
$precision_results = Log1pPerformanceAnalyzer::precisionTest($test_ranges);

foreach ($precision_results as $range => $result) {
    echo "{$range}:\n";
    echo "  サンプル数: " . $result['count'] . "\n";
    echo "  最大誤差: " . sprintf("%.2e", $result['max_error']) . "\n";
    echo "  平均誤差: " . sprintf("%.2e", $result['average_error']) . "\n";
    echo "  精度向上: " . (is_finite($result['relative_precision_gain']) ? 
                        round($result['relative_precision_gain'], 2) . "桁" : "無限大") . "\n\n";
}

// パフォーマンステスト
echo "パフォーマンステスト結果:\n";
$perf_results = Log1pPerformanceAnalyzer::performanceTest(50000);

echo "実行回数: " . number_format($perf_results['iterations']) . "\n";
echo "log1p実行時間: " . round($perf_results['log1p_time'] * 1000, 2) . " ms\n";
echo "log(1+x)実行時間: " . round($perf_results['log_time'] * 1000, 2) . " ms\n";
echo "log1p処理速度: " . number_format($perf_results['log1p_per_second'], 0) . " ops/sec\n";
echo "log(1+x)処理速度: " . number_format($perf_results['log_per_second'], 0) . " ops/sec\n";
echo "速度比: " . round($perf_results['performance_ratio'], 3) . "\n\n";

// メモリ使用量テスト
echo "メモリ使用量テスト結果:\n";
$memory_results = Log1pPerformanceAnalyzer::memoryUsageTest(5000);

echo "配列サイズ: " . number_format($memory_results['array_size']) . "\n";
echo "log1p処理時間: " . round($memory_results['log1p_time'] * 1000, 2) . " ms\n";
echo "log(1+x)処理時間: " . round($memory_results['log_time'] * 1000, 2) . " ms\n";
echo "log1pスループット: " . number_format($memory_results['throughput_log1p'], 0) . " items/sec\n";
echo "log(1+x)スループット: " . number_format($memory_results['throughput_log'], 0) . " items/sec\n";
?>

実用的なユーティリティクラス

<?php
class Log1pUtilities {
    /**
     * 複数の値に対するlog1pの合成
     * log1p(x₁) + log1p(x₂) + ... を効率的に計算
     */
    public static function sumLog1p($values) {
        $result = 0;
        $valid_count = 0;
        
        foreach ($values as $value) {
            if (is_numeric($value) && $value > -1) {
                $result += log1p($value);
                $valid_count++;
            }
        }
        
        return [
            'sum' => $result,
            'count' => $valid_count,
            'average' => $valid_count > 0 ? $result / $valid_count : 0
        ];
    }
    
    /**
     * 累積log1p計算
     */
    public static function cumulativeLog1p($values) {
        $cumulative = [];
        $running_sum = 0;
        
        foreach ($values as $value) {
            if (is_numeric($value) && $value > -1) {
                $running_sum += log1p($value);
            }
            $cumulative[] = $running_sum;
        }
        
        return $cumulative;
    }
    
    /**
     * 移動平均log1p
     */
    public static function movingAverageLog1p($values, $window_size) {
        if ($window_size <= 0 || $window_size > count($values)) {
            throw new InvalidArgumentException("無効なウィンドウサイズです");
        }
        
        $moving_averages = [];
        
        for ($i = 0; $i <= count($values) - $window_size; $i++) {
            $window = array_slice($values, $i, $window_size);
            $window_sum = 0;
            $valid_count = 0;
            
            foreach ($window as $value) {
                if (is_numeric($value) && $value > -1) {
                    $window_sum += log1p($value);
                    $valid_count++;
                }
            }
            
            $moving_averages[] = $valid_count > 0 ? $window_sum / $valid_count : 0;
        }
        
        return $moving_averages;
    }
    
    /**
     * 条件付きlog1p計算
     */
    public static function conditionalLog1p($values, $condition_func) {
        $results = [];
        
        foreach ($values as $value) {
            if ($condition_func($value) && is_numeric($value) && $value > -1) {
                $results[] = log1p($value);
            } else {
                $results[] = null;
            }
        }
        
        return $results;
    }
    
    /**
     * log1pの逆関数(expm1相当)
     */
    public static function inverseLog1p($log_values) {
        $results = [];
        
        foreach ($log_values as $log_value) {
            if (is_numeric($log_value)) {
                // log1p(x) = y なら x = expm1(y)
                $results[] = expm1($log_value);
            } else {
                $results[] = null;
            }
        }
        
        return $results;
    }
    
    /**
     * 統計情報付きのlog1p計算
     */
    public static function statisticalLog1p($values) {
        $log_values = [];
        $valid_values = [];
        
        foreach ($values as $value) {
            if (is_numeric($value) && $value > -1) {
                $log_val = log1p($value);
                $log_values[] = $log_val;
                $valid_values[] = $value;
            }
        }
        
        if (empty($log_values)) {
            return null;
        }
        
        sort($log_values);
        $count = count($log_values);
        
        return [
            'original_values' => $valid_values,
            'log_values' => $log_values,
            'count' => $count,
            'min' => $log_values[0],
            'max' => $log_values[$count - 1],
            'mean' => array_sum($log_values) / $count,
            'median' => $count % 2 === 0 ? 
                      ($log_values[$count/2 - 1] + $log_values[$count/2]) / 2 : 
                      $log_values[floor($count/2)],
            'std_dev' => self::calculateStdDev($log_values)
        ];
    }
    
    private static function calculateStdDev($values) {
        $mean = array_sum($values) / count($values);
        $variance = 0;
        
        foreach ($values as $value) {
            $variance += pow($value - $mean, 2);
        }
        
        return sqrt($variance / count($values));
    }
}

// 使用例
echo "\nユーティリティ関数の使用例:\n\n";

$sample_data = [0.01, 0.05, 0.1, 0.02, 0.08, 0.03, 0.12, 0.06];

// 合計log1p
$sum_result = Log1pUtilities::sumLog1p($sample_data);
echo "合計log1p:\n";
echo "  合計値: " . round($sum_result['sum'], 6) . "\n";
echo "  有効数: " . $sum_result['count'] . "\n";
echo "  平均値: " . round($sum_result['average'], 6) . "\n\n";

// 累積log1p
$cumulative = Log1pUtilities::cumulativeLog1p($sample_data);
echo "累積log1p: [" . implode(', ', array_map(function($x) { return round($x, 4); }, $cumulative)) . "]\n\n";

// 移動平均
$moving_avg = Log1pUtilities::movingAverageLog1p($sample_data, 3);
echo "移動平均log1p (ウィンドウ=3): [" . 
     implode(', ', array_map(function($x) { return round($x, 4); }, $moving_avg)) . "]\n\n";

// 条件付き計算(0.05以上の値のみ)
$conditional = Log1pUtilities::conditionalLog1p($sample_data, function($x) { return $x >= 0.05; });
$valid_conditional = array_filter($conditional, function($x) { return $x !== null; });
echo "条件付きlog1p (≥0.05): [" . 
     implode(', ', array_map(function($x) { return round($x, 4); }, $valid_conditional)) . "]\n\n";

// 統計情報
$stats = Log1pUtilities::statisticalLog1p($sample_data);
if ($stats) {
    echo "統計情報:\n";
    echo "  有効数: " . $stats['count'] . "\n";
    echo "  最小値: " . round($stats['min'], 6) . "\n";
    echo "  最大値: " . round($stats['max'], 6) . "\n";
    echo "  平均値: " . round($stats['mean'], 6) . "\n";
    echo "  中央値: " . round($stats['median'], 6) . "\n";
    echo "  標準偏差: " . round($stats['std_dev'], 6) . "\n";
}
?>

まとめ

log1p関数は、log(1 + x)を高精度で計算する特殊な数学関数で、以下のような特徴と利点があります:

主要な特徴

  • 高精度計算: 特に小さなx値に対してlog(1 + x)よりも数値的に安定
  • 数値的安定性: 浮動小数点演算における精度損失を最小化
  • 効率性: 専用の最適化されたアルゴリズムを使用

主要な活用分野

  1. 金融計算: 複利計算、リスク分析、対数リターンの計算
  2. 統計分析: 確率計算、対数オッズ、情報理論
  3. 機械学習: ロジスティック回帰、交差エントロピー損失
  4. 数値解析: 高精度数値計算、科学計算
  5. データサイエンス: 小さな変化量の分析

重要なポイント

  1. 定義域: x > -1でのみ定義
  2. 精度: 小さなx値でlog(1 + x)より高精度
  3. 用途: 金融・統計・機械学習分野での専門計算
  4. パフォーマンス: 専用最適化により高速

log1p関数は、高精度な数値計算が要求される専門分野において、計算の精度と安定性を大幅に向上させる重要なツールです。特に金融工学や機械学習において、小さな値の対数計算が必要な場面では必須の関数といえるでしょう。


この記事がPHPでの高精度数値計算の参考になれば幸いです。log1p関数の特性を理解して、より正確で信頼性の高い計算を実現してください。

タイトルとURLをコピーしました