データサイエンス用の数理プログラム


```
import numpy as np
from sklearn.feature_extraction.text import CountVectorizer
from numpy.linalg import norm

# サンプルテキストデータ
documents = [
    "the cat sat on the mat",
    "the dog barked at the cat",
    "the cat chased the mouse",
    "the dog chased the ball"
]

# BoWベクトルの生成
vectorizer = CountVectorizer()
bow_matrix = vectorizer.fit_transform(documents).toarray()

# 単語のリストとBoW行列の表示
print("Vocabulary:", vectorizer.get_feature_names_out())
print("\nBoW Matrix:")
print(bow_matrix)

# 2つの文書間の内積とコサイン類似度を計算する関数
def cosine_similarity(vec1, vec2):
    dot_product = np.dot(vec1, vec2)
    norm1 = norm(vec1)
    norm2 = norm(vec2)
    cosine_sim = dot_product / (norm1 * norm2)
    return dot_product, cosine_sim

# 文書1と文書2のBoWベクトルを取り出す
vec1 = bow_matrix[0]  # "the cat sat on the mat"
vec2 = bow_matrix[1]  # "the dog barked at the cat"

# 内積とコサイン類似度の計算
dot_product, cosine_sim = cosine_similarity(vec1, vec2)

print(f"\nDot product between doc1 and doc2: {dot_product}")
print(f"Cosine similarity between doc1 and doc2: {cosine_sim}")

# 複数の文書間のコサイン類似度行列を計算
def cosine_similarity_matrix(matrix):
    similarity_matrix = np.dot(matrix, matrix.T) / (norm(matrix, axis=1).reshape(-1, 1) * norm(matrix, axis=1))
    return similarity_matrix

# コサイン類似度行列の計算と表示
cosine_sim_matrix = cosine_similarity_matrix(bow_matrix)
print("\nCosine Similarity Matrix:")
print(cosine_sim_matrix)

```

```
import numpy as np
import matplotlib.pyplot as plt

# データセットの作成（単純な線形データ）
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# モデルの初期パラメータ
theta = np.random.randn(2, 1)  # theta0（バイアス項）とtheta1（重み）

# 学習率と繰り返し数
learning_rate = 0.01
n_iterations = 1000
m = len(X)

# Xにバイアス項（1）を追加
X_b = np.c_[np.ones((m, 1)), X]

# 勾配降下法
for iteration in range(n_iterations):
    # 予測値を計算
    y_pred = X_b.dot(theta)
    
    # 損失関数の勾配（全微分）を計算
    gradients = 2/m * X_b.T.dot(y_pred - y)
    
    # パラメータの更新
    theta -= learning_rate * gradients

# 最終的なパラメータを出力
print(f"Final parameters: theta0 = {theta[0][0]}, theta1 = {theta[1][0]}")

# 結果のプロット
plt.scatter(X, y, color='blue', label='Data')
plt.plot(X, X_b.dot(theta), color='red', label='Linear Regression')
plt.xlabel('X')
plt.ylabel('y')
plt.title('Linear Regression with Gradient Descent')
plt.legend()
plt.show()


```


```
import autograd.numpy as np
from autograd import jacobian

# ベクトル値関数の定義（例：f(x) = [x1^2 + x2, x1 * x2]）
def vector_valued_function(x):
    return np.array([x[0]**2 + x[1], x[0] * x[1]])

# ヤコビ行列を計算する関数を生成
jacobian_matrix = jacobian(vector_valued_function)

# サンプル入力
x = np.array([3.0, 2.0])

# ヤコビ行列を計算
J = jacobian_matrix(x)

# 結果を出力
print("Jacobian matrix at x = [3.0, 2.0]:")
print(J)


```


```
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
from scipy.stats import norm

# 標準正規分布の確率密度関数（PDF）
def standard_normal_pdf(x):
    return norm.pdf(x)

# リーマン積分: 数値積分を近似
def riemann_integration(f, a, b, n=1000):
    x = np.linspace(a, b, n)
    dx = (b - a) / n
    return np.sum(f(x) * dx)

# scipyを使った数値積分（より精度の高い方法）
def scipy_integration(f, a, b):
    result, error = integrate.quad(f, a, b)
    return result

# 積分範囲
a, b = -5, 5

# リーマン積分の計算
riemann_result = riemann_integration(standard_normal_pdf, a, b)
print(f"Riemann Approximation: {riemann_result}")

# scipyの数値積分の計算（実質ルベーグ積分に相当）
scipy_result = scipy_integration(standard_normal_pdf, a, b)
print(f"SciPy Integration (Lebesgue-like): {scipy_result}")

# PDFのプロット
x = np.linspace(-5, 5, 1000)
y = standard_normal_pdf(x)

plt.plot(x, y, label='Standard Normal PDF', color='blue')
plt.fill_between(x, y, alpha=0.2, color='blue')
plt.title('Standard Normal Distribution PDF')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.show()

```

```
import autograd.numpy as np
from autograd import grad, hessian
import matplotlib.pyplot as plt

# 目的関数の定義
def function(x):
    return x**2 - 4*x + 4

# 一階および二階導関数の計算
gradient = grad(function)
hessian_function = hessian(function)

# ニュートン法による最小値探索
def newton_method(f, grad_f, hessian_f, x_init, tol=1e-6, max_iter=100):
    x = x_init
    for i in range(max_iter):
        grad_val = grad_f(x)
        hessian_val = hessian_f(x)
        
        # 更新式：x_{n+1} = x_n - f'(x) / f''(x)
        step = grad_val / hessian_val
        
        # xの更新
        x_new = x - step
        
        # 収束条件
        if np.abs(x_new - x) < tol:
            print(f"Converged after {i+1} iterations.")
            break
        
        x = x_new
    
    return x

# 初期値の設定
x_init = 0.0

# ニュートン法で最小値を探す
min_x = newton_method(function, gradient, hessian_function, x_init)

print(f"The minimum value is found at x = {min_x}")

# 関数のプロット
x_vals = np.linspace(-2, 6, 400)
y_vals = function(x_vals)

plt.plot(x_vals, y_vals, label='f(x) = x^2 - 4x + 4')
plt.scatter(min_x, function(min_x), color='red', label=f'Minimum at x = {min_x}', zorder=5)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Newton\'s Method for Minimization')
plt.legend()
plt.grid(True)
plt.show()


```


