科学技術計算.jlのソースを表示

== はじめに ==

本書は、科学技術計算に携わる研究者、技術者、学生のための実践的なハンドブックです。基礎的な数学から最新の機械学習手法まで、幅広い内容を網羅しています。各章には具体的な計算例と[[Julia]]によるプログラミング例を収録し、実務での即戦力となる知識の習得を目指しています。

== 第1章 基礎数学 ==

=== 1.1 代数学の基礎 ===

数学的な記述の基礎となる代数学について解説します。特に、科学技術計算で頻繁に現れる多項式、行列、ベクトルの取り扱いに重点を置きます。

==== 1.1.1 集合と写像 ====

集合論の基本概念から始めましょう。集合とは、ものの集まりを表す最も基本的な概念です。

; 集合 : 集合とは、明確な基準によって、ものの属する・属さないが判定できる、ものの集まりです。

例えば、以下のような集合が考えられます：
* 自然数の集合 <math>\mathbb{N} = \{1, 2, 3, \dots\}</math>
* 整数の集合 <math>\mathbb{Z} = \{\dots, -2, -1, 0, 1, 2, \dots\}</math>
* 実数の集合 <math>\mathbb{R}</math>

; 集合の演算
集合 <math>A</math>, <math>B</math> に対して、以下の基本的な演算が定義されます：

* '''和集合''' <math>A \cup B</math> ： <math>A</math> または <math>B</math> のいずれかに属する要素の集合
* '''積集合''' <math>A \cap B</math> ： <math>A</math> と <math>B</math> の両方に属する要素の集合
* '''差集合''' <math>A \setminus B</math> ： <math>A</math> に属し、<math>B</math> に属さない要素の集合

==== 1.1.2 群・環・体 ====

代数的構造の基本概念である群、環、体について説明します。

; 群 : 集合 <math>G</math> と二項演算 <math>\cdot</math> が与えられているとき、以下の4つの条件を満たすとき、<math>(G, \cdot)</math> を'''群'''といいます：

# '''結合法則'''：任意の <math>a, b, c \in G</math> に対して <math>(a \cdot b) \cdot c = a \cdot (b \cdot c)</math>
# '''単位元の存在'''：<math>e \in G</math> が存在して、任意の <math>a \in G</math> に対して<math>a \cdot e = e \cdot a = a</math>
# '''逆元の存在'''：任意の <math>a \in G</math> に対して、 <math>a \cdot b = b \cdot a = e</math> となる <math>b \in G</math> が存在
# '''閉じている'''：任意の <math>a, b \in G</math> に対して <math>a \cdot b \in G</math>

; 実践例：
:<syntaxhighlight lang=julia copy>
using LinearAlgebra

# 回転行列の群の例
rotation_matrix(theta) =[cos(theta) -sin(theta); sin(theta) cos(theta)]

# 2つの回転の合成
theta1 = π/4  # 45度
theta2 = π/3  # 60度
R1 = rotation_matrix(theta1)
R2 = rotation_matrix(theta2)

# 結合法則の確認
R3 = rotation_matrix(theta1 + theta2)
println(isapprox(R1 * R2, R3))  # true
</syntaxhighlight>

=== 1.2 微分積分学 ===

==== 1.2.1 微分の基礎 ====

微分は、関数の局所的な変化率を表す重要な概念です。科学技術計算において、物理現象の記述や最適化問題の解法に不可欠です。

; 導関数の定義
関数 <math> f(x) </math> の <math> x </math> における導関数 <math> f'(x) </math> は以下のように定義されます：

<math> f'(x) = \lim_{{h \to 0}} \frac{f(x+h) - f(x)}{h} </math>

; 偏微分
多変数関数 <math> f(x_1, x_2, \dots, x_n) </math> に対する <math> x_i </math> についての偏微分は、他の変数を定数とみなして行う微分です：

<math> \frac{\partial f}{\partial x_i} = \lim_{{h \to 0}} \frac{f(x_1, \dots, x_i + h, \dots, x_n) - f(x_1, \dots, x_i, \dots, x_n)}{h} </math>

; 実践例：
:<syntaxhighlight lang=julia copy>
using ForwardDiff

# 数値微分の例
f(x) = x^2 * sin(x)

# ForwardDiffを使用した数値微分
x = 2.0
dx = ForwardDiff.derivative(f, x)
println("f'($x) ≈ $dx")

# 中心差分法による独自実装
numerical_derivative(f, x, h=1e-6) = (f(x + h) - f(x - h)) / (2h)

dx_custom = numerical_derivative(f, x)
println("独自実装による f'($x) ≈ $dx_custom")
</syntaxhighlight>

==== 1.2.2 積分の基礎 ====

積分は、関数の累積的な効果や面積、体積を計算する際に用いられる基本的な演算です。

; 定積分の定義
関数 <math> f(x) </math> の区間 <math>[a, b]</math> における定積分は以下のように定義されます：

<math> \int_{a}^{b} f(x) \,dx = \lim_{{n \to \infty}} \sum_{i=1}^{n} f(x_i) \Delta x </math>

ここで、<math> \Delta x = \frac{b - a}{n} </math>であり、<math> x_i </math> は各小区間の代表点です。

; 数値積分の実装
:<syntaxhighlight lang=julia copy>
function trapezoidal_integration(f, a, b, n=1000)
    """台形則による数値積分"""
    x = range(a, stop=b, length=n)
    dx = (b - a) / (n - 1)
    return dx * (sum(f.(x)) - (f(a) + f(b))/2)
end

# 例：sin(x)の0からπまでの積分（結果は2になるはず）
g(x) = sin(x)

result = trapezoidal_integration(g, 0, π)
println("∫[0→π] sin(x)dx ≈ $result")

# より高度な方法：QuadGKの使用
using QuadGK
result_quadgk, error = quadgk(g, 0, π)
println("QuadGK使用: ∫[0→π] sin(x)dx ≈ $result_quadgk")
</syntaxhighlight>

==== 1.2.3 応用例：最適化問題 ====

微分積分学の重要な応用例として、最適化問題を考えましょう。特に、勾配降下法は機械学習でも頻繁に使用される手法です。

:<syntaxhighlight lang=julia copy>
function gradient_descent(f, df, x0, learning_rate=0.1, max_iter=100, tol=1e-6)
    """
    勾配降下法による最適化
    f: 目的関数
    df: fの導関数
    x0: 初期値
    """
    x = x0
    history = [x]

    for i in 1:max_iter
        gradient = df(x)
        x_new = x - learning_rate * gradient

        if abs(f(x_new) - f(x)) < tol
            break
        end

        x = x_new
        push!(history, x)
    end

    return x, history
end

# 例：x²の最小値を求める
f(x) = x^2
df(x) = 2x

x_min, history = gradient_descent(f, df, 2.0)
println("x²の最小値は x = $x_min で達成されます")
</syntaxhighlight>

=== 1.3 線形代数 ===

==== 1.3.1 ベクトルと行列の基礎 ====

線形代数は、科学技術計算の中で最も頻繁に使用される数学的道具の一つです。特に、大規模なデータ処理や機械学習において重要な役割を果たします。

==== 1.3.1 ベクトルと行列の基本演算 ====

ベクトルと行列は、多次元データを扱う上で基本となる数学的オブジェクトです。

; ベクトルの基本演算
* 内積：a·b = Σᵢ aᵢbᵢ
* ノルム：||a|| = √(Σᵢ aᵢ²)
* 外積（3次元の場合）：a×b = (a₂b₃-a₃b₂, a₃b₁-a₁b₃, a₁b₂-a₂b₁)

:<syntaxhighlight lang=julia copy>
using LinearAlgebra

# ベクトル演算の例
a = [1, 2, 3]
b = [4, 5, 6]

# 内積
dot_product = dot(a, b)
println("内積: $dot_product")

# ノルム（L2ノルム）
norm_a = norm(a)
println("ベクトルaのノルム: $norm_a")

# 外積
cross_product = cross(a, b)
println("外積: $cross_product")
</syntaxhighlight>

; 行列の基本演算
* 行列の和差：(A±B)ᵢⱼ = Aᵢⱼ±Bᵢⱼ
* 行列の積：(AB)ᵢⱼ = Σₖ AᵢₖBₖⱼ
* 転置：(Aᵀ)ᵢⱼ = Aⱼᵢ

:<syntaxhighlight lang=julia copy>
# 行列演算の例
A = [1 2; 3 4]
B = [5 6; 7 8]

# 行列の和
C = A + B
println("行列の和:\n", C)

# 行列の積
D = A * B
println("行列の積:\n", D)

# 転置
AT = transpose(A)
println("Aの転置:\n", AT)
</syntaxhighlight>

==== 1.3.2 線形方程式系の解法 ====

線形方程式系 <math>Ax = b</math> の解法は、科学技術計算における基本的な問題の一つです。

