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

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)