```
import sympy as sp

# 変数とラグランジュ乗数の定義
x, y, lamb = sp.symbols('x y lambda')

# 目的関数 f(x, y)
f = x**2 + y**2

# 制約条件 g(x, y) = 0
g = x + y - 1

# ラグランジュ関数の定義
L = f + lamb * g

# ラグランジュ関数の偏微分を計算
L_x = sp.diff(L, x)  # ∂L/∂x
L_y = sp.diff(L, y)  # ∂L/∂y
L_lamb = sp.diff(L, lamb)  # ∂L/∂λ

# 方程式のセットアップ
equations = [L_x, L_y, L_lamb]

# 方程式を解く
solution = sp.solve(equations, (x, y, lamb))

# 結果を表示
print(f"Optimal values: x = {solution[x]}, y = {solution[y]}, lambda = {solution[lamb]}")


```


```

import numpy as np
import matplotlib.pyplot as plt

# ネットワークの層数と各層のサイズ
n_layers = 50
layer_size = 100

# 初期化方法の設定（小さな乱数を用いる）
np.random.seed(42)
weights = [np.random.randn(layer_size, layer_size) * 0.01 for _ in range(n_layers)]

# 入力ベクトルの初期化
x = np.random.randn(layer_size)

# 勾配の初期化
gradients = []

# 順伝播と逆伝播のシミュレーション
for i, W in enumerate(weights):
    # 順伝播：シンプルに線形変換を行います
    x = W @ x  # 重みWを使って次の層へ

    # 活性化関数を適用（例: ReLU）
    x = np.maximum(0, x)

    # 勾配の計算（バックプロパゲーションのシミュレーション）
    grad = np.random.randn(layer_size)  # 出力側からのランダムな勾配を設定

    for j in range(i, -1, -1):
        # 重みの勾配を計算（ランダムな方向の影響を模倣）
        grad = weights[j].T @ grad

    # 勾配の大きさを記録
    gradients.append(np.linalg.norm(grad))

# 勾配の大きさをプロット
plt.plot(gradients, label="Gradient Norm")
plt.yscale('log')
plt.xlabel('Layer')
plt.ylabel('Gradient Magnitude (log scale)')
plt.title('Vanishing and Exploding Gradient in Deep Network')
plt.legend()
plt.grid(True)
plt.show()

```

```
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import convolve

# パラメータの設定
alpha = 0.01  # 熱伝導率
x = np.linspace(-10, 10, 400)  # 空間範囲
dx = x[1] - x[0]

# 初期温度分布 u0(x)
def initial_temperature(x):
    return np.exp(-x**2)

u0 = initial_temperature(x)

# ガウスカーネル関数 k(x, t)
def gaussian_kernel(x, t, alpha):
    return (1 / np.sqrt(4 * np.pi * alpha * t)) * np.exp(-x**2 / (4 * alpha * t))

# 時間 t での温度分布を計算
def heat_distribution(u0, x, t, alpha):
    kernel = gaussian_kernel(x, t, alpha)
    # 畳み込みを計算（境界を切り取るため 'same' を指定）
    u_t = convolve(u0, kernel, mode='same') * dx
    return u_t

# 時刻 t のリストを指定
times = [0.1, 0.5, 1.0, 2.0]

# 初期分布のプロット
plt.plot(x, u0, label='Initial: t=0', color='black')

# 各時間ステップごとに温度分布を計算してプロット
for t in times:
    u_t = heat_distribution(u0, x, t, alpha)
    plt.plot(x, u_t, label=f't={t}')

plt.xlabel('Position (x)')
plt.ylabel('Temperature (u)')
plt.title('Heat Propagation over Time')
plt.legend()
plt.grid(True)
plt.show()


```


```
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

# サンプルデータを生成（2次元データ）
np.random.seed(42)
n_samples = 100
mean = [0, 0]
cov = [[3, 1], [1, 2]]  # 共分散行列
data = np.random.multivariate_normal(mean, cov, n_samples)

# PCAを適用して潜在変数を取得
pca = PCA(n_components=1)  # 1次元に削減
latent_variables = pca.fit_transform(data)
reconstructed_data = pca.inverse_transform(latent_variables)

# 主成分ベクトル（固有ベクトル）を取得
principal_components = pca.components_

# 元のデータと主成分方向をプロット
plt.figure(figsize=(8, 6))
plt.scatter(data[:, 0], data[:, 1], alpha=0.6, label='Original Data')
plt.scatter(reconstructed_data[:, 0], reconstructed_data[:, 1], alpha=0.6, label='Reconstructed Data', color='orange')

# 主成分の方向をプロット
for i, (length, vector) in enumerate(zip(pca.explained_variance_, principal_components)):
    v = vector * 3 * np.sqrt(length)  # 固有値に基づいてベクトルの長さを調整
    plt.plot([0, v[0]], [0, v[1]], label=f'Principal Component {i+1}', linewidth=2)

plt.xlabel('X1')
plt.ylabel('X2')
plt.title('PCA: Original and Reconstructed Data')
plt.grid(True)
plt.legend()
plt.axis('equal')
plt.show()

# 主成分方向と潜在変数を表示
print("Principal Components:")
print(principal_components)
print("\nLatent Variables (First 5 samples):")
print(latent_variables[:5])


```


```
import numpy as np
import matplotlib.pyplot as plt

# シンプルなニューラルネットワークモデル（1層の線形モデル）
class SimpleNeuralNetwork:
    def __init__(self):
        # 重みとバイアスの初期化（ランダムな小さな値）
        self.weight = np.random.randn()
        self.bias = np.random.randn()

    # 順伝播（予測関数）
    def predict(self, x):
        return self.weight * x + self.bias

    # 損失関数（平均二乗誤差）
    def compute_loss(self, x, y):
        y_pred = self.predict(x)
        return np.mean((y - y_pred) ** 2)

    # 勾配を計算
    def compute_gradients(self, x, y):
        y_pred = self.predict(x)
        # 重みとバイアスに関する勾配の計算
        grad_weight = -2 * np.mean(x * (y - y_pred))
        grad_bias = -2 * np.mean(y - y_pred)
        return grad_weight, grad_bias

    # パラメータの更新
    def update_parameters(self, grad_weight, grad_bias, learning_rate):
        self.weight -= learning_rate * grad_weight
        self.bias -= learning_rate * grad_bias

# データの生成（線形データ）
np.random.seed(42)
x_train = np.linspace(-1, 1, 100)
y_train = 3 * x_train + 2 + np.random.randn(*x_train.shape) * 0.1  # y = 3x + 2 にノイズを加えたデータ

# モデルのインスタンス化
model = SimpleNeuralNetwork()

# ハイパーパラメータ
learning_rate = 0.1
n_epochs = 50

# 損失の履歴を保存
loss_history = []

# 勾配降下法によるトレーニング
for epoch in range(n_epochs):
    # 勾配の計算
    grad_weight, grad_bias = model.compute_gradients(x_train, y_train)
    
    # パラメータの更新
    model.update_parameters(grad_weight, grad_bias, learning_rate)
    
    # 損失の計算
    loss = model.compute_loss(x_train, y_train)
    loss_history.append(loss)
    print(f'Epoch {epoch+1}/{n_epochs}, Loss: {loss}')

# 学習結果の表示
plt.figure(figsize=(10, 5))

# 損失関数の履歴をプロット
plt.subplot(1, 2, 1)
plt.plot(range(n_epochs), loss_history, label='Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Loss over Epochs')
plt.grid(True)
plt.legend()

# 学習結果のプロット（予測線）
plt.subplot(1, 2, 2)
plt.scatter(x_train, y_train, label='Data', color='blue', alpha=0.6)
plt.plot(x_train, model.predict(x_train), color='red', label='Learned Model')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Data and Learned Linear Model')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# 学習されたパラメータを表示
print(f"Learned weight: {model.weight}")
print(f"Learned bias: {model.bias}")


```