; ガウス消去法
:<syntaxhighlight lang=julia>
function gaussian_elimination(A, b)
    """
    ガウス消去法による線形方程式系の解法
    A: n×n係数行列
    b: n次元ベクトル
    """
    n = size(A, 1)
    Ab = [A b]  # 拡大係数行列

    # 前進消去
    for i in 1:n
        # ピボット選択
        pivot = argmax(abs.(Ab[i:n, i])) + i - 1
        Ab[i, :], Ab[pivot, :] = Ab[pivot, :], Ab[i, :]
        for j in i+1:n
            factor = Ab[j, i] / Ab[i, i]
            Ab[j, :] -= factor * Ab[i, :]
        end
    end

    # 後退代入
    x = zeros(n)
    for i in n:-1:1
        x[i] = (Ab[i, end] - sum(Ab[i, i+1:n] .* x[i+1:n])) / Ab[i, i]
    end

    return x
end

# 例題
A = [2 1 -1; -3 -1 2; -2 1 2]
b = [8, -11, -3]
x = gaussian_elimination(A, b)
println("解: ", x)
</syntaxhighlight>

==== 1.3.3 固有値と固有ベクトル ====

固有値問題は、振動解析や主成分分析など、多くの応用場面で現れます。

:<syntaxhighlight lang=julia>
function power_method(A; max_iter=1000, tol=1e-10)
    """
    べき乗法による最大固有値と対応する固有ベクトルの計算
    """
    n = size(A, 1)
    x = rand(n)
    x /= norm(x)
    
    for _ in 1:max_iter
        x_new = A * x
        x_new /= norm(x_new)
        if norm(x - x_new) < tol
            break
        end
        x = x_new
    end

    λ = (x' * A * x) / (x' * x)
    return λ, x
end

# 例題
A = [4 -2; -2 5]
eval_power, evec_power = power_method(A)
println("べき乗法による最大固有値: ", eval_power)
println("対応する固有ベクトル: ", evec_power)
</syntaxhighlight>

==== 1.3.4 応用例：画像の主成分分析 ====

線形代数の重要な応用例として、画像の次元削減を行う主成分分析を実装してみましょう。

:<syntaxhighlight lang=julia>
function pca_compression(image, n_components)
    """
    PCAによる画像圧縮
    image: グレースケール画像（2次元配列）
    n_components: 使用する主成分の数
    """
    mean_img = mean(image, dims=2)
    centered = image .- mean_img
    
    # 共分散行列の計算
    cov_matrix = centered * centered' / size(centered, 2)
    
    # 固有値分解
    eigenvalues, eigenvectors = eigen(cov_matrix)
    idx = sortperm(eigenvalues, rev=true)
    eigenvectors = eigenvectors[:, idx]
    
    # 次元削減
    components = eigenvectors[:, 1:n_components]
    projected = components' * centered
    reconstructed = components * projected .+ mean_img
    
    return reconstructed
end

# 例：ランダムなパターンの圧縮
using Random
Random.seed!(42)
original = rand(100, 1000)  # 100×1000のランダム画像
compressed = pca_compression(original, 10)

# 圧縮率と再構成誤差の計算
compression_ratio = length(original) / (length(compressed) + size(compressed, 2))
error = mean((original - compressed) .^ 2)
println("圧縮率: ", compression_ratio)
println("平均二乗誤差: ", error)
</syntaxhighlight>

=== 1.3.2 線形方程式系の解法 ===

線形方程式系 <math>Ax = b</math> の解法は、科学技術計算における基本的な問題の一つです。

==== ガウス消去法 ====

:<syntaxhighlight lang=julia copy>
function gaussian_elimination(A, b)
    """
    ガウス消去法による線形方程式系の解法
    A: n×n係数行列
    b: n次元ベクトル
    """
    n = size(A, 1)
    # 拡大係数行列を作成
    Ab = hcat(A, b)

    # 前進消去
    for i in 1:n
        # ピボット選択
        pivot = argmax(abs.(Ab[i:end, i])) + i - 1
        if pivot != i
            Ab[[i, pivot], :] = Ab[[pivot, i], :]
        end

        for j in (i + 1):n
            factor = Ab[j, i] / Ab[i, i]
            Ab[j, :] -= factor * Ab[i, :]
        end
    end

    # 後退代入
    x = zeros(n)
    for i in n:-1:1
        x[i] = (Ab[i, end] - dot(Ab[i, i+1:end-1], x[i+1:end])) / Ab[i, i]
    end

    return x
end

# 例題
A = [2.0 1.0 -1.0;
     -3.0 -1.0 2.0;
     -2.0 1.0 2.0]
b = [8.0, -11.0, -3.0]

x = gaussian_elimination(A, b)
println("解: ", x)

# 検証：Juliaのバックスラッシュ演算子を使用
x_julia = A \ b
println("Juliaのバックスラッシュ演算子による解: ", x_julia)
</syntaxhighlight>

=== 1.3.3 固有値と固有ベクトル ===

固有値問題は、振動解析や主成分分析など、多くの応用場面で現れます。

==== べき乗法 ====

:<syntaxhighlight lang=julia copy>
function power_method(A, max_iter=1000, tol=1e-10)
    """
    べき乗法による最大固有値と対応する固有ベクトルの計算
    """
    n = size(A, 1)
    x = rand(n)
    x = x / norm(x)

    for _ in 1:max_iter
        x_new = A * x
        x_new = x_new / norm(x_new)

        if all(isapprox.(x, x_new, rtol=tol))
            break
        end
        x = x_new
    end

    eigenvalue = dot(x, A * x) / dot(x, x)
    return eigenvalue, x
end

# 例：対称行列の固有値・固有ベクトル
A = [4.0 -2.0;
     -2.0 5.0]

# べき乗法による計算
eval_power, evec_power = power_method(A)
println("べき乗法による最大固有値: ", eval_power)
println("対応する固有ベクトル: ", evec_power)

# Juliaの固有値分解による計算（検証用）
evals, evecs = eigen(A)
println("\nJuliaのeigen関数による固有値: ", evals)
println("対応する固有ベクトル: ", evecs)
</syntaxhighlight>

=== 1.3.4 応用例：画像の主成分分析 ===

線形代数の重要な応用例として、画像の次元削減を行う主成分分析を実装してみましょう。

:<syntaxhighlight lang=julia copy>
function pca_compression(image, n_components)
    """
    PCAによる画像圧縮
    image: グレースケール画像（2次元配列）
    n_components: 使用する主成分の数
    """
    # 平均を引く
    mean = mean(image, dims=2)
    centered = image .- mean

    # 共分散行列の計算
    cov = centered * centered' / size(centered, 2)

    # 固有値分解
    eigenvalues, eigenvectors = eigen(cov)
    idx = sortperm(eigenvalues, rev=true)
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # 次元削減
    components = eigenvectors[:, 1:n_components]
    projected = components' * centered
    reconstructed = components * projected .+ mean

    return reconstructed
end

# 例：ランダムなパターンの圧縮
using Random
Random.seed!(42)
original = rand(100, 1000)  # 100×1000のランダム画像
compressed = pca_compression(original, 10)

# 圧縮率と再構成誤差の計算
compression_ratio = length(original) / (length(compressed) + size(compressed, 2))
error = mean((original .- compressed).^2)
println("圧縮率: ", compression_ratio)
println("平均二乗誤差: ", error)
</syntaxhighlight>

=== 1.4 複素解析 ===

複素解析は、科学技術計算において信号処理、制御理論、量子力学など、様々な分野で重要な役割を果たします。

==== 1.4.1 複素数の基礎 ====

:<syntaxhighlight lang=julia copy>
# 複素数の生成と基本演算
z1 = 2 + 3im
z2 = 1 - 2im

# 四則演算
println("和: ", z1 + z2)
println("差: ", z1 - z2)
println("積: ", z1 * z2)
println("商: ", z1 / z2)

# 極形式への変換
r = abs(z1)      # 絶対値（モジュラス）
theta = angle(z1)  # 偏角（アーギュメント）
println("絶対値: ", r)
println("偏角: ", theta, " rad")

# 極形式から直交形式への変換
z3 = r * (cos(theta) + im * sin(theta))
println("極形式から再構成: ", z3)
</syntaxhighlight>

==== 1.4.2 複素関数と解析性 ====

:<syntaxhighlight lang=julia copy>
function check_analytic(u, v, x, y, h=1e-6)
    """
    コーシー・リーマンの方程式を数値的に確認する
    """
    # 偏微分を数値的に計算
    du_dx = (u(x + h, y) - u(x - h, y)) / (2 * h)
    du_dy = (u(x, y + h) - u(x, y - h)) / (2 * h)
    dv_dx = (v(x + h, y) - v(x - h, y)) / (2 * h)
    dv_dy = (v(x, y + h) - v(x, y - h)) / (2 * h)

    println("∂u/∂x ≈ ", du_dx, ", ∂v/∂y ≈ ", dv_dy)
    println("∂u/∂y ≈ ", du_dy, ", -∂v/∂x ≈ ", -dv_dx)

    return isapprox(du_dx, dv_dy) && isapprox(du_dy, -dv_dx)