```
import numpy as np
import matplotlib.pyplot as plt

# シンプルなニューラルネットワークモデル
class SimpleNeuralNetwork:
    def __init__(self):
        self.weight = np.random.randn()
        self.bias = np.random.randn()

    def predict(self, x):
        return self.weight * x + self.bias

    def compute_loss(self, x, y):
        y_pred = self.predict(x)
        return np.mean((y - y_pred) ** 2)

    def compute_gradients(self, x, y):
        y_pred = self.predict(x)
        grad_weight = -2 * np.mean(x * (y - y_pred))
        grad_bias = -2 * np.mean(y - y_pred)
        return grad_weight, grad_bias

    def update_parameters(self, grad_weight, grad_bias, learning_rate):
        self.weight -= learning_rate * grad_weight
        self.bias -= learning_rate * grad_bias

# データの生成（大規模なデータセット）
np.random.seed(42)
x_train = np.linspace(-1, 1, 10000)  # 10000サンプル
y_train = 3 * x_train + 2 + np.random.randn(*x_train.shape) * 0.1

# モデルのインスタンス化
model = SimpleNeuralNetwork()

# ハイパーパラメータ
learning_rate = 0.01
n_epochs = 50
batch_size = 32  # ミニバッチサイズ

# 損失の履歴を保存
loss_history = []

# ミニバッチ勾配降下法によるトレーニング
for epoch in range(n_epochs):
    # シャッフルしてバッチごとに処理
    indices = np.random.permutation(len(x_train))
    x_shuffled = x_train[indices]
    y_shuffled = y_train[indices]
    
    for start in range(0, len(x_train), batch_size):
        end = start + batch_size
        x_batch = x_shuffled[start:end]
        y_batch = y_shuffled[start:end]
        
        # 勾配の計算
        grad_weight, grad_bias = model.compute_gradients(x_batch, y_batch)
        
        # パラメータの更新
        model.update_parameters(grad_weight, grad_bias, learning_rate)
    
    # 全データでの損失を計算して記録
    loss = model.compute_loss(x_train, y_train)
    loss_history.append(loss)
    print(f'Epoch {epoch+1}/{n_epochs}, Loss: {loss}')

# 学習結果の表示
plt.plot(range(n_epochs), loss_history, label='Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Loss over Epochs (Mini-Batch Gradient Descent)')
plt.grid(True)
plt.legend()
plt.show()


```


```
import torch
import torch.nn as nn
import torch.optim as optim

# モデルの定義
class SimpleLinearModel(nn.Module):
    def __init__(self):
        super(SimpleLinearModel, self).__init__()
        self.linear = nn.Linear(1, 1)  # 入力1次元、出力1次元の線形層

    def forward(self, x):
        return self.linear(x)

# データの生成（大規模なデータセット）
x_train = torch.linspace(-1, 1, 10000).view(-1, 1)  # 10000サンプル
y_train = 3 * x_train + 2 + torch.randn(x_train.size()) * 0.1

# モデルのインスタンス化
model = SimpleLinearModel()

# 損失関数（平均二乗誤差）
criterion = nn.MSELoss()

# 最適化手法（勾配降下法）
optimizer = optim.SGD(model.parameters(), lr=0.01)

# トレーニング
n_epochs = 50
batch_size = 32
loss_history = []

for epoch in range(n_epochs):
    permutation = torch.randperm(x_train.size(0))
    
    for i in range(0, x_train.size(0), batch_size):
        indices = permutation[i:i + batch_size]
        x_batch, y_batch = x_train[indices], y_train[indices]
        
        # 順伝播
        y_pred = model(x_batch)
        loss = criterion(y_pred, y_batch)
        
        # 勾配をゼロに戻してから逆伝播とパラメータの更新
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
    
    # 全データでの損失を計算して記録
    with torch.no_grad():
        y_pred_full = model(x_train)
        full_loss = criterion(y_pred_full, y_train).item()
        loss_history.append(full_loss)
        print(f'Epoch {epoch+1}/{n_epochs}, Loss: {full_loss}')

# 学習結果の表示
plt.plot(range(n_epochs), loss_history, label='Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Loss over Epochs (PyTorch)')
plt.grid(True)
plt.legend()
plt.show()


```



```
import numpy as np

# 活性化関数とその導関数（ReLUを使用）
def relu(x):
    return np.maximum(0, x)

def relu_derivative(x):
    return np.where(x > 0, 1, 0)

# シグモイド関数とその導関数（出力層の活性化関数として使用）
def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def sigmoid_derivative(x):
    return sigmoid(x) * (1 - sigmoid(x))

# 損失関数（平均二乗誤差）
def compute_loss(y_true, y_pred):
    return np.mean((y_true - y_pred) ** 2)

# ニューラルネットワークの初期化
np.random.seed(42)
input_size = 2
hidden_size = 3
output_size = 1

# 重みとバイアスの初期化
W1 = np.random.randn(input_size, hidden_size)
b1 = np.random.randn(hidden_size)
W2 = np.random.randn(hidden_size, output_size)
b2 = np.random.randn(output_size)

# 順伝播
def forward(x):
    z1 = x @ W1 + b1  # 入力から隠れ層への線形変換
    a1 = relu(z1)     # 活性化関数（ReLU）
    z2 = a1 @ W2 + b2  # 隠れ層から出力層への線形変換
    a2 = sigmoid(z2)   # 出力層の活性化関数（シグモイド）
    return z1, a1, z2, a2

# 逆伝播
def backward(x, y_true, z1, a1, z2, a2, learning_rate=0.01):
    global W1, b1, W2, b2

    # 出力層の勾配
    dL_da2 = -(y_true - a2)  # 損失関数の導関数（MSEの導関数）
    da2_dz2 = sigmoid_derivative(z2)  # シグモイドの導関数
    dL_dz2 = dL_da2 * da2_dz2

    # 隠れ層から出力層への勾配
    dL_dW2 = a1.T @ dL_dz2  # W2に対する勾配
    dL_db2 = np.sum(dL_dz2, axis=0)  # b2に対する勾配

    # 隠れ層の勾配
    dL_da1 = dL_dz2 @ W2.T
    da1_dz1 = relu_derivative(z1)
    dL_dz1 = dL_da1 * da1_dz1

    # 入力層から隠れ層への勾配
    dL_dW1 = x.T @ dL_dz1  # W1に対する勾配
    dL_db1 = np.sum(dL_dz1, axis=0)  # b1に対する勾配

    # パラメータの更新
    W1 -= learning_rate * dL_dW1
    b1 -= learning_rate * dL_db1
    W2 -= learning_rate * dL_dW2
    b2 -= learning_rate * dL_db2

# サンプルデータ（XOR問題）
x_train = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
y_train = np.array([[0], [1], [1], [0]])  # XORのターゲット

# トレーニング
n_epochs = 10000
learning_rate = 0.1

for epoch in range(n_epochs):
    z1, a1, z2, a2 = forward(x_train)
    loss = compute_loss(y_train, a2)
    backward(x_train, y_train, z1, a1, z2, a2, learning_rate)

    if epoch % 1000 == 0:
        print(f'Epoch {epoch}, Loss: {loss}')

# 学習結果の表示
print("Trained weights and biases:")
print(f"W1: {W1}")
print(f"b1: {b1}")
print(f"W2: {W2}")
print(f"b2: {b2}")

# 最終的な予測
_, _, _, y_pred = forward(x_train)
print("Predictions after training:")
print(y_pred)



```



```

import numpy as np
import matplotlib.pyplot as plt

# シンプルな線形回帰モデル
class LinearRegression:
    def __init__(self):
        self.weight = np.random.randn()
        self.bias = np.random.randn()
        # Adam用の初期パラメータ
        self.m_w = 0  # 勾配の一次モーメント（重み用）
        self.v_w = 0  # 勾配の二次モーメント（重み用）
        self.m_b = 0  # 勾配の一次モーメント（バイアス用）
        self.v_b = 0  # 勾配の二次モーメント（バイアス用）

    def predict(self, x):
        return self.weight * x + self.bias

    def compute_loss(self, x, y):
        y_pred = self.predict(x)
        return np.mean((y - y_pred) ** 2)

    def compute_gradients(self, x, y):
        y_pred = self.predict(x)
        grad_weight = -2 * np.mean(x * (y - y_pred))
        grad_bias = -2 * np.mean(y - y_pred)
        return grad_weight, grad_bias

    def update_parameters(self, grad_weight, grad_bias, learning_rate, beta1, beta2, epsilon, t):
        # モーメントの更新
        self.m_w = beta1 * self.m_w + (1 - beta1) * grad_weight
        self.v_w = beta2 * self.v_w + (1 - beta2) * (grad_weight ** 2)
        self.m_b = beta1 * self.m_b + (1 - beta1) * grad_bias
        self.v_b = beta2 * self.v_b + (1 - beta2) * (grad_bias ** 2)

        # バイアス補正
        m_w_hat = self.m_w / (1 - beta1 ** t)
        v_w_hat = self.v_w / (1 - beta2 ** t)
        m_b_hat = self.m_b / (1 - beta1 ** t)
        v_b_hat = self.v_b / (1 - beta2 ** t)

        # パラメータの更新
        self.weight -= learning_rate * m_w_hat / (np.sqrt(v_w_hat) + epsilon)
        self.bias -= learning_rate * m_b_hat / (np.sqrt(v_b_hat) + epsilon)

# データの生成
np.random.seed(42)
x_train = np.linspace(-1, 1, 100)
y_train = 3 * x_train + 2 + np.random.randn(*x_train.shape) * 0.1

# モデルのインスタンス化
model = LinearRegression()

# ハイパーパラメータ
learning_rate = 0.01
beta1 = 0.9
beta2 = 0.999
epsilon = 1e-8
n_epochs = 1000

# 損失の履歴を保存
loss_history = []

# Adamによるトレーニング
for epoch in range(1, n_epochs + 1):
    grad_weight, grad_bias = model.compute_gradients(x_train, y_train)
    model.update_parameters(grad_weight, grad_bias, learning_rate, beta1, beta2, epsilon, epoch)
    
    # 損失の計算
    loss = model.compute_loss(x_train, y_train)
    loss_history.append(loss)
    
    if epoch % 100 == 0:
        print(f'Epoch {epoch}, Loss: {loss}')

# 学習結果の表示
plt.figure(figsize=(10, 5))

# 損失関数の履歴をプロット
plt.subplot(1, 2, 1)
plt.plot(range(n_epochs), loss_history, label='Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Loss over Epochs (Adam)')
plt.grid(True)
plt.legend()

# 学習結果のプロット（予測線）
plt.subplot(1, 2, 2)
plt.scatter(x_train, y_train, label='Data', color='blue', alpha=0.6)
plt.plot(x_train, model.predict(x_train), color='red', label='Learned Model')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Data and Learned Linear Model')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# 学習されたパラメータを表示
print(f"Learned weight: {model.weight}")
print(f"Learned bias: {model.bias}")


```