end

# 例：f(z) = z² = (x + yi)² の場合
u(x, y) = x^2 - y^2  # 実部
v(x, y) = 2*x*y      # 虚部

# z = 1 + i での解析性を確認
is_analytic = check_analytic(u, v, 1.0, 1.0)
println("関数は解析的か: ", is_analytic)
</syntaxhighlight>

==== 1.4.3 複素積分 ====

:<syntaxhighlight lang=julia copy>
function residue_theorem_example()
    """
    留数定理を用いた積分の例：
    ∮ (1/(z²+1)) dz を単位円周上で計算
    """
    # 極点: z = ±i
    # 単位円内の極: z = i
    # 留数: 1/(2i)

    using QuadGK

    integrand(theta) = begin
        # z = e^(iθ) とパラメータ化
        z = exp(im * theta)
        return 1 / (z^2 + 1) * im * z
    end

    # 数値積分（0 から 2π）
    result_numerical, _ = quadgk(integrand, 0, 2π)

    # 留数定理による結果
    result_residue = 2π * im * (1/(2im))

    println("数値積分結果: ", real(result_numerical))
    println("留数定理による結果: ", real(result_residue))
end

residue_theorem_example()
</syntaxhighlight>

==== 1.4.4 応用例：FFTを用いた信号処理 ====

:<syntaxhighlight lang=julia copy>
function signal_processing_example()
    """
    FFTを用いた信号のフィルタリング例
    """
    # サンプリングパラメータ
    t = LinRange(0, 1, 1000)
    dt = t[2] - t[1]
    freq = fftfreq(length(t), 1/dt)

    # 元の信号（2つの周波数成分）
    signal = sin.(2π*10*t) + 0.5*sin.(2π*50*t)

    # ノイズ追加
    noisy_signal = signal + randn(length(t)) * 0.2

    # FFTの実行
    fft_result = fft(noisy_signal)

    # ローパスフィルタの適用
    cutoff_freq = 30
    mask = abs.(freq) .< cutoff_freq
    filtered_fft = fft_result .* mask

    # 逆FFTで信号を復元
    filtered_signal = real(ifft(filtered_fft))

    return t, signal, noisy_signal, filtered_signal
end

# 信号処理の実行と結果の表示
t, orig, noisy, filtered = signal_processing_example()
println("信号処理完了")
</syntaxhighlight>

=== 1.5 確率・統計 ===

確率・統計は、データ分析、機械学習、実験データの解析など、幅広い応用場面で必要となる基礎知識です。

==== 1.5.1 確率の基礎 ====

:<syntaxhighlight lang=julia copy>
using Distributions

# サイコロを振る実験のシミュレーション
function dice_experiment(n_trials=10000)
    """
    公平なサイコロを振る実験のシミュレーション
    """
    results = rand(1:6, n_trials)
    probabilities = counts(results, 1:6) / n_trials

    println("実験的確率:")
    for face in 1:6
        println("目$face: ", probabilities[face])
    end

    # カイ二乗検定で一様性を確認
    chi2_stat = sum((counts(results, 1:6) .- n_trials/6).^2 / (n_trials/6))
    p_value = ccdf(Chisq(5), chi2_stat)

    println("\nカイ二乗検定:")
    println("統計量: ", chi2_stat)
    println("p値: ", p_value)
end

dice_experiment()
</syntaxhighlight>

==== 1.5.2 確率分布と統計量 ====

:<syntaxhighlight lang=julia copy>
function distribution_examples()
    """
    代表的な確率分布のサンプリングと可視化
    """
    using Random
    Random.seed!(42)
    n_samples = 1000

    # 正規分布
    normal_samples = rand(Normal(0, 1), n_samples)

    # 指数分布
    exponential_samples = rand(Exponential(1), n_samples)

    # 二項分布
    binomial_samples = rand(Binomial(10, 0.3), n_samples)

    # 基本統計量の計算
    function compute_statistics(samples, name)
        println("\n$nameの統計量:")
        println("平均: ", mean(samples))
        println("分散: ", var(samples))
        println("標準偏差: ", std(samples))
        println("中央値: ", median(samples))
    end

    compute_statistics(normal_samples, "正規分布")
    compute_statistics(exponential_samples, "指数分布")
    compute_statistics(binomial_samples, "二項分布")
end

distribution_examples()
</syntaxhighlight>

==== 1.5.3 統計的推定と検定 ====

:<syntaxhighlight lang=julia copy>
function statistical_inference_example()
    """
    統計的推定と検定の例
    """
    using Distributions
    using HypothesisTests

    # データ生成
    Random.seed!(42)
    population_mean = 100
    population_std = 15
    sample_size = 30

    sample = rand(Normal(population_mean, population_std), sample_size)

    # 平均の信頼区間
    confidence_level = 0.95
    t_dist = TDist(sample_size - 1)
    t_stat = quantile(t_dist, (1 + confidence_level) / 2)
    sample_mean = mean(sample)
    sample_std = std(sample)
    margin_error = t_stat * (sample_std / sqrt(sample_size))

    ci_lower = sample_mean - margin_error
    ci_upper = sample_mean + margin_error

    println("$confidence_level%信頼区間:")
    println("[", ci_lower, ", ", ci_upper, "]")

    # 仮説検定
    null_mean = 95  # 帰無仮説の平均値
    test_result = OneSampleTTest(sample, null_mean)

    println("\n一標本のt検定:")
    println("帰無仮説: μ = ", null_mean)
    println("t統計量: ", test_result.t)
    println("p値: ", pvalue(test_result))
end

statistical_inference_example()
</syntaxhighlight>

==== 1.5.4 回帰分析 ====

:<syntaxhighlight lang=julia copy>
using GLM
using DataFrames

function regression_analysis_example()
    """
    回帰分析の実装例
    """
    # データ生成
    Random.seed!(42)
    n_samples = 100

    # 説明変数
    X1 = rand(n_samples)
    X2 = rand(n_samples)

    # 目的変数（ノイズを含む）
    y = 2*X1 + 3*X2 + randn(n_samples) * 0.1

    # データフレームの作成
    data = DataFrame(X1=X1, X2=X2, y=y)

    # モデルの学習
    model = lm(@formula(y ~ X1 + X2), data)

    # 結果の表示
    println("回帰係数:")
    println(coef(model))

    println("\n決定係数 R²: ", r2(model))
    println("平均二乗誤差: ", deviance(model) / n_samples)
end

regression_analysis_example()
</syntaxhighlight>

==== 1.5.5 応用例：ブートストラップ法による信頼区間推定 ====

:<syntaxhighlight lang=julia copy>
using Random
using Statistics

function bootstrap_example()
    """
    ブートストラップ法による中央値の信頼区間推定
    """
    # データ生成
    Random.seed!(42)
    data = randn(100)
    data = exp.(data) # 対数正規分布に従う乱数

    # ブートストラップサンプリング
    n_bootstrap = 1000
    bootstrap_medians = zeros(n_bootstrap)

    for i in 1:n_bootstrap
        bootstrap_sample = rand(data, length(data)) # 重複ありランダムサンプリング
        bootstrap_medians[i] = median(bootstrap_sample)
    end

    # 信頼区間の計算（95%）
    ci_lower = quantile(bootstrap_medians, 0.025)
    ci_upper = quantile(bootstrap_medians, 0.975)

    println("中央値のブートストラップ信頼区間:")
    println("点推定値: ", median(data))
    println("95%信頼区間: [", ci_lower, ", ", ci_upper, "]")
end

bootstrap_example()
</syntaxhighlight>

== 第2章 数値計算の基礎 ==

=== 2.1 数値表現と誤差 ===

計算機における数値計算では、有限の精度で実数を表現する必要があり、そこから生じる誤差を適切に理解し管理することが重要です。

==== 2.1.1 浮動小数点数の表現 ====

IEEE 754規格による浮動小数点数の表現について理解しましょう。

:<syntaxhighlight lang=julia copy>
using Printf

function floating_point_example()
    """
    浮動小数点数の特性を確認する例
    """
    println("機械イプシロン: ", eps(Float64))
    println("最小の正規化数: ", typemin(Float64))
    println("最大の有限数: ", typemax(Float64))

    # 丸め誤差の例
    x = 0.1 + 0.2
    println("\n0.1 + 0.2 = ", x)
    println("0.3との差: ", x - 0.3)

    # 非正規化数の例
    x = Float64(1.0)
    for _ in 1:1074
        x /= 2
        @printf("x = %.2e\n", x)
    end
    @printf("最小の非正規化数: %.2e\n", x)