```
import numpy as np
import matplotlib.pyplot as plt

# シンプルな線形回帰モデル
class LinearRegression:
    def __init__(self):
        self.weight = np.random.randn()
        self.bias = np.random.randn()

        # Momentum用の初期パラメータ
        self.velocity_w = 0
        self.velocity_b = 0

        # RMSProp用の初期パラメータ
        self.s_w = 0
        self.s_b = 0

        # Adam用の初期パラメータ
        self.m_w = 0
        self.v_w = 0
        self.m_b = 0
        self.v_b = 0

    def predict(self, x):
        return self.weight * x + self.bias

    def compute_loss(self, x, y):
        y_pred = self.predict(x)
        return np.mean((y - y_pred) ** 2)

    def compute_gradients(self, x, y):
        y_pred = self.predict(x)
        grad_weight = -2 * np.mean(x * (y - y_pred))
        grad_bias = -2 * np.mean(y - y_pred)
        return grad_weight, grad_bias

    def update_parameters_momentum(self, grad_weight, grad_bias, learning_rate, beta=0.9):
        # Momentumの更新
        self.velocity_w = beta * self.velocity_w + (1 - beta) * grad_weight
        self.velocity_b = beta * self.velocity_b + (1 - beta) * grad_bias
        self.weight -= learning_rate * self.velocity_w
        self.bias -= learning_rate * self.velocity_b

    def update_parameters_rmsprop(self, grad_weight, grad_bias, learning_rate, beta=0.9, epsilon=1e-8):
        # RMSPropの更新
        self.s_w = beta * self.s_w + (1 - beta) * (grad_weight ** 2)
        self.s_b = beta * self.s_b + (1 - beta) * (grad_bias ** 2)
        self.weight -= learning_rate * grad_weight / (np.sqrt(self.s_w) + epsilon)
        self.bias -= learning_rate * grad_bias / (np.sqrt(self.s_b) + epsilon)

    def update_parameters_adam(self, grad_weight, grad_bias, learning_rate, beta1=0.9, beta2=0.999, epsilon=1e-8, t=1):
        # Adamの更新
        self.m_w = beta1 * self.m_w + (1 - beta1) * grad_weight
        self.v_w = beta2 * self.v_w + (1 - beta2) * (grad_weight ** 2)
        self.m_b = beta1 * self.m_b + (1 - beta1) * grad_bias
        self.v_b = beta2 * self.v_b + (1 - beta2) * (grad_bias ** 2)

        # バイアス補正
        m_w_hat = self.m_w / (1 - beta1 ** t)
        v_w_hat = self.v_w / (1 - beta2 ** t)
        m_b_hat = self.m_b / (1 - beta1 ** t)
        v_b_hat = self.v_b / (1 - beta2 ** t)

        self.weight -= learning_rate * m_w_hat / (np.sqrt(v_w_hat) + epsilon)
        self.bias -= learning_rate * m_b_hat / (np.sqrt(v_b_hat) + epsilon)

# データの生成
np.random.seed(42)
x_train = np.linspace(-1, 1, 100)
y_train = 3 * x_train + 2 + np.random.randn(*x_train.shape) * 0.1

# モデルのインスタンス化
model_momentum = LinearRegression()
model_rmsprop = LinearRegression()
model_adam = LinearRegression()

# ハイパーパラメータ
learning_rate = 0.01
n_epochs = 100
batch_size = 16
beta1 = 0.9
beta2 = 0.999
epsilon = 1e-8

# 損失の履歴を保存
loss_history_momentum = []
loss_history_rmsprop = []
loss_history_adam = []

# 各手法のトレーニング
for epoch in range(1, n_epochs + 1):
    # シャッフルしてミニバッチを作成
    indices = np.random.permutation(len(x_train))
    x_shuffled = x_train[indices]
    y_shuffled = y_train[indices]
    
    for start in range(0, len(x_train), batch_size):
        end = start + batch_size
        x_batch = x_shuffled[start:end]
        y_batch = y_shuffled[start:end]
        
        # 勾配の計算
        grad_weight_m, grad_bias_m = model_momentum.compute_gradients(x_batch, y_batch)
        grad_weight_r, grad_bias_r = model_rmsprop.compute_gradients(x_batch, y_batch)
        grad_weight_a, grad_bias_a = model_adam.compute_gradients(x_batch, y_batch)
        
        # 各手法のパラメータの更新
        model_momentum.update_parameters_momentum(grad_weight_m, grad_bias_m, learning_rate, beta=0.9)
        model_rmsprop.update_parameters_rmsprop(grad_weight_r, grad_bias_r, learning_rate, beta=0.9, epsilon=epsilon)
        model_adam.update_parameters_adam(grad_weight_a, grad_bias_a, learning_rate, beta1=beta1, beta2=beta2, epsilon=epsilon, t=epoch)
    
    # 全データでの損失を計算して記録
    loss_m = model_momentum.compute_loss(x_train, y_train)
    loss_r = model_rmsprop.compute_loss(x_train, y_train)
    loss_a = model_adam.compute_loss(x_train, y_train)
    loss_history_momentum.append(loss_m)
    loss_history_rmsprop.append(loss_r)
    loss_history_adam.append(loss_a)
    
    if epoch % 10 == 0:
        print(f'Epoch {epoch}, Loss (Momentum): {loss_m}, Loss (RMSProp): {loss_r}, Loss (Adam): {loss_a}')

# 学習結果の表示
plt.figure(figsize=(12, 5))

# 損失関数の履歴をプロット
plt.subplot(1, 2, 1)
plt.plot(range(n_epochs), loss_history_momentum, label='Momentum')
plt.plot(range(n_epochs), loss_history_rmsprop, label='RMSProp', linestyle='--')
plt.plot(range(n_epochs), loss_history_adam, label='Adam', linestyle='-.')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.title('Loss over Epochs')
plt.grid(True)
plt.legend()

# 学習結果のプロット（予測線）
plt.subplot(1, 2, 2)
plt.scatter(x_train, y_train, label='Data', color='blue', alpha=0.6)
plt.plot(x_train, model_momentum.predict(x_train), color='red', label='Momentum')
plt.plot(x_train, model_rmsprop.predict(x_train), color='green', linestyle='--', label='RMSProp')
plt.plot(x_train, model_adam.predict(x_train), color='orange', linestyle='-.', label='Adam')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Data and Learned Models')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# 学習されたパラメータを表示
print(f"Momentum - Learned weight: {model_momentum.weight}, bias: {model_momentum.bias}")
print(f"RMSProp - Learned weight: {model_rmsprop.weight}, bias: {model_rmsprop.bias}")
print(f"Adam - Learned weight: {model_adam.weight}, bias: {model_adam.bias}")



```


```
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 正規分布のパラメータ
mu_1, sigma_1 = 0, 1  # 分布 P: 平均0、標準偏差1
mu_2, sigma_2 = 1, 2  # 分布 Q: 平均1、標準偏差2

# KLダイバージェンスの計算
def kl_divergence_normal(mu1, sigma1, mu2, sigma2):
    return np.log(sigma2 / sigma1) + (sigma1 ** 2 + (mu1 - mu2) ** 2) / (2 * sigma2 ** 2) - 0.5

kl_div = kl_divergence_normal(mu_1, sigma_1, mu_2, sigma_2)
print(f"KL Divergence between N({mu_1}, {sigma_1**2}) and N({mu_2}, {sigma_2**2}): {kl_div}")

# 可視化
x = np.linspace(-5, 5, 1000)
p = norm.pdf(x, mu_1, sigma_1)  # 分布 P の確率密度関数
q = norm.pdf(x, mu_2, sigma_2)  # 分布 Q の確率密度関数

plt.plot(x, p, label=f'P: N({mu_1}, {sigma_1**2})', color='blue')
plt.plot(x, q, label=f'Q: N({mu_2}, {sigma_2**2})', color='red')
plt.fill_between(x, p, q, where=(p > q), color='blue', alpha=0.2, interpolate=True, label='P > Q')
plt.fill_between(x, p, q, where=(p < q), color='red', alpha=0.2, interpolate=True, label='Q > P')
plt.title('Probability Densities P and Q')
plt.xlabel('x')
plt.ylabel('Density')
plt.legend()
plt.grid(True)
plt.show()


```



```
import numpy as np
import matplotlib.pyplot as plt

# サンプルデータの生成
np.random.seed(42)
X = np.random.multivariate_normal([2, 3], [[1, 0.8], [0.8, 1]], size=100)

# 分散共分散行列の計算
cov_matrix = np.cov(X.T)  # 行ごとの共分散を計算するために転置
print("Covariance Matrix:")
print(cov_matrix)

# 分散共分散行列の固有値と固有ベクトルを求める（対角化）
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

print("\nEigenvalues:")
print(eigenvalues)
print("\nEigenvectors (Orthogonal Matrix):")
print(eigenvectors)

# 固有ベクトル行列（直交行列）の確認：Q^T * Q = I
Q = eigenvectors
Q_T_Q = Q.T @ Q
print("\nQ^T * Q:")
print(Q_T_Q)  # 単位行列に近ければ、直交性が確認できる

# データを新しい基底（固有ベクトル）に射影
X_projected = X @ Q

# 射影されたデータの分散共分散行列の対角成分のみが残ることを確認
projected_cov_matrix = np.cov(X_projected.T)
print("\nCovariance Matrix of Projected Data:")
print(projected_cov_matrix)

# 可視化
plt.figure(figsize=(8, 6))

# 元のデータのプロット
plt.scatter(X[:, 0], X[:, 1], alpha=0.5, label='Original Data', color='blue')

# 固有ベクトルの方向をプロット
origin = np.mean(X, axis=0)
for eigenvalue, eigenvector in zip(eigenvalues, eigenvectors.T):
    plt.quiver(*origin, *(eigenvector * 2), scale=5, color='red', label='Eigenvector', angles='xy', scale_units='xy')

# 射影されたデータのプロット
plt.scatter(X_projected[:, 0], X_projected[:, 1], alpha=0.5, label='Projected Data', color='green')

plt.title('Data Projection Using Eigenvectors (Principal Axes)')
plt.xlabel('X1')
plt.ylabel('X2')
plt.legend()
plt.grid(True)
plt.axis('equal')
plt.show()


```


```

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, hessian, Matrix, diff

# 変数の定義
x, y = symbols('x y')

# 関数の定義（例1: 極小点を持つ関数）
f = x**2 + y**2

# 関数の定義（例2: 鞍点を持つ関数）
g = x**2 - y**2

# ヘッセ行列の計算
H_f = hessian(f, (x, y))
H_g = hessian(g, (x, y))

print("Hessian Matrix of f(x, y) = x^2 + y^2:")
print(H_f)

print("\nHessian Matrix of g(x, y) = x^2 - y^2:")
print(H_g)

# ヘッセ行列を用いた解析
# ここでは (0, 0) を臨界点として解析
H_f_eval = H_f.subs({x: 0, y: 0})
H_g_eval = H_g.subs({x: 0, y: 0})

print("\nEvaluated Hessian at (0, 0) for f(x, y):")
print(H_f_eval)

print("\nEvaluated Hessian at (0, 0) for g(x, y):")
print(H_g_eval)

# 固有値を求める
eigenvalues_f = np.linalg.eigvals(np.array(H_f_eval).astype(np.float64))
eigenvalues_g = np.linalg.eigvals(np.array(H_g_eval).astype(np.float64))

print("\nEigenvalues of Hessian for f(x, y) at (0, 0):")
print(eigenvalues_f)

print("\nEigenvalues of Hessian for g(x, y) at (0, 0):")
print(eigenvalues_g)

# 極小、極大、鞍点の判定
def classify_critical_point(eigenvalues):
    if all(e > 0 for e in eigenvalues):
        return "Local minimum"
    elif all(e < 0 for e in eigenvalues):
        return "Local maximum"
    else:
        return "Saddle point"

print("\nClassification of critical point for f(x, y):")
print(classify_critical_point(eigenvalues_f))

print("\nClassification of critical point for g(x, y):")
print(classify_critical_point(eigenvalues_g))



```





```
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

# サンプルデータの生成（3次元データ）
np.random.seed(42)
mean = [0, 0, 0]
cov = [[1, 0.8, 0.5], [0.8, 1, 0.3], [0.5, 0.3, 1]]  # 共分散行列
data = np.random.multivariate_normal(mean, cov, 100)

# PCAの適用
pca = PCA(n_components=3)
pca.fit(data)

# 主成分の情報を取得
principal_components = pca.components_  # 主成分の方向（固有ベクトル）
explained_variance = pca.explained_variance_  # 各主成分の固有値（分散を表す）
explained_variance_ratio = pca.explained_variance_ratio_  # 各主成分の分散比率

print("Principal Components (Directions):")
print(principal_components)
print("\nExplained Variance (Eigenvalues):")
print(explained_variance)
print("\nExplained Variance Ratio:")
print(explained_variance_ratio)

# データを主成分に変換（射影）
data_transformed = pca.transform(data)

# 各主成分ごとの情報量の可視化
plt.figure(figsize=(8, 6))
plt.bar(range(1, len(explained_variance_ratio) + 1), explained_variance_ratio, alpha=0.7)
plt.xlabel('Principal Component')
plt.ylabel('Explained Variance Ratio')
plt.title('Explained Variance by Principal Components')
plt.grid(True)
plt.show()

# 2次元に次元削減したデータをプロット
plt.figure(figsize=(8, 6))
plt.scatter(data_transformed[:, 0], data_transformed[:, 1], alpha=0.6, color='blue', label='Projected Data')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.title('Data Projected onto First Two Principal Components')
plt.grid(True)
plt.legend()
plt.show()




```