end

floating_point_example()
</syntaxhighlight>

==== 2.1.2 数値誤差の種類と対策 ====

:<syntaxhighlight lang=julia copy>
function error_analysis_example()
    """
    数値計算における主な誤差の例と対策
    """
    # 桁落ち（相対的に大きな数からの小さな数の減算）
    function catastrophic_cancellation()
        x = 1.0
        delta = 1e-16
        bad_result = (x + delta) - x
        println("桁落ちの例:")
        println("期待値: $delta")
        println("実際の結果: $bad_result")

        # 対策: 数式の変形
        better_result = delta
        println("改善後の結果: $better_result")
    end

    # 丸め誤差の蓄積
    function error_accumulation()
        sum1 = sum([0.1 for _ in 1:10])
        sum2 = 0.1 * 10
        println("\n丸め誤差の蓄積:")
        println("繰り返し加算: $sum1")
        println("乗算: $sum2")
        println("差: $(sum1 - sum2)")

        # 対策: Kahan summation algorithm
        function kahan_sum(numbers)
            s = 0.0
            c = 0.0  # 補償項
            for num in numbers
                y = num - c
                t = s + y
                c = (t - s) - y
                s = t
            end
            return s
        end

        sum3 = kahan_sum([0.1 for _ in 1:10])
        println("Kahan summation結果: $sum3")
    end

    catastrophic_cancellation()
    error_accumulation()
end

error_analysis_example()
</syntaxhighlight>

==== 2.1.3 条件数と数値安定性 ====

:<syntaxhighlight lang=julia copy>
function condition_number_example()
    """
    条件数と数値安定性の関係を示す例
    """
    # 行列の条件数と連立方程式の解の精度
    function matrix_condition()
        # 良条件の行列
        A1 = [2.0 1.0; 1.0 2.0]

        # 悪条件の行列
        A2 = [2.0 1.999; 1.999 2.0]

        b = [1.0, 1.0]

        println("条件数の比較:")
        println("良条件行列の条件数: $(cond(A1))")
        println("悪条件行列の条件数: $(cond(A2))")

        # 摂動を加えた場合の解の変化
        delta = 1e-10
        b_perturbed = b .+ delta

        x1 = A1 \ b
        x1_perturbed = A1 \ b_perturbed

        x2 = A2 \ b
        x2_perturbed = A2 \ b_perturbed

        println("\n摂動に対する解の変化:")
        println("良条件行列:")
        println("相対誤差: $(norm(x1_perturbed - x1) / norm(x1))")
        println("悪条件行列:")
        println("相対誤差: $(norm(x2_perturbed - x2) / norm(x2))")
    end

    matrix_condition()
end

condition_number_example()
</syntaxhighlight>

==== 2.1.4 実践的な対策と推奨事項 ====

:<syntaxhighlight lang=julia copy>
function numerical_best_practices()
    """
    数値計算における推奨プラクティス
    """
    # 1. スケーリング
    function scaling_example()
        # 悪い例
        x = [1e-8, 1e8]
        println("スケーリング前の数値範囲: $(minimum(x)) to $(maximum(x))")

        # 良い例
        x_scaled = x ./ maximum(abs.(x))
        println("スケーリング後の数値範囲: $(minimum(x_scaled)) to $(maximum(x_scaled))")
    end

    # 2. 数値的に安定なアルゴリズムの選択
    function stable_algorithm_example()
        # 悪い例：不安定な二次方程式の解法
        function unstable_quadratic(a, b, c)
            x1 = (-b + sqrt(b^2 - 4*a*c)) / (2*a)
            x2 = (-b - sqrt(b^2 - 4*a*c)) / (2*a)
            return x1, x2
        end

        # 良い例：数値的に安定な解法
        function stable_quadratic(a, b, c)
            if b >= 0
                x1 = (-b - sqrt(b^2 - 4*a*c)) / (2*a)
                x2 = c / (a*x1)
            else
                x1 = (-b + sqrt(b^2 - 4*a*c)) / (2*a)
                x2 = c / (a*x1)
            end
            return x1, x2
        end

        # 比較
        a, b, c = 1.0, 1e10, 1.0
        println("\n二次方程式の解法比較:")
        println("不安定な方法: $(unstable_quadratic(a, b, c))")
        println("安定な方法: $(stable_quadratic(a, b, c))")
    end

    scaling_example()
    stable_algorithm_example()
end

numerical_best_practices()
</syntaxhighlight>

=== 2.2 方程式の解法 ===

==== 2.2.1 非線形方程式の数値解法 ====

非線形方程式 f(x) = 0 の解を求めるための代表的な手法を実装します。

:<syntaxhighlight lang=julia copy>
function newton_method(f, df, x0; tol=1e-6, max_iter=100)
    """
    ニュートン法による非線形方程式の求解

    Parameters:
        f: 対象となる関数
        df: fの導関数
        x0: 初期値
        tol: 収束判定の閾値
        max_iter: 最大反復回数

    Returns:
        (解, 反復回数)
    """
    x = x0
    for i in 1:max_iter
        fx = f(x)
        if abs(fx) < tol
            return x, i
        end

        dfx = df(x)
        if abs(dfx) < 1e-10  # 導関数がゼロに近い場合
            error("導関数がゼロに近いため計算を中断します")
        end

        x_new = x - fx/dfx
        if abs(x_new - x) < tol
            return x_new, i + 1
        end
        x = x_new
    end

    error("最大反復回数に達しました")
end

function secant_method(f, x0, x1; tol=1e-6, max_iter=100)
    """
    割線法による非線形方程式の求解
    """
    for i in 1:max_iter
        f0, f1 = f(x0), f(x1)
        if abs(f1) < tol
            return x1, i
        end

        if abs(f1 - f0) < 1e-10
            error("割線の傾きがゼロに近いため計算を中断します")
        end

        x_new = x1 - f1 * (x1 - x0)/(f1 - f0)
        if abs(x_new - x1) < tol
            return x_new, i + 1
        end

        x0, x1 = x1, x_new
    end

    error("最大反復回数に達しました")
end

# 方程式を解く例
function equation_solving_example()
    """
    非線形方程式を解く例
    """
    # 例1: x³ - x - 1 = 0
    f1(x) = x^3 - x - 1
    df1(x) = 3*x^2 - 1

    try
        solution_newton, iterations = newton_method(f1, df1, 1.0)
        println("ニュートン法による解:")
        println("x ≈ $solution_newton")
        println("反復回数: $iterations")
        println("残差: $(abs(f1(solution_newton)))")
    catch e
        println("エラー: $e")
    end

    # 例2: cos(x) = x
    f2(x) = cos(x) - x

    try
        solution_secant, iterations = secant_method(f2, 0.0, 1.0)
        println("\n割線法による解:")
        println("x ≈ $solution_secant")
        println("反復回数: $iterations")
        println("残差: $(abs(f2(solution_secant)))")
    catch e
        println("エラー: $e")
    end
end

equation_solving_example()
</syntaxhighlight>

==== 2.2.2 連立方程式の解法 ====

:<syntaxhighlight lang=julia copy>
function linear_system_solvers()
    """
    連立一次方程式の各種解法の実装と比較
    """
    function gauss_seidel(A, b, x0; tol=1e-6, max_iter=1000)
        """
        ガウス・ザイデル法による連立一次方程式の求解
        """
        n = length(b)
        x = copy(x0)

        for it in 1:max_iter
            x_new = copy(x)
            for i in 1:n
                s1 = dot(A[i, 1:i-1], x_new[1:i-1])
                s2 = dot(A[i, i+1:end], x[i+1:end])
                x_new[i] = (b[i] - s1 - s2) / A[i, i]
            end

            if all(abs.(x - x_new) .< tol)
                return x_new, it + 1
            end
            x = x_new
        end

        error("最大反復回数に達しました")
    end

    # 例題: 3×3の連立方程式
    A = [4.0 -1.0 0.0; -1.0 4.0 -1.0; 0.0 -1.0 4.0]
    b = [1.0, 5.0, 0.0]
    x0 = zeros(length(b))

    # 1. ガウス・ザイデル法
    try
        solution_gs, iterations = gauss_seidel(A, b, x0)
        println("ガウス・ザイデル法による解:")
        println("x = $solution_gs")
        println("反復回数: $iterations")
        println("残差ノルム: $(norm(A * solution_gs - b))")
    catch e
        println("エラー: $e")
    end

    # 2. Juliaの直接解法との比較
    solution_julia = A \ b
    println("\nJulia線形ソルバーによる解:")
    println("x = $solution_julia")
    println("残差ノルム: $(norm(A * solution_julia - b))")

    # 3. 条件数と精度の関係
    cond_A = cond(A)
    println("\n行列の条件数: $cond_A")
    println("予想される相対誤差の上限: $(cond_A * eps(Float64))")