```
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.cross_decomposition import CCA

# サンプルデータの生成（3次元データ）
np.random.seed(42)
mean = [0, 0, 0]
cov = [[1, 0.8, 0.5], [0.8, 1, 0.3], [0.5, 0.3, 1]]  # 共分散行列
data = np.random.multivariate_normal(mean, cov, 100)

# PCAの適用と分散共分散行列の表示
pca = PCA(n_components=2)
pca.fit(data)
data_pca = pca.transform(data)
cov_matrix = np.cov(data.T)

print("Covariance Matrix:")
print(cov_matrix)

print("\nPCA Components (Eigenvectors):")
print(pca.components_)

print("\nExplained Variance Ratio (PCA):")
print(pca.explained_variance_ratio_)

# 特異値分解（SVD）の適用
U, S, Vt = np.linalg.svd(data)
print("\nSingular Values (SVD):")
print(S)
print("\nLeft Singular Vectors (U):")
print(U[:5])  # 先頭5個のみ表示
print("\nRight Singular Vectors (Vt):")
print(Vt)

# CCAの適用（例として2次元データセットを生成して適用）
X = np.random.rand(100, 2)  # 2つのデータセットを用意
Y = X + np.random.normal(0, 0.1, X.shape)  # Xにノイズを加えたデータ
cca = CCA(n_components=2)
X_c, Y_c = cca.fit_transform(X, Y)

print("\nCanonical Correlations (CCA):")
print(np.corrcoef(X_c.T, Y_c.T))

# 可視化: PCA後のデータのプロット
plt.figure(figsize=(12, 5))

# PCAによる次元削減後のデータプロット
plt.subplot(1, 3, 1)
plt.scatter(data_pca[:, 0], data_pca[:, 1], alpha=0.7, label='PCA Projected Data')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('PCA Projection')
plt.grid(True)
plt.legend()

# SVDの特異値をプロット
plt.subplot(1, 3, 2)
plt.bar(range(1, len(S) + 1), S, alpha=0.7)
plt.xlabel('Component')
plt.ylabel('Singular Value')
plt.title('Singular Values (SVD)')
plt.grid(True)

# CCAによる相関のプロット
plt.subplot(1, 3, 3)
plt.scatter(X_c[:, 0], Y_c[:, 0], alpha=0.7, label='CCA Component 1')
plt.scatter(X_c[:, 1], Y_c[:, 1], alpha=0.7, label='CCA Component 2', color='orange')
plt.xlabel('X_c')
plt.ylabel('Y_c')
plt.title('Canonical Correlation Analysis')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()




```




```

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import TruncatedSVD
from sklearn.metrics import mean_squared_error

# ランダム行列の生成
np.random.seed(42)
A = np.random.randn(100, 50)  # 100x50のランダム行列

# 特異値分解（SVD）の計算
U, S, Vt = np.linalg.svd(A, full_matrices=False)
print("Singular values (first 10):", S[:10])

# 特異値の可視化
plt.figure(figsize=(8, 5))
plt.plot(S, 'bo-', label='Singular Values')
plt.xlabel('Index')
plt.ylabel('Singular Value')
plt.title('Singular Values of Random Matrix')
plt.grid(True)
plt.legend()
plt.show()

# 次元削減のためにTruncatedSVDを使用
svd = TruncatedSVD(n_components=10)
A_reduced = svd.fit_transform(A)  # 次元削減された行列（100x10）

# 次元削減後の再構成
A_reconstructed = svd.inverse_transform(A_reduced)

# 再構成誤差（損失関数としてのMSE）
mse_loss = mean_squared_error(A, A_reconstructed)
print(f"Reconstruction MSE Loss: {mse_loss}")

# データの生成（回帰モデルの例）
x = np.linspace(-1, 1, 100).reshape(-1, 1)
y = 3 * x.squeeze() + 2 + np.random.randn(100) * 0.1  # 線形データにノイズを追加

# 重み行列の生成と特異値解析（深層モデルを模擬）
W = np.random.randn(10, 1)
U_W, S_W, Vt_W = np.linalg.svd(W, full_matrices=False)
print("\nSingular values of weight matrix W:", S_W)

# 学習後の予測値（ここでは単純に線形変換を模擬）
y_pred = A_reduced @ W

# 損失関数の計算（MSE）
mse_pred = mean_squared_error(y, y_pred)
print(f"Prediction MSE Loss: {mse_pred}")

# 予測と実際の値の可視化
plt.figure(figsize=(8, 5))
plt.scatter(x, y, label='Actual', color='blue', alpha=0.5)
plt.scatter(x, y_pred, label='Predicted', color='red', alpha=0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Actual vs Predicted Values')
plt.legend()
plt.grid(True)
plt.show()


```




```
import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

# シンプルなテキストエンコーダー（BoWエンコーダーとして実装）
class TextEncoder(nn.Module):
    def __init__(self, vocab_size, embedding_dim):
        super(TextEncoder, self).__init__()
        self.embedding = nn.Embedding(vocab_size, embedding_dim)
        self.fc = nn.Linear(embedding_dim, 32)

    def forward(self, x):
        embedded = self.embedding(x).mean(dim=1)  # 各単語の埋め込みの平均を取る
        encoded = torch.relu(self.fc(embedded))
        return encoded

# 数値データエンコーダー
class NumericEncoder(nn.Module):
    def __init__(self, input_dim):
        super(NumericEncoder, self).__init__()
        self.fc = nn.Linear(input_dim, 32)

    def forward(self, x):
        encoded = torch.relu(self.fc(x))
        return encoded

# マルチモーダル統合モデル
class MultimodalModel(nn.Module):
    def __init__(self, vocab_size, embedding_dim, numeric_input_dim):
        super(MultimodalModel, self).__init__()
        self.text_encoder = TextEncoder(vocab_size, embedding_dim)
        self.numeric_encoder = NumericEncoder(numeric_input_dim)
        self.fc = nn.Linear(64, 2)  # テキストと数値の統合後の分類

    def forward(self, text_input, numeric_input):
        text_features = self.text_encoder(text_input)
        numeric_features = self.numeric_encoder(numeric_input)
        combined = torch.cat((text_features, numeric_features), dim=1)
        output = self.fc(combined)
        return output

# ダミーデータの生成
np.random.seed(42)
vocab_size = 100
embedding_dim = 50
numeric_input_dim = 10

# テキストデータ（語彙サイズが100までのインデックスのダミーデータ）
text_data = np.random.randint(0, vocab_size, (1000, 10))
numeric_data = np.random.randn(1000, numeric_input_dim)
labels = np.random.randint(0, 2, 1000)  # 2クラスのラベル

# データの前処理
scaler = StandardScaler()
numeric_data = scaler.fit_transform(numeric_data)
text_data = torch.tensor(text_data, dtype=torch.long)
numeric_data = torch.tensor(numeric_data, dtype=torch.float32)
labels = torch.tensor(labels, dtype=torch.long)

# データの分割
text_train, text_test, num_train, num_test, y_train, y_test = train_test_split(
    text_data, numeric_data, labels, test_size=0.2, random_state=42)

# モデルのインスタンス化
model = MultimodalModel(vocab_size, embedding_dim, numeric_input_dim)

# 損失関数とオプティマイザの定義
criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(model.parameters(), lr=0.001)

# トレーニングループ
n_epochs = 10
for epoch in range(n_epochs):
    model.train()
    optimizer.zero_grad()
    outputs = model(text_train, num_train)
    loss = criterion(outputs, y_train)
    loss.backward()
    optimizer.step()

    # 精度の計算
    model.eval()
    with torch.no_grad():
        test_outputs = model(text_test, num_test)
        test_preds = torch.argmax(test_outputs, dim=1)
        accuracy = accuracy_score(y_test.numpy(), test_preds.numpy())
    
    print(f'Epoch {epoch+1}/{n_epochs}, Loss: {loss.item()}, Test Accuracy: {accuracy}')

print("Training complete.")



```