end

linear_system_solvers()
</syntaxhighlight>

==== 2.2.3 最適化問題の解法 ====

:<syntaxhighlight lang=julia copy>
function optimization_examples()
    """
    最適化問題の解法例
    """
    using Optim

    # 1. 最急降下法の実装
    function gradient_descent(f, grad_f, x0; learning_rate=0.1, tol=1e-6, max_iter=1000)
        """
        最急降下法による最適化
        """
        x = copy(x0)
        f_value = f(x)

        for i in 1:max_iter
            grad = grad_f(x)
            if norm(grad) < tol
                return x, f_value
            end

            x_new = x - learning_rate * grad
            f_value_new = f(x_new)

            if abs(f_value_new - f_value) < tol
                return x_new, f_value_new
            end

            x = x_new
            f_value = f_value_new
        end

        error("最大反復回数に達しました")
    end

    # 例題1: ロゼンブロック関数の最小化
    function rosenbrock(x)
        return (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
    end

    function rosenbrock_grad(x)
        return [-2*(1 - x[1]) - 400*x[1]*(x[2] - x[1]^2),
                200*(x[2] - x[1]^2)]
    end

    x0 = [-1.0, 1.0]

    # 最急降下法による解法
    try
        solution_gd, min_value = gradient_descent(rosenbrock, rosenbrock_grad, x0)
        println("最急降下法による最適化結果:")
        println("x = $solution_gd")
        println("f(x) = $min_value")
    catch e
        println("エラー: $e")
    end

    # Optim.jlの最適化ソルバーとの比較
    result = optimize(rosenbrock, x0, BFGS(), autodiff=:forward)
    println("\nOptim.jl BFGSによる最適化結果:")
    println("x = $(result.minimizer)")
    println("f(x) = $(result.minimum)")
    println("収束: $(Optim.converged(result))")
    println("反復回数: $(Optim.iterations(result))")
end

optimization_examples()
</syntaxhighlight>

=== 2.3 数値微分・積分 ===

==== 2.3.1 数値微分の手法 ====

:<syntaxhighlight lang=julia copy>
struct DerivativeResult
    value::Float64
    error_estimate::Float64
    step_size::Float64
end

function numerical_derivative(f::Function, x::Float64; method::String="central")
    eps = eps(Float64)
    h = sqrt(eps) * (1 + abs(x))  # 初期ステップサイズ
    
    if method == "forward"
        derivative = (f(x + h) - f(x)) / h
        error_est = abs(h)
        
    elseif method == "central"
        derivative = (f(x + h) - f(x - h)) / (2 * h)
        error_est = abs(h)^2
        
    elseif method == "richardson"
        # リチャードソン外挿
        h2 = h / 2
        d1 = (f(x + h) - f(x - h)) / (2 * h)
        d2 = (f(x + h2) - f(x - h2)) / (2 * h2)
        derivative = (4 * d2 - d1) / 3  # 4次精度の近似
        error_est = abs(d2 - d1) / 3
        
    else
        error("Unknown method")
    end
    
    return DerivativeResult(derivative, error_est, h)
end

function higher_derivatives(f::Function, x::Float64, order::Int=2)
    if order < 1
        error("Order must be positive")
    end
    
    h = cbrt(eps(Float64)) * (1 + abs(x))
    
    if order == 1
        return numerical_derivative(f, x, method="central").value
    elseif order == 2
        # 中心差分による2階導関数
        return (f(x + h) - 2*f(x) + f(x - h)) / (h * h)
    else
        # 再帰的に高次導関数を計算
        first_derivative(y) = numerical_derivative(f, y, method="central").value
        return higher_derivatives(first_derivative, x, order - 1)
    end
end

# 数値微分の例
function derivative_examples()
    # テスト関数
    f(x) = sin(x) * exp(-x/2)
    # 解析的な導関数
    df(x) = exp(-x/2) * (cos(x) - 0.5*sin(x))
    
    x = 1.0
    
    # 各手法での計算結果
    methods = ["forward", "central", "richardson"]
    println("数値微分の比較:")
    for method in methods
        result = numerical_derivative(f, x, method=method)
        println("\n$(method)差分法:")
        println("計算値: $(result.value)")
        println("真値との誤差: $(abs(result.value - df(x)))")
        println("推定誤差: $(result.error_estimate)")
        println("使用ステップサイズ: $(result.step_size)")
    end
    
    # 高次導関数の計算
    println("\n高次導関数:")
    for order in 1:3
        deriv = higher_derivatives(f, x, order)
        println("$(order)階導関数: $(deriv)")
    end
end

derivative_examples()
</syntaxhighlight>

==== 2.3.2 数値積分の手法 ====

:<syntaxhighlight lang=julia copy>
struct IntegrationResult
    value::Float64
    error_estimate::Float64
    n_evaluations::Int
end

function adaptive_simpson(f::Function, a::Float64, b::Float64; tol::Float64=1e-8, max_depth::Int=50)
    function simpson_rule(f, a, b)
        c = (a + b) / 2
        h = (b - a) / 6
        return h * (f(a) + 4*f(c) + f(b))
    end
    
    function _adaptive_simpson_recurse(a, b, fa, fb, fc, integral, tol, depth)
        c = (a + b) / 2
        h = (b - a) / 12
        d = (a + c) / 2
        e = (c + b) / 2
        fd = f(d)
        fe = f(e)
        
        left_integral = h * (fa + 4*fd + fc)
        right_integral = h * (fc + 4*fe + fb)
        whole_integral = left_integral + right_integral
        
        if depth >= max_depth
            return whole_integral, abs(whole_integral - integral)
        end
        
        if abs(whole_integral - integral) <= 15 * tol
            return whole_integral, abs(whole_integral - integral) / 15
        end
        
        left_result, left_error = _adaptive_simpson_recurse(a, c, fa, fc, fd, left_integral, tol/2, depth + 1)
        right_result, right_error = _adaptive_simpson_recurse(c, b, fc, fb, fe, right_integral, tol/2, depth + 1)
        
        return left_result + right_result, left_error + right_error
    end
    
    # 初期評価
    c = (a + b) / 2
    fa, fb, fc = f(a), f(b), f(c)
    initial_integral = simpson_rule(f, a, b)
    
    result, error = _adaptive_simpson_recurse(a, b, fa, fb, fc, initial_integral, tol, 0)
    
    return IntegrationResult(result, error, max_depth)
end

function gauss_quadrature(f::Function, a::Float64, b::Float64; n::Int=5)
    # ガウス・ルジャンドル求積の重みと節点（5点法）
    x5 = [0.0, ±0.5384693101056831, ±0.9061798459386640]
    w5 = [0.5688888888888889, 0.4786286704993665, 0.2369268850561891]
    
    # 区間変換
    mid = (b + a) / 2
    half_length = (b - a) / 2
    x = mid .+ half_length .* x5
    
    # 積分計算
    result = half_length * sum(w5 .* f.(x))
    
    # 誤差推定（簡易的）
    error_est = abs(result) * 1e-12
    
    return IntegrationResult(result, error_est, 5)
end

# 数値積分の例
function integration_examples()
    # テスト関数
    f(x) = exp(-x^2) * sin(2x)
    
    # 積分区間
    a, b = 0, 2
    
    # 適応的シンプソン法
    result_simpson = adaptive_simpson(f, a, b)
    println("適応的シンプソン法:")
    println("積分値: $(result_simpson.value)")
    println("推定誤差: $(result_simpson.error_estimate)")
    
    # ガウス求積法
    result_gauss = gauss_quadrature(f, a, b)
    println("\nガウス求積法:")
    println("積分値: $(result_gauss.value)")
    println("推定誤差: $(result_gauss.error_estimate)")
    
    # QuadGKとの比較
    using QuadGK
    result_quadgk, error_quadgk = quadgk(f, a, b)
    println("\nQuadGK:")
    println("積分値: $(result_quadgk)")
    println("推定誤差: $(error_quadgk)")
end

integration_examples()
</syntaxhighlight>

=== 2.4 線形方程式の解法 ===

==== 2.4.1 直接法 ====

:<syntaxhighlight lang=julia copy>
function lu_decomposition(A::Matrix{Float64})
    n = size(A, 1)
    L = zeros(n, n)
    U = zeros(n, n)
    
    for i in 1:n
        # L の対角成分
        L[i, i] = 1.0
        
        # U の i 行目
        for j in i:n
            sum_u = sum(L[i, k] * U[k, j] for k in 1:i-1)
            U[i, j] = A[i, j] - sum_u
        end
        
        # L の i 列目
        for j in i+1:n
            sum_l = sum(L[j, k] * U[k, i] for k in 1:i-1)
            L[j, i] = (A[j, i] - sum_l) / U[i, i]
        end
    end
    
    return L, U
end

function solve_lu(L::Matrix{Float64}, U::Matrix{Float64}, b::Vector{Float64})
    n = length(b)
    # 前進代入 (Ly = b)
    y = zeros(n)
    for i in 1:n
        y[i] = b[i] - dot(L[i, 1:i-1], y[1:i-1])
    end
    
    # 後退代入 (Ux = y)
    x = zeros(n)
    for i in n:-1:1
        x[i] = (y[i] - dot(U[i, i+1:n], x[i+1:n])) / U[i, i]
    end
    
    return x
end
</syntaxhighlight>

==== 2.4.2 反復法 ====

:<syntaxhighlight lang=julia copy>
function conjugate_gradient(A::Matrix{Float64}, b::Vector{Float64}; x0::Union{Vector{Float64}, Nothing}=nothing, tol::Float64=1e-10, max_iter::Int=1000)
    n = length(b)
    if x0 === nothing
        x = zeros(n)
    else
        x = copy(x0)
    end
    
    r = b - A * x
    p = copy(r)
    r_norm_sq = dot(r, r)
    
    for i in 1:max_iter
        Ap = A * p
        alpha = r_norm_sq / dot(p, Ap)
        x += alpha * p
        r -= alpha * Ap
        
        r_norm_sq_new = dot(r, r)
        if sqrt(r_norm_sq_new) < tol
            return x, i + 1
        end
        
        beta = r_norm_sq_new / r_norm_sq
        r_norm_sq = r_norm_sq_new
        p = r + beta * p
    end
    
    error("最大反復回数に達しました")
end
</syntaxhighlight>

==== 2.4.3 前処理付き反復法 ====

:<syntaxhighlight lang=julia copy>
function incomplete_cholesky(A::Matrix{Float64})
    n = size(A, 1)
    L = zeros(n, n)
    
    for i in 1:n
        for j in 1:i
            s = A[i, j]
            for k in 1:j-1
                s -= L[i, k] * L[j, k]
            end
            if i == j
                if s <= 0
                    error("行列が正定値でない可能性があります")
                end
                L[i, j] = sqrt(s)
            else
                L[i, j] = s / L[j, j]
            end
        end
    end
    
    return L
end

function preconditioned_cg(A::Matrix{Float64}, b::Vector{Float64}, M::Matrix{Float64}; x0::Union{Vector{Float64}, Nothing}=nothing, tol::Float64=1e-10, max_iter::Int=1000)
    n = length(b)
    if x0 === nothing
        x = zeros(n)
    else
        x = copy(x0)
    end
    
    r = b - A * x
    z = M \ r
    p = copy(z)
    rz = dot(r, z)
    
    for i in 1:max_iter
        Ap = A * p
        alpha = rz / dot(p, Ap)
        x += alpha * p
        r -= alpha * Ap
        
        if norm(r) < tol
            return x, i + 1
        end
        
        z = M \ r
        rz_new = dot(r, z)
        beta = rz_new / rz
        rz = rz_new
        p = z + beta * p
    end
    
    error("最大反復回数に達しました")
end
</syntaxhighlight>

=== 2.5 固有値計算 ===

==== 2.5.1 べき乗法 ====

:<syntaxhighlight lang=julia copy>
function power_method(A::Matrix{Float64}; tol::Float64=1e-10, max_iter::Int=1000)
    n = size(A, 1)
    x = rand(n)
    x = x / norm(x)
    
    lambda_old = 0
    
    for i in 1:max_iter
        # 行列とベクトルの積
        y = A * x
        
        # レイリー商による固有値の計算
        lambda_new = dot(x, y)
        
        # ベクトルの正規化
        x = y / norm(y)
        
        # 収束判定
        if abs(lambda_new - lambda_old) < tol
            return lambda_new, x
        end
        
        lambda_old = lambda_new
    end
    
    error("最大反復回数に達しました")
end
</syntaxhighlight>

== 第3章 信号処理 ==

=== 3.1 フーリエ解析 ===

フーリエ解析は信号処理の基礎となる重要な手法です。時間領域の信号を周波数領域で解析することで、信号の特性を理解し、効果的な処理が可能になります。

==== 3.1.1 離散フーリエ変換（DFT）の基礎 ====

:<syntaxhighlight lang=julia copy>
using Plots

function dft(x::AbstractVector{ComplexF64})
    """
    離散フーリエ変換の直接計算による実装

    Parameters:
        x: 入力信号

    Returns:
        X: 周波数領域の信号
    """
    N = length(x)
    X = zeros(ComplexF64, N)

    for k in 0:N-1
        for n in 0:N-1
            X[k+1] += x[n+1] * exp(-2im * pi * k * n / N)
        end
    end

    return X
end

function idft(X::AbstractVector{ComplexF64})
    """
    逆離散フーリエ変換の直接計算による実装

    Parameters:
        X: 周波数領域の信号

    Returns:
        x: 時間領域の信号
    """
    N = length(X)
    x = zeros(ComplexF64, N)

    for n in 0:N-1
        for k in 0:N-1
            x[n+1] += X[k+1] * exp(2im * pi * k * n / N)
        end
    end

    return x / N
end

function analyze_signal(t::AbstractVector{Float64}, x::AbstractVector{Float64})
    """
    信号のスペクトル解析を行う

    Parameters:
        t: 時間配列
        x: 信号配列

    Returns:
        freq: 周波数配列
        spectrum: パワースペクトル
    """
    N = length(t)
    dt = t[2] - t[1]  # サンプリング間隔

    # FFTの計算
    X = fft(x)

    # 周波数軸の生成
    freq = fftfreq(N, 1/dt)

    # パワースペクトルの計算
    spectrum = abs2.(X)

    return freq, spectrum
end
</syntaxhighlight>

==== 3.1.2 窓関数と周波数漏れ ====

:<syntaxhighlight lang=julia copy>
function apply_window(x::AbstractVector{Float64}, window_type::String="hanning")
    """
    信号に窓関数を適用する

    Parameters:
        x: 入力信号
        window_type: 窓関数の種類 ('hanning', 'hamming', 'blackman')

    Returns:
        x_windowed: 窓関数適用後の信号
    """
    N = length(x)

    if window_type == "hanning"
        window = 0.5 * (1 .- cos.(2 * pi * (0:N-1) / (N - 1)))
    elseif window_type == "hamming"
        window = 0.54 .- 0.46 * cos.(2 * pi * (0:N-1) / (N - 1))
    elseif window_type == "blackman"
        window = (0.42 .- 0.5 * cos.(2 * pi * (0:N-1) / (N - 1)) .+
                  0.08 * cos.(4 * pi * (0:N-1) / (N - 1)))
    else
        error("未対応の窓関数です")
    end

    return x .* window
end

function spectral_analysis_with_window(t::AbstractVector{Float64},
                                      x::AbstractVector{Float64},
                                      window_type::String="hanning")
    """
    窓関数を適用してスペクトル解析を行う

    Parameters:
        t: 時間配列
        x: 信号配列
        window_type: 窓関数の種類

    Returns:
        freq: 周波数配列
        spectrum_raw: 生のパワースペクトル
        spectrum_windowed: 窓関数適用後のパワースペクトル
    """
    # 生信号のスペクトル解析
    freq, spectrum_raw = analyze_signal(t, x)

    # 窓関数の適用
    x_windowed = apply_window(x, window_type)

    # 窓関数適用後の信号のスペクトル解析
    _, spectrum_windowed = analyze_signal(t, x_windowed)

    return freq, spectrum_raw, spectrum_windowed
end
</syntaxhighlight>

==== 3.1.3 短時間フーリエ変換（STFT） ====

:<syntaxhighlight lang=julia copy>
function stft(x::AbstractVector{Float64},
              window_size::Int64,
              hop_size::Int64,
              window_type::String="hanning")
    """
    短時間フーリエ変換の実装

    Parameters:
        x: 入力信号
        window_size: 窓サイズ
        hop_size: フレームシフト量
        window_type: 窓関数の種類

    Returns:
        times: 時間フレーム
        frequencies: 周波数ビン
        spectrogram: スペクトログラム
    """
    # 窓関数の生成
    window = apply_window(ones(window_size), window_type)

    # フレーム数の計算
    n_frames = 1 + (length(x) - window_size) ÷ hop_size

    # スペクトログラムの初期化
    spectrogram = zeros(ComplexF64, window_size ÷ 2 + 1, n_frames)

    # STFTの計算
    for i in 0:n_frames-1
        start = i * hop_size
        frame = x[start+1:start + window_size] .* window
        spectrum = rfft(frame)
        spectrogram[:, i+1] = spectrum
    end

    # 時間軸と周波数軸の生成
    times = (0:n_frames-1) * hop_size
    frequencies = rfftfreq(window_size, 1.0)

    return times, frequencies, spectrogram
end

function istft(spectrogram::AbstractMatrix{ComplexF64},
              window_size::Int64,
              hop_size::Int64,
              window_type::String="hanning")
    """
    逆短時間フーリエ変換の実装

    Parameters:
        spectrogram: スペクトログラム
        window_size: 窓サイズ
        hop_size: フレームシフト量
        window_type: 窓関数の種類

    Returns:
        x_reconstructed: 再構成された信号
    """
    # 窓関数の生成
    window = apply_window(ones(window_size), window_type)

    n_frames = size(spectrogram, 2)
    expected_signal_len = hop_size * (n_frames - 1) + window_size

    # 信号の再構成
    x_reconstructed = zeros(expected_signal_len)
    window_sum = zeros(expected_signal_len)

    for i in 0:n_frames-1
        start = i * hop_size
        frame = irfft(spectrogram[:, i+1], window_size)
        x_reconstructed[start+1:start + window_size] += frame .* window
        window_sum[start+1:start + window_size] += window .^ 2
    end

    # オーバーラップアド処理
    non_zero = window_sum .> 1e-10
    x_reconstructed[non_zero] ./= window_sum[non_zero]

    return x_reconstructed
end
</syntaxhighlight>

=== 3.2 ディジタルフィルタ ===

==== 3.2.1 FIRフィルタの設計 ====

:<syntaxhighlight lang=julia copy>
function design_fir_filter(cutoff::Float64,
                           filter_order::Int64,
                           filter_type::String="lowpass")
    """
    窓関数法によるFIRフィルタの設計

    Parameters:
        cutoff: カットオフ周波数 (0 < cutoff < 0.5)
        filter_order: フィルタの次数
        filter_type: フィルタの種類 ('lowpass' or 'highpass')

    Returns:
        h: フィルタ係数
    """
    if !(0 < cutoff < 0.5)
        error("カットオフ周波数は0と0.5の間である必要があります")
    end

    # インパルス応答の計算
    n = 0:filter_order
    omega_c = 2 * pi * cutoff

    if filter_type == "lowpass"
        h = sinc.(omega_c * ((n .- filter_order/2) / pi)
    elseif filter_type == "highpass"
        h = sinc.(n .- filter_order/2) .- sinc.(omega_c * (n .- filter_order/2) / pi)
    else
        error("未対応のフィルタ種類です")
    end

    # ハミング窓の適用
    window = 0.54 .- 0.46 * cos.(2 * pi * n / filter_order)
    h = h .* window

    # 正規化
    h = h / sum(h)

    return h
end

function apply_fir_filter(x::AbstractVector{Float64}, h::AbstractVector{Float64})
    """
    FIRフィルタの適用

    Parameters:
        x: 入力信号
        h: フィルタ係数

    Returns:
        y: フィルタリング後の信号
    """
    return conv(x, h, "same")
end
</syntaxhighlight>

==== 3.2.2 IIRフィルタの設計 ====

:<syntaxhighlight lang=julia copy>
using DSP

function design_iir_butterworth(cutoff::Float64,
                                 order::Int64,
                                 filter_type::String="lowpass")
    """
    バターワースIIRフィルタの設計

    Parameters:
        cutoff: カットオフ周波数 (0 < cutoff < 0.5)
        order: フィルタの次数
        filter_type: フィルタの種類

    Returns:
        b: 分子の係数
        a: 分母の係数
    """
    nyquist = 0.5
    cutoff_normalized = cutoff / nyquist

    if filter_type == "lowpass"
        b, a = butter(order, cutoff_normalized, lowpass=true)
    elseif filter_type == "highpass"
        b, a = butter(order, cutoff_normalized, highpass=true)
    else
        error("未対応のフィルタ種類です")
    end

    return b, a
end

function apply_iir_filter(x::AbstractVector{Float64},
                          b::AbstractVector{Float64},
                          a::AbstractVector{Float64},
                          initial_conditions::Union{AbstractVector{Float64}, Nothing}=nothing)
    """
    IIRフィルタの適用

    Parameters:
        x: 入力信号
        b: 分子の係数
        a: 分母の係数
        initial_conditions: 初期状態

    Returns:
        y: フィルタリング後の信号
    """
    if initial_conditions === nothing
        y = filt(b, a, x)
    else
        y = filt(b, a, x, zi=initial_conditions)
    end
    return y
end
</syntaxhighlight>

=== 3.3 画像処理の基礎 ===

==== 3.3.1 空間フィルタリング ====

:<syntaxhighlight lang=julia copy>
using Images

function spatial_filter(image::AbstractMatrix{Float64},
                        kernel::AbstractMatrix{Float64})
    """
    2次元畳み込みによる空間フィルタリング

    Parameters:
        image: 入力画像
        kernel: フィルタカーネル

    Returns:
        filtered_image: フィルタリング後の画像
    """
    # パディングサイズの計算
    pad_h = size(kernel, 1) ÷ 2
    pad_w = size(kernel, 2) ÷ 2

    # パディング処理
    padded_image = padarray(image, pad_h, pad_w, method=:reflect)

    # 出力画像の初期化
    filtered_image = zeros(Float64, size(image))

    # 畳み込み演算
    for i in 1:size(image, 1)
        for j in 1:size(image, 2)
            window = padded_image[i:i + size(kernel, 1) - 1, j:j + size(kernel, 2) - 1]
            filtered_image[i, j] = sum(window .* kernel)
        end
    end

    return filtered_image
end

function create_gaussian_kernel(size::Int64, sigma::Float64)
    """
    ガウシアンフィルタカーネルの生成

    Parameters:
        size: カーネルサイズ
        sigma: 標準偏差

    Returns:
        kernel: ガウシアンカーネル
    """
    x, y = meshgrid(-size ÷ 2:size ÷ 2, -size ÷ 2:size ÷ 2)

    kernel = exp.(-(x.^2 + y.^2) / (2 * sigma^2))
    return kernel / sum(kernel)
end
</syntaxhighlight>

==== 3.3.2 エッジ検出 ====

:<syntaxhighlight lang=julia copy>
function sobel_edge_detection(image::AbstractMatrix{Float64})
    """
    Sobelオペレータによるエッジ検出

    Parameters:
        image: 入力画像（グレースケール）

    Returns:
        magnitude: エッジの強度
        direction: エッジの方向
        edges: 二値化されたエッジ画像
    """
    # Sobelカーネルの定義
    kernel_x = [-1 0 1;
                -2 0 2;
                -1 0 1]

    kernel_y = [-1 -2 -1;
                0 0 0;
                1 2 1]

    # 勾配の計算
    gradient_x = spatial_filter(image, kernel_x)
    gradient_y = spatial_filter(image, kernel_y)

    # エッジの強度と方向の計算
    magnitude = sqrt.(gradient_x.^2 + gradient_y.^2)
    direction = atan.(gradient_y, gradient_x)

    # 閾値処理によるエッジの二値化
    threshold = 0.5 * maximum(magnitude)
    edges = magnitude .> threshold

    return magnitude, direction, edges
end

function canny_edge_detection(image::AbstractMatrix{Float64},
                              low_threshold::Float64=0.1,
                              high_threshold::Float64=0.3)
    """
    Cannyエッジ検出器の実装

    Parameters:
        image: 入力画像
        low_threshold: 低閾値（最大値に対する比率）
        high_threshold: 高閾値（最大値に対する比率）

    Returns:
        edges: 検出されたエッジ
    """
    # ステップ1: ガウシアンフィルタによるノイズ除去
    smoothed = spatial_filter(image, create_gaussian_kernel(5, 1.4))

    # ステップ2: 勾配の計算
    magnitude, direction, _ = sobel_edge_detection(smoothed)

    # ステップ3: 非最大値抑制
    suppressed = non_maximum_suppression(magnitude, direction)

    # ステップ4: ヒステリシス閾値処理
    high = high_threshold * maximum(suppressed)
    low = low_threshold * maximum(suppressed)

    return hysteresis_thresholding(suppressed, low, high)
end

function non_maximum_suppression(magnitude::AbstractMatrix{Float64},
                                  direction::AbstractMatrix{Float64})
    """
    非最大値抑制の実装
    """
    height, width = size(magnitude)
    suppressed = zeros(Float64, size(magnitude))

    # 勾配方向を0°、45°、90°、135°に量子化
    direction = rad2deg.(direction) .% 180

    for i in 2:height - 1
        for j in 2:width - 1
            if direction[i, j] < 22.5 || direction[i, j] >= 157.5
                neighbors = [magnitude[i, j - 1], magnitude[i, j + 1]]
            elseif 22.5 <= direction[i, j] < 67.5
                neighbors = [magnitude[i - 1, j - 1], magnitude[i + 1, j + 1]]
            elseif 67.5 <= direction[i, j] < 112.5
                neighbors = [magnitude[i - 1, j], magnitude[i + 1, j]]
            else
                neighbors = [magnitude[i - 1, j + 1], magnitude[i + 1, j - 1]]
            end

            if magnitude[i, j] >= maximum(neighbors)
                suppressed[i, j] = magnitude[i, j]
            end
        end
    end

    return suppressed
end
</syntaxhighlight>

==== 3.3.3 特徴点検出 ====

:<syntaxhighlight lang=julia copy>
function harris_corner_detector(image::AbstractMatrix{Float64},
                                 k::Float64=0.04, threshold::Float64=0.01)
    """
    Harrisコーナー検出器の実装

    Parameters:
        image: 入力画像
        k: Harrisレスポンス関数のパラメータ
        threshold: コーナー検出の閾値

    Returns:
        corners: コーナー点の位置（True/False）
    """
    # 勾配の計算
    dx = spatial_filter(image, [-1 0 1])
    dy = spatial_filter(image, [-1 0 1]')

    # 構造テンソルの要素
    Ixx = spatial_filter(dx.^2, create_gaussian_kernel(3, 1.0))
    Iyy = spatial_filter(dy.^2, create_gaussian_kernel(3, 1.0))
    Ixy = spatial_filter(dx .* dy, create_gaussian_kernel(3, 1.0))

    # Harrisレスポンスの計算
    det = Ixx .* Iyy .- Ixy.^2
    trace = Ixx .+ Iyy
    response = det .- k * trace.^2

    # 閾値処理と非最大値抑制
    corners = zeros(Bool, size(image))
    response_max = maximum(response)
    threshold_abs = threshold * response_max

    for i in 2:size(image, 1) - 1
        for j in 2:size(image, 2) - 1
            if response[i, j] > threshold_abs
                window = response[i - 1:i + 1, j - 1:j + 1]
                if response[i, j] == maximum(window)
                    corners[i, j] = true
                end
            end
        end
    end

    return corners
end
</syntaxhighlight>

=== 3.4 時系列解析 ===

==== 3.4.1 自己相関分析 ====

:<syntaxhighlight lang=julia copy>
function autocorrelation(x::AbstractVector{Float64},
                       max_lag::Union{Int64, Nothing}=nothing)
    """
    時系列データの自己相関関数を計算

    Parameters:
        x: 入力時系列データ
        max_lag: 最大ラグ（デフォルトはデータ長の半分）

    Returns:
        acf: 自己相関関数
    """
    if max_lag === nothing
        max_lag = length(x) ÷ 2
    end

    # 平均を引く
    x_centered = x .- mean(x)

    # 自己相関の計算
    acf = zeros(Float64, max_lag + 1)
    variance = var(x_centered)

    for lag in 0:max_lag
        acf[lag + 1] = sum(x_centered[lag + 1:end] .* x_centered[1:end - lag])
    end

    # 正規化
    return acf / (length(x) * variance)
end

function partial_autocorrelation(x::AbstractVector{Float64},
                                 max_lag::Union{Int64, Nothing}=nothing)
    """
    偏自己相関関数の計算（Durbin-Levinsonアルゴリズム）

    Parameters:
        x: 入力時系列データ
        max_lag: 最大ラグ

    Returns:
        pacf: 偏自己相関関数
    """
    if max_lag === nothing
        max_lag = length(x) ÷ 2
    end

    # 自己相関関数の計算
    acf = autocorrelation(x, max_lag)

    # 偏自己相関係数の初期化
    pacf = zeros(Float64, max_lag + 1)
    pacf[1] = 1.0

    # Durbin-Levinsonアルゴリズム
    for k in 1:max_lag
        phi = zeros(Float64, k)
        phi[k] = acf[k + 1]

        for j in 1:k - 1
            phi[j] = pacf[j]
        end

        for j in 1:k - 1
            phi[j] -= pacf[k] * pacf[k - j]
        end

        pacf[k] = phi[k]
    end

    return pacf
end
</syntaxhighlight>

==== 3.4.2 スペクトル推定 ====

:<syntaxhighlight lang=julia copy>
function welch_periodogram(x::AbstractVector{Float64},
                           window_size::Union{Int64, Nothing}=nothing,
                           overlap::Float64=0.5,
                           window_type::String="hanning")
    """
    Welch法によるパワースペクトル密度の推定

    Parameters:
        x: 入力時系列データ
        window_size: 窓サイズ
        overlap: オーバーラップ率（0～1）
        window_type: 窓関数の種類

    Returns:
        frequencies: 周波数配列
        psd: パワースペクトル密度
    """
    if window_size === nothing
        window_size = length(x) ÷ 8
    end

    # オーバーラップするサンプル数
    hop_size = Int64(floor(window_size * (1 - overlap)))

    # 窓関数の適用
    window = apply_window(ones(window_size), window_type)

    # セグメントの数の計算
    n_segments = (length(x) - window_size) ÷ hop_size + 1

    # スペクトルの初期化
    spectrum = zeros(Float64, window_size ÷ 2 + 1)

    # 各セグメントのスペクトルの計算と平均化
    for i in 0:n_segments - 1
        start = i * hop_size
        segment = x[start + 1:start + window_size] .* window
        segment_spectrum = abs2.(rfft(segment))
        spectrum += segment_spectrum
    end

    # 正規化
    spectrum /= n_segments

    # 周波数軸の生成
    frequencies = rfftfreq(window_size, 1.0)

    return frequencies, spectrum
end
</syntaxhighlight>

=== 3.5 ウェーブレット変換 ===

==== 3.5.1 連続ウェーブレット変換 ====

:<syntaxhighlight lang=julia copy>
function morlet_wavelet(t::AbstractVector{Float64}, omega0::Float64=6.0)
    """
    モルレーウェーブレットの生成

    Parameters:
        t: 時間配列
        omega0: 中心周波数

    Returns:
        psi: ウェーブレット関数
    """
    return pi^(-0.25) * exp.(1im * omega0 * t) .* exp.(-t.^2 / 2)
end

function cwt(x::AbstractVector{Float64},
             scales::AbstractVector{Float64},
             dt::Float64=1.0)
    """
    連続ウェーブレット変換の実装

    Parameters:
        x: 入力信号
        scales: スケールパラメータの配列
        dt: サンプリング間隔

    Returns:
        W: ウェーブレット係数（時間-スケール表現）
    """
    n = length(x)
    W = zeros(ComplexF64, length(scales), n)

    # 各スケールでのウェーブレット変換
    for i in 1:length(scales)
        scale = scales[i]
        # ウェーブレットの生成
        t = (-n ÷ 2 + 1:n ÷ 2) * dt
        wavelet = morlet_wavelet(t / scale)

        # 畳み込み（時間領域で計算）
        x_pad = padarray(x, n ÷ 2, method=:reflect)
        convolution = conv(x_pad, reverse(wavelet))
        W[i, :] = convolution[n ÷ 2 + 1:end - n ÷ 2]
    end

    return W
end
</syntaxhighlight>

----

<!--
== 第1章 基礎数学 ==
=== 代数学の基礎 ===
=== 微分積分学 ===
=== 線形代数 ===
=== 複素解析 ===
=== 確率・統計 ===

== 第2章 数値計算の基礎 ==
=== 数値表現と誤差 ===
=== 方程式の解法 ===
=== 数値微分・積分 ===
=== 線形方程式の解法 ===
=== 固有値計算 ===
=== 最適化手法 ===

== 第3章 信号処理 ==
=== フーリエ解析 ===
=== ディジタルフィルタ ===
=== 画像処理の基礎 ===
=== 時系列解析 ===
=== ウェーブレット変換 ===

== 第4章 シミュレーション技法 ==
=== 常微分方程式の解法 ===
=== 偏微分方程式の解法 ===
=== モンテカルロ法 ===
=== 分子動力学法 ===
=== 有限要素法 ===

== 第5章 データ解析と機械学習 ==
=== 回帰分析 ===
=== 主成分分析 ===
=== クラスタリング ===
=== ニューラルネットワーク ===
=== サポートベクターマシン ===

== 第6章 プログラミング実践 ==
=== アルゴリズムとデータ構造 ===
=== 科学技術計算言語（Python, MATLAB等） ===
=== 並列計算 ===
=== GPU計算 ===
=== 可視化技法 ===

== 第7章 応用例 ==
=== 物理学での応用 ===
=== 工学での応用 ===
=== 化学・生物学での応用 ===
=== 金融工学での応用 ===
=== 環境科学での応用 ===

附録
=== 数学公式集 ===
=== 物理定数表 ===
=== 単位換算表 ===
=== 参考文献 ===
=== 索引 ===
-->

----

== 下位階層のページ ==
{{特別:前方一致ページ一覧/{{PAGENAME}}/}}

{{DEFAULTSORT:かかくきしゆつけいさん JL}}
[[Category:科学技術計算|JL]]
[[Category:Julia]]