```
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

# データの生成
np.random.seed(42)
n_samples = 100
X = np.random.rand(n_samples, 3)  # 3つの独立変数（X1, X2, X3）
beta_true = np.array([2, -3, 4])  # 真の回帰係数
y = X @ beta_true + np.random.randn(n_samples) * 0.5  # ノイズを加えた従属変数

# データの分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 最小二乗法によるパラメータの推定
X_train_b = np.c_[np.ones((X_train.shape[0], 1)), X_train]  # 切片用の1を追加
beta_hat = np.linalg.inv(X_train_b.T @ X_train_b) @ X_train_b.T @ y_train

print("Estimated coefficients (beta_hat):", beta_hat)

# 予測値の計算
X_test_b = np.c_[np.ones((X_test.shape[0], 1)), X_test]
y_pred = X_test_b @ beta_hat

# MSEの計算
mse = mean_squared_error(y_test, y_pred)
print("Mean Squared Error:", mse)

# Scikit-Learnによる重回帰の実装（簡易版）
model = LinearRegression()
model.fit(X_train, y_train)
y_pred_sklearn = model.predict(X_test)
mse_sklearn = mean_squared_error(y_test, y_pred_sklearn)

print("\nScikit-Learn Estimated coefficients (including intercept):", np.append(model.intercept_, model.coef_))
print("Scikit-Learn Mean Squared Error:", mse_sklearn)

# 結果のプロット（予測値と実際の値の比較）
plt.figure(figsize=(8, 5))
plt.scatter(y_test, y_pred, label='Manual OLS Prediction', color='blue', alpha=0.6)
plt.scatter(y_test, y_pred_sklearn, label='Scikit-Learn Prediction', color='red', alpha=0.6)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], color='green', linestyle='--', label='Ideal Fit')
plt.xlabel('Actual values')
plt.ylabel('Predicted values')
plt.title('Actual vs Predicted')
plt.legend()
plt.grid(True)
plt.show()




```




```
import numpy as np
import matplotlib.pyplot as plt

# データの生成
np.random.seed(42)
n_samples = 100
n_features = 3
X = np.random.randn(n_samples, n_features)  # 説明変数（デザイン行列）
beta_true = np.array([2, -3, 4])  # 真の回帰係数
y = X @ beta_true + np.random.randn(n_samples) * 0.5  # ノイズを加えた目的変数

# 擬似逆行列を用いた最小二乗推定量の計算
X_pseudo_inv = np.linalg.pinv(X)  # 擬似逆行列の計算
beta_hat = X_pseudo_inv @ y  # 最小二乗推定量の計算

print("True coefficients (beta_true):", beta_true)
print("Estimated coefficients (beta_hat):", beta_hat)

# 予測値の計算
y_pred = X @ beta_hat

# 予測値と実際の値のプロット
plt.figure(figsize=(8, 5))
plt.scatter(y, y_pred, alpha=0.7, color='blue', label='Predicted vs Actual')
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', label='Ideal Fit')
plt.xlabel('Actual values')
plt.ylabel('Predicted values')
plt.title('Predicted vs Actual Values using Pseudo-Inverse')
plt.grid(True)
plt.legend()
plt.show()

# MSE（平均二乗誤差）の計算
mse = np.mean((y - y_pred) ** 2)
print("Mean Squared Error (MSE):", mse)




```





```


import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# データ生成（重回帰分析用）
np.random.seed(42)
n_samples = 100
X = np.random.randn(n_samples, 3)  # 3つの独立変数
true_beta = np.array([2, -1, 3])  # 真の回帰係数
y = X @ true_beta + np.random.randn(n_samples) * 0.5  # ノイズ付きの従属変数

# SVDを使って擬似逆行列を計算し、重回帰分析の推定量を求める
U, S, Vt = np.linalg.svd(X, full_matrices=False)
Sigma_pinv = np.diag(1 / S)  # Σ^+ を計算
X_pinv = Vt.T @ Sigma_pinv @ U.T  # 擬似逆行列 X^+ = V Σ^+ U^T
beta_hat = X_pinv @ y

print("Estimated coefficients (beta_hat using SVD):", beta_hat)

# NumPyを用いた最小二乗推定量との比較
beta_hat_np = np.linalg.pinv(X) @ y
print("\nEstimated coefficients (beta_hat using NumPy's pinv):", beta_hat_np)

# 多変量正規分布の生成
mean = np.array([0, 0])  # 平均ベクトル
cov = np.array([[1, 0.5], [0.5, 1]])  # 共分散行列
data = np.random.multivariate_normal(mean, cov, 500)

# 多変量正規分布のプロット
plt.figure(figsize=(6, 6))
plt.scatter(data[:, 0], data[:, 1], alpha=0.5, label='Samples from Multivariate Normal')
plt.xlabel('X1')
plt.ylabel('X2')
plt.title('Samples from Multivariate Normal Distribution')
plt.grid(True)
plt.axis('equal')
plt.legend()
plt.show()

# 変分自由エネルギーの計算
# ここでは、簡単な例として、データ点の分布と多変量正規分布のKLダイバージェンスを計算します
mu_q = np.mean(data, axis=0)  # q(x) の平均（データの平均）
cov_q = np.cov(data.T)  # q(x) の共分散行列（データの共分散）

# p(x) は与えられた平均と共分散を持つ分布
p = multivariate_normal(mean=mean, cov=cov)

# 変分自由エネルギーとして、データ点に対する対数尤度を計算（簡単な例）
log_likelihood = np.mean(p.logpdf(data))
kl_divergence = np.trace(np.linalg.inv(cov) @ cov_q) + (mu_q - mean).T @ np.linalg.inv(cov) @ (mu_q - mean) - data.shape[1] + np.log(np.linalg.det(cov) / np.linalg.det(cov_q))

print("\nLog-Likelihood of data under p(x):", log_likelihood)
print("KL Divergence between q(x) and p(x):", kl_divergence)

# 変分自由エネルギーの近似として、ELBO (Evidence Lower Bound) を計算
elbo = log_likelihood - kl_divergence
print("\nApproximate Variational Free Energy (ELBO):", elbo)


```





```




```





```




```