データ分析のための数理モデル入門 本質をとらえた分析のためにを読んだあとのPythonコード

import sympy as sp

# Define the variable and the cost function
P = sp.Symbol('P')
C = 2000 / P**0.3 + 0.05 * P**2 + 500

# Calculate the derivative of the cost function
dC_dP = sp.diff(C, P)

# Solve for dC/dP = 0 to find the critical points
critical_points = sp.solve(dC_dP, P)

# Display the critical points and evaluate the cost at these points
critical_points = [point.evalf() for point in critical_points if point.is_real and point > 0]
optimal_costs = [(P_val, C.subs(P, P_val).evalf()) for P_val in critical_points]

optimal_costs

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

# Set random seed for reproducibility
np.random.seed(0)

# Parameters for normal distribution
mu = 100    # Mean
sigma = 15  # Standard deviation

# Generate 10,000 samples from the normal distribution
samples = np.random.normal(mu, sigma, 10000)

# Plot the histogram of the samples
plt.figure(figsize=(10, 6))
plt.hist(samples, bins=50, density=True, alpha=0.6, color='g')

# Plot the probability density function (PDF) for comparison
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = stats.norm.pdf(x, mu, sigma)
plt.plot(x, p, 'k', linewidth=2)
title = f"Normal Distribution (mean = {mu}, std = {sigma})"
plt.title(title)
plt.xlabel('Value')
plt.ylabel('Probability Density')
plt.grid()
plt.show()

# Central Limit Theorem: Sum of 10 normal distributions
sums = np.sum(np.random.normal(mu, sigma, (10000, 10)), axis=1)

# Plot the histogram of the summed samples
plt.figure(figsize=(10, 6))
plt.hist(sums, bins=50, density=True, alpha=0.6, color='b')

# Plot the PDF of the resulting normal distribution
mu_sum = np.mean(sums)
sigma_sum = np.std(sums)
p_sum = stats.norm.pdf(x, mu_sum, sigma_sum)
plt.plot(x, p_sum, 'k', linewidth=2)
title = f"Sum of 10 Independent Normal Distributions (mean = {mu_sum:.2f}, std = {sigma_sum:.2f})"
plt.title(title)
plt.xlabel('Value')
plt.ylabel('Probability Density')
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from scipy.stats import pearsonr

# Create synthetic datasets for demonstration
# Dataset 1
np.random.seed(0)
X1 = np.random.rand(50, 1)
Y1 = 0.5 * X1 + np.random.normal(0, 0.05, (50, 1))  # Linear relationship with some noise

# Dataset 2
X2 = np.random.rand(50, 1)
Y2 = 0.5 * X2 + np.random.normal(0, 0.15, (50, 1))  # Similar relationship but more noise

# Fit linear regression models
model1 = LinearRegression().fit(X1, Y1)
model2 = LinearRegression().fit(X2, Y2)

# Predict using the models
Y1_pred = model1.predict(X1)
Y2_pred = model2.predict(X2)

# Calculate R-squared values
r2_1 = r2_score(Y1, Y1_pred)
r2_2 = r2_score(Y2, Y2_pred)

# Calculate Pearson correlation coefficients
pearson_r1, _ = pearsonr(X1.flatten(), Y1.flatten())
pearson_r2, _ = pearsonr(X2.flatten(), Y2.flatten())

# Print results
print(f"Dataset 1: R^2 = {r2_1:.2f}, Pearson r = {pearson_r1:.2f}")
print(f"Dataset 2: R^2 = {r2_2:.2f}, Pearson r = {pearson_r2:.2f}")

# Plot the datasets and regression lines
plt.figure(figsize=(14, 6))

# Plot Dataset 1
plt.subplot(1, 2, 1)
plt.scatter(X1, Y1, color='blue', label='Data Points')
plt.plot(X1, Y1_pred, color='red', label=f'Linear Fit (R^2 = {r2_1:.2f}, r = {pearson_r1:.2f})')
plt.title("Dataset 1: Small Variance")
plt.xlabel("X")
plt.ylabel("Y")
plt.legend()
plt.grid()

# Plot Dataset 2
plt.subplot(1, 2, 2)
plt.scatter(X2, Y2, color='blue', label='Data Points')
plt.plot(X2, Y2_pred, color='red', label=f'Linear Fit (R^2 = {r2_2:.2f}, r = {pearson_r2:.2f})')
plt.title("Dataset 2: Large Variance")
plt.xlabel("X")
plt.ylabel("Y")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose

# Set random seed for reproducibility
np.random.seed(42)

# Create a date range
date_range = pd.date_range(start='2020-01-01', periods=200, freq='D')

# Generate trend component
trend = np.linspace(start=10, stop=50, num=200)

# Generate seasonal component (sinusoidal wave for seasonality)
seasonal = 10 * np.sin(np.linspace(0, 3 * np.pi, 200))

# Generate noise component
noise = np.random.normal(0, 3, 200)

# Combine the components to form the time series data
time_series = trend + seasonal + noise

# Create a pandas DataFrame
data = pd.DataFrame({'Date': date_range, 'Value': time_series})
data.set_index('Date', inplace=True)

# Plot the original time series
plt.figure(figsize=(12, 6))
plt.plot(data, label='Original Time Series')
plt.title('Original Time Series (Trend + Seasonality + Noise)')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.grid()
plt.show()

# Decompose the time series into trend, seasonality, and noise
decomposition = seasonal_decompose(data['Value'], model='additive', period=30)  # Period can be adjusted based on seasonality

# Plot the decomposition
fig = decomposition.plot()
fig.set_size_inches(12, 10)
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller

# サンプルデータの作成（非定常な時系列データ）
np.random.seed(0)
date_range = pd.date_range(start='2020-01-01', periods=200, freq='D')
trend = np.linspace(10, 50, 200)  # 明確なトレンド成分
seasonal = 10 * np.sin(np.linspace(0, 4 * np.pi, 200))  # 季節性成分
noise = np.random.normal(0, 3, 200)  # ノイズ成分
data = trend + seasonal + noise

# データフレームの作成
df = pd.DataFrame({'Value': data}, index=date_range)

# オリジナルの時系列データをプロット
plt.figure(figsize=(12, 6))
plt.plot(df, label='Original Time Series')
plt.title('Original Time Series with Trend and Seasonality')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.grid()
plt.show()

# 1. トレンドを取り除く（1次差分を取る）
diff_data = df.diff().dropna()

plt.figure(figsize=(12, 6))
plt.plot(diff_data, label='1st Order Differenced Data')
plt.title('1st Order Differenced Time Series')
plt.xlabel('Date')
plt.ylabel('Differenced Value')
plt.legend()
plt.grid()
plt.show()

# 2. 季節性を取り除く（周期30日と仮定して差分を取る）
seasonal_diff_data = df.diff(30).dropna()

plt.figure(figsize=(12, 6))
plt.plot(seasonal_diff_data, label='Seasonally Differenced Data (Lag=30)')
plt.title('Seasonally Differenced Time Series')
plt.xlabel('Date')
plt.ylabel('Seasonally Differenced Value')
plt.legend()
plt.grid()
plt.show()

# 3. 対数変換（分散の安定化）
log_data = np.log(df)

plt.figure(figsize=(12, 6))
plt.plot(log_data, label='Log Transformed Data')
plt.title('Log Transformed Time Series')
plt.xlabel('Date')
plt.ylabel('Log Value')
plt.legend()
plt.grid()
plt.show()

# 4. 定常性の確認（ADF検定）
# ADF検定の関数
def adf_test(timeseries):
    result = adfuller(timeseries)
    print('ADF Statistic: %f' % result[0])
    print('p-value: %f' % result[1])
    print('Critical Values:')
    for key, value in result[4].items():
        print('\t%s: %.3f' % (key, value))

print("Original Data ADF Test:")
adf_test(df['Value'])

print("\nDifferenced Data ADF Test:")
adf_test(diff_data['Value'])

print("\nSeasonally Differenced Data ADF Test:")
adf_test(seasonal_diff_data['Value'])

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# サンプルデータの生成（非定常な時系列データ）
np.random.seed(42)
date_range = pd.date_range(start='2020-01-01', periods=200, freq='D')
trend = np.linspace(10, 50, 200)  # 緩やかなトレンド成分
seasonal = 10 * np.sin(np.linspace(0, 4 * np.pi, 200))  # 季節性成分
noise = np.random.normal(0, 3, 200)  # ノイズ成分
data = trend + seasonal + noise

# データフレームの作成
df = pd.DataFrame({'Value': data}, index=date_range)

# 1. 元データのプロット
plt.figure(figsize=(12, 6))
plt.plot(df, label='Original Time Series')
plt.title('Original Time Series with Trend and Seasonality')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.grid()
plt.show()

# 2. ARIMA モデルのフィッティングと予測
# 非定常性を取り除く（1次差分を取る）
diff_data = df.diff().dropna()

# ACFおよびPACFプロット
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
plot_acf(diff_data, ax=ax1, title='ACF (Autocorrelation)')
plot_pacf(diff_data, ax=ax2, title='PACF (Partial Autocorrelation)')
plt.show()

# ARIMAモデルの構築とフィッティング
arima_model = ARIMA(df['Value'], order=(1, 1, 1))  # (p, d, q)
arima_result = arima_model.fit()

# 予測
arima_forecast = arima_result.forecast(steps=30)  # 30日間の予測

# 3. SARIMA モデルのフィッティングと予測
# SARIMAモデルの構築とフィッティング（p, d, q, x [P, D, Q, m]）
sarima_model = SARIMAX(df['Value'], order=(1, 1, 1), seasonal_order=(1, 1, 1, 30))
sarima_result = sarima_model.fit()

# 予測
sarima_forecast = sarima_result.forecast(steps=30)  # 30日間の予測

# 4. 予測結果の可視化
plt.figure(figsize=(12, 6))
plt.plot(df, label='Original Time Series')
plt.plot(arima_forecast, label='ARIMA Forecast', color='orange')
plt.plot(sarima_forecast, label='SARIMA Forecast', color='green')
plt.title('ARIMA and SARIMA Forecast')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.grid()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from filterpy.kalman import KalmanFilter

# サンプルデータの作成
np.random.seed(42)
n = 100  # データ数
dt = 1  # 時間ステップ

# 状態変数 (位置と速度)
x = np.zeros((2, n))
x[0, 0] = 0  # 初期位置
x[1, 0] = 1  # 初期速度

# 状態遷移行列
G = np.array([[1, dt],
              [0, 1]])  # 位置と速度の関係を表す

# 観測行列
F = np.array([[1, 0]])  # 観測は位置のみ

# 状態ノイズと観測ノイズの標準偏差
sigma_w = 0.1  # 状態ノイズの標準偏差
sigma_v = 1.0  # 観測ノイズの標準偏差

# 状態ノイズと観測ノイズの生成
w = np.random.normal(0, sigma_w, (2, n))  # 状態ノイズ
v = np.random.normal(0, sigma_v, n)  # 観測ノイズ

# 状態の生成
for t in range(1, n):
    x[:, t] = G @ x[:, t-1] + w[:, t]

# 観測データの生成
y = F @ x + v

# 観測データのプロット
plt.figure(figsize=(12, 6))
plt.plot(range(n), y[0], 'ro', markersize=3, label='Observed Data')
plt.title('Generated Observed Data with Noise')
plt.xlabel('Time')
plt.ylabel('Position')
plt.legend()
plt.grid()
plt.show()

# カルマンフィルターの設定
kf = KalmanFilter(dim_x=2, dim_z=1)
kf.F = G  # 状態遷移行列
kf.H = F  # 観測行列
kf.R = np.array([[sigma_v**2]])  # 観測ノイズの分散行列
kf.Q = np.array([[sigma_w**2, 0], [0, sigma_w**2]])  # 状態ノイズの分散行列
kf.x = np.array([0, 0])  # 状態の初期値
kf.P = np.eye(2)  # 初期共分散行列

# カルマンフィルターによる状態推定
filtered_state_means = np.zeros((2, n))  # フィルタされた状態の保存用
for t in range(n):
    kf.predict()  # 予測
    kf.update(y[:, t])  # 観測データで更新
    filtered_state_means[:, t] = kf.x  # フィルタされた状態を保存

# 状態推定結果のプロット
plt.figure(figsize=(12, 6))
plt.plot(range(n), x[0], 'b-', label='True Position')  # 真の位置
plt.plot(range(n), y[0], 'ro', markersize=3, label='Observed Data')  # 観測データ
plt.plot(range(n), filtered_state_means[0], 'g-', label='Filtered Position')  # フィルタされた位置
plt.title('Kalman Filter State Estimation')
plt.xlabel('Time')
plt.ylabel('Position')
plt.legend()
plt.grid()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf

# サンプル時系列データの生成
np.random.seed(42)
n = 200
t = np.arange(n)
# 周期成分とノイズを含む時系列データ
y = 10 * np.sin(2 * np.pi * t / 24) + np.random.normal(0, 2, n)

# データフレームに変換
df = pd.DataFrame({'Value': y}, index=pd.date_range(start='2020-01-01', periods=n, freq='D'))

# 自己相関のプロット
plt.figure(figsize=(12, 6))
plt.plot(df, label='Time Series Data')
plt.title('Time Series Data')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.grid()
plt.show()

# 自己相関のプロット
plot_acf(df['Value'], lags=50)
plt.title('Autocorrelation Function (ACF)')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import dct
from scipy.signal import stft, get_window
import pywt  # ウェーブレット変換用のライブラリ

# サンプル信号の生成
t = np.linspace(0, 1, 500, endpoint=False)  # 500サンプルの時間軸
signal = np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 80 * t)  # 周波数50Hzと80Hzの正弦波

# 1. 離散フーリエ変換 (DFT)
dft = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(signal), t[1] - t[0])

# 2. 高速フーリエ変換 (FFT) -> FFTとDFTは同じ関数を使う
fft = dft

# 3. 窓付きフーリエ変換 (Hann Window)
window = get_window('hann', len(signal))
windowed_signal = signal * window
fft_windowed = np.fft.fft(windowed_signal)

# 4. 短時間フーリエ変換 (STFT)
f_stft, t_stft, Zxx = stft(signal, fs=1.0, nperseg=100)

# 5. 離散コサイン変換 (DCT)
dct_transformed = dct(signal, type=2, norm='ortho')

# 6. ウェーブレット変換 (Wavelet Transform)
wavelet_coeffs, freqs_wt = pywt.cwt(signal, scales=np.arange(1, 128), wavelet='morl', sampling_period=t[1] - t[0])

# プロット
fig, axs = plt.subplots(3, 2, figsize=(12, 12))

# 1. DFT
axs[0, 0].plot(frequencies, np.abs(dft))
axs[0, 0].set_title("Discrete Fourier Transform (DFT)")
axs[0, 0].set_xlabel("Frequency (Hz)")
axs[0, 0].set_ylabel("Magnitude")
axs[0, 0].set_xlim(0, 150)

# 2. FFT
axs[0, 1].plot(frequencies, np.abs(fft))
axs[0, 1].set_title("Fast Fourier Transform (FFT)")
axs[0, 1].set_xlabel("Frequency (Hz)")
axs[0, 1].set_ylabel("Magnitude")
axs[0, 1].set_xlim(0, 150)

# 3. 窓付きフーリエ変換
axs[1, 0].plot(frequencies, np.abs(fft_windowed))
axs[1, 0].set_title("Windowed Fourier Transform (Hann Window)")
axs[1, 0].set_xlabel("Frequency (Hz)")
axs[1, 0].set_ylabel("Magnitude")
axs[1, 0].set_xlim(0, 150)

# 4. 短時間フーリエ変換 (STFT)
axs[1, 1].pcolormesh(t_stft, f_stft, np.abs(Zxx), shading='gouraud')
axs[1, 1].set_title("Short-Time Fourier Transform (STFT)")
axs[1, 1].set_xlabel("Time (s)")
axs[1, 1].set_ylabel("Frequency (Hz)")
axs[1, 1].set_ylim(0, 150)

# 5. 離散コサイン変換 (DCT)
axs[2, 0].plot(dct_transformed)
axs[2, 0].set_title("Discrete Cosine Transform (DCT)")
axs[2, 0].set_xlabel("Index")
axs[2, 0].set_ylabel("Magnitude")

# 6. ウェーブレット変換
axs[2, 1].imshow(np.abs(wavelet_coeffs), extent=[0, 1, 1, 128], cmap='PRGn', aspect='auto',
                 vmax=np.abs(wavelet_coeffs).max(), vmin=-np.abs(wavelet_coeffs).max())
axs[2, 1].set_title("Wavelet Transform (CWT)")
axs[2, 1].set_xlabel("Time (s)")
axs[2, 1].set_ylabel("Scale")

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft, fftshift
from scipy.signal import hilbert
from scipy.fft import fft2, fftn, ifftn
from scipy.linalg import expm
import cmath

# サンプル信号を作成
t = np.linspace(0, 1, 400)
f1 = 5  # 5 Hz
f2 = 20  # 20 Hz
signal = np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)

# 1. 標準のフーリエ変換と逆フーリエ変換
fft_signal = fft(signal)
ifft_signal = ifft(fft_signal)

# 2. ハルフーリエ変換 (Hartley Transform)
hartley_transform = np.real(fft_signal) - np.imag(fft_signal)

# 3. 2D フーリエ変換
# サンプル2Dデータ作成
X, Y = np.meshgrid(t, t)
Z = np.sin(2 * np.pi * f1 * X) + np.sin(2 * np.pi * f2 * Y)
fft2_signal = fft2(Z)

# 4. フェーズボコーダー
analytic_signal = hilbert(signal)
amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))

# プロット
fig, axs = plt.subplots(4, 2, figsize=(15, 10))

# 1. 標準のフーリエ変換と逆フーリエ変換
axs[0, 0].plot(t, signal)
axs[0, 0].set_title('Original Signal')
axs[0, 1].plot(t, np.abs(fft_signal))
axs[0, 1].set_title('FFT of Signal')

# 2. ハルフーリエ変換
axs[1, 0].plot(t, hartley_transform)
axs[1, 0].set_title('Hartley Transform of Signal')

# 3. 2Dフーリエ変換
axs[1, 1].imshow(np.log(np.abs(fftshift(fft2_signal))), extent=[0, 1, 0, 1], cmap='jet')
axs[1, 1].set_title('2D Fourier Transform of 2D Signal')

# 4. フェーズボコーダー
axs[2, 0].plot(t, amplitude_envelope)
axs[2, 0].set_title('Amplitude Envelope of Signal')
axs[2, 1].plot(t, instantaneous_phase)
axs[2, 1].set_title('Instantaneous Phase of Signal')

# 5. ラプラス変換とゼータ変換は、一般的にシンボリックな解や信号の解析に使うため、簡単なプロットではなく、解析ツールが必要です。
axs[3, 0].text(0.5, 0.5, "Laplace and Z-Transform\nVisual representation not supported", fontsize=14, ha='center')
axs[3, 0].set_axis_off()
axs[3, 1].text(0.5, 0.5, "Laplace and Z-Transform\nVisual representation not supported", fontsize=14, ha='center')
axs[3, 1].set_axis_off()

plt.tight_layout()
plt.show()

# 必要なライブラリをインポート
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score

# データセットの読み込み
iris = load_iris()
X = iris.data  # 特徴量
y = iris.target  # ラベル

# 訓練データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# ロジスティック回帰モデルの作成
model = LogisticRegression(max_iter=200)
model.fit(X_train, y_train)

# テストデータでの予測
y_pred = model.predict(X_test)

# 精度の計算
accuracy = accuracy_score(y_test, y_pred)
print(f'Accuracy: {accuracy:.2f}')

# 必要なライブラリをインポート
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

# データセットの読み込み
boston = load_boston()
X = boston.data  # 特徴量
y = boston.target  # ラベル（住宅価格）

# 訓練データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 線形回帰モデルの作成
model = LinearRegression()
model.fit(X_train, y_train)

# テストデータでの予測
y_pred = model.predict(X_test)

# 平均二乗誤差 (MSE) の計算
mse = mean_squared_error(y_test, y_pred)
print(f'Mean Squared Error: {mse:.2f}')

# 必要なライブラリをインポート
from sklearn.datasets import load_iris
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt

# データセットの読み込み
iris = load_iris()
X = iris.data  # 特徴量のみを使用（ラベルは使用しない）

# K-means モデルの作成（クラスタ数を3に指定）
kmeans = KMeans(n_clusters=3, random_state=42)
kmeans.fit(X)

# 各データポイントのクラスタラベルを取得
labels = kmeans.labels_

# クラスタリング結果を可視化
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis')
plt.xlabel(iris.feature_names[0])
plt.ylabel(iris.feature_names[1])
plt.title('K-means Clustering of Iris Dataset')
plt.show()

import numpy as np
from scipy.linalg import expm
import matplotlib.pyplot as plt

# 係数行列Aを定義
A = np.array([[1, 2], [3, 4]])

# 固有値と固有ベクトルを求める
eigenvalues, eigenvectors = np.linalg.eig(A)

# 固有値と固有ベクトルを表示
print("Eigenvalues:", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

# 初期条件ベクトルx0
x0 = np.array([1, 0])

# 行列指数関数を使って時間tでの解を求める関数
def solve_system(A, x0, t):
    exp_At = expm(A * t)  # 行列の指数関数
    return np.dot(exp_At, x0)

# tの範囲を指定
t_values = np.linspace(0, 5, 100)

# 解を格納するリスト
x1_values = []
x2_values = []

# 各時刻tに対する解を計算
for t in t_values:
    x_t = solve_system(A, x0, t)
    x1_values.append(x_t[0])
    x2_values.append(x_t[1])

# 結果をプロット
plt.plot(t_values, x1_values, label="x1(t)")
plt.plot(t_values, x2_values, label="x2(t)")
plt.xlabel('Time t')
plt.ylabel('Solution x(t)')
plt.legend()
plt.title("Solution of the System of Differential Equations")
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

# 状態遷移行列を定義
P = np.array([[0.1, 0.6, 0.3], 
              [0.4, 0.3, 0.3], 
              [0.2, 0.5, 0.3]])

# ノード（状態）の定義
states = ['State 1', 'State 2', 'State 3']

# NetworkX グラフを作成
G = nx.DiGraph()  # 有向グラフ

# ノードを追加
for state in states:
    G.add_node(state)

# エッジ（遷移）とその重み（遷移確率）を追加
for i, origin in enumerate(states):
    for j, destination in enumerate(states):
        probability = P[i, j]
        if probability > 0:  # 遷移確率が0でなければエッジを追加
            G.add_edge(origin, destination, weight=probability)

# レイアウトの設定
pos = nx.spring_layout(G)

# ノードとエッジを描画
nx.draw_networkx_nodes(G, pos, node_size=2000, node_color='lightblue')
nx.draw_networkx_labels(G, pos, font_size=12, font_weight='bold')

# エッジを描画（遷移確率をエッジのラベルとして表示）
edge_labels = {(origin, destination): f'{G[origin][destination]["weight"]:.2f}' 
               for origin, destination in G.edges}
nx.draw_networkx_edges(G, pos, arrowstyle='->', arrowsize=15)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=10)

# グラフの表示
plt.title('State Transition Diagram')
plt.axis('off')  # 軸を非表示にする
plt.show()

import numpy as np

# 1. 状態遷移行列を定義
P = np.array([[0.5, 0.3, 0.2], 
              [0.2, 0.6, 0.2], 
              [0.1, 0.3, 0.6]])

# 2. 定常状態の確率分布を求めるための方程式
# 定常状態πは πP = π を満たすので、転置行列を使って方程式を設定する
# πP - π = 0 ⇔ (P^T - I)π = 0, ただしπの和は1（∑π = 1）という制約を追加
I = np.eye(len(P))  # 単位行列
A = (P.T - I)[:-1]  # 最後の行を削除しランクを保つ
b = np.zeros(len(P)-1)

# 制約条件 π1 + π2 + π3 = 1 を追加
A = np.vstack([A, np.ones(len(P))])
b = np.append(b, 1)

# 3. 定常状態の確率分布を解く（連立方程式を解く）
steady_state = np.linalg.solve(A, b)
print("定常状態の確率分布: ", steady_state)

# 4. 各状態の稼働率を計算（ここでは稼働率 = 定常状態の確率分布）
# 状態の稼働率を表示
for i, rate in enumerate(steady_state):
    print(f"状態 {i + 1} の稼働率: {rate:.2f}")

# 5. 定常状態の可視化（棒グラフ）
import matplotlib.pyplot as plt

plt.bar(range(1, len(steady_state)+1), steady_state, color='skyblue')
plt.xlabel('State')
plt.ylabel('Steady State Probability')
plt.title('Steady State Probability Distribution')
plt.xticks(range(1, len(steady_state)+1))
plt.grid(axis='y')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from tensorflow.keras.datasets import mnist
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Dense

# 1. データセットの読み込みと前処理
(x_train, _), (x_test, _) = mnist.load_data()  # ラベルは不要なので "_" を使用
x_train = x_train.astype('float32') / 255.0  # 正規化
x_test = x_test.astype('float32') / 255.0    # 正規化
x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:])))  # (60000, 784)
x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:])))      # (10000, 784)

# 2. オートエンコーダーのモデル構築
# エンコーダー部分
input_img = Input(shape=(784,))  # 28x28ピクセル = 784次元の入力
encoded = Dense(64, activation='relu')(input_img)  # エンコーディング: 64次元に圧縮

# デコーダー部分
decoded = Dense(784, activation='sigmoid')(encoded)  # 再構築: 784次元に戻す

# 3. オートエンコーダーモデルの定義
autoencoder = Model(input_img, decoded)  # 入力から出力までのモデルを定義

# 4. モデルのコンパイル
autoencoder.compile(optimizer='adam', loss='binary_crossentropy')

# 5. モデルの訓練
autoencoder.fit(x_train, x_train,
                epochs=50,
                batch_size=256,
                shuffle=True,
                validation_data=(x_test, x_test))  # 入力と出力が同じ（自己符号化）

# 6. テストデータを用いた復元の確認
decoded_imgs = autoencoder.predict(x_test)

# 7. 元の画像と復元された画像を表示
n = 10  # 表示する画像の数
plt.figure(figsize=(20, 4))
for i in range(n):
    # 元の画像
    ax = plt.subplot(2, n, i + 1)
    plt.imshow(x_test[i].reshape(28, 28), cmap='gray')
    plt.title("Original")
    plt.axis('off')

    # 復元された画像
    ax = plt.subplot(2, n, i + 1 + n)
    plt.imshow(decoded_imgs[i].reshape(28, 28), cmap='gray')
    plt.title("Reconstructed")
    plt.axis('off')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from tensorflow.keras.datasets import mnist
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Dense

# 1. データセットの読み込みと前処理
(x_train, _), (x_test, _) = mnist.load_data()  # ラベルは不要なので "_" を使用
x_train = x_train.astype('float32') / 255.0  # 正規化
x_test = x_test.astype('float32') / 255.0    # 正規化
x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:])))  # (60000, 784)
x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:])))      # (10000, 784)

# 2. オートエンコーダーのモデル構築
# エンコーダー部分
input_img = Input(shape=(784,))  # 28x28ピクセル = 784次元の入力
encoded = Dense(64, activation='relu')(input_img)  # エンコーディング: 64次元に圧縮

# デコーダー部分
decoded = Dense(784, activation='sigmoid')(encoded)  # 再構築: 784次元に戻す

# 3. オートエンコーダーモデルの定義
autoencoder = Model(input_img, decoded)  # 入力から出力までのモデルを定義

# 4. モデルのコンパイル
autoencoder.compile(optimizer='adam', loss='binary_crossentropy')

# 5. モデルの訓練
autoencoder.fit(x_train, x_train,
                epochs=50,
                batch_size=256,
                shuffle=True,
                validation_data=(x_test, x_test))  # 入力と出力が同じ（自己符号化）

# 6. テストデータを用いた復元の確認
decoded_imgs = autoencoder.predict(x_test)

# 7. 元の画像と復元された画像を表示
n = 10  # 表示する画像の数
plt.figure(figsize=(20, 4))
for i in range(n):
    # 元の画像
    ax = plt.subplot(2, n, i + 1)
    plt.imshow(x_test[i].reshape(28, 28), cmap='gray')
    plt.title("Original")
    plt.axis('off')

    # 復元された画像
    ax = plt.subplot(2, n, i + 1 + n)
    plt.imshow(decoded_imgs[i].reshape(28, 28), cmap='gray')
    plt.title("Reconstructed")
    plt.axis('off')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture

# 1. 多次元時系列データの生成
# 例として、2つの異なる正規分布から生成された2次元時系列データを作成
np.random.seed(42)
n_samples = 500

# クラスタ1のデータ (平均 [0, 0], 分散 [1, 1])
mean1 = [0, 0]
cov1 = [[1, 0.5], [0.5, 1]]  # 共分散行列
data1 = np.random.multivariate_normal(mean1, cov1, n_samples)

# クラスタ2のデータ (平均 [3, 3], 分散 [1, 1])
mean2 = [3, 3]
cov2 = [[1, -0.7], [-0.7, 1]]
data2 = np.random.multivariate_normal(mean2, cov2, n_samples)

# データを結合して、ラベルを作成（実際のクラスタラベル）
X = np.vstack((data1, data2))
true_labels = np.array([0] * n_samples + [1] * n_samples)

# 2. ガウス混合モデル（GMM）の定義とフィッティング
gmm = GaussianMixture(n_components=2, covariance_type='full', random_state=42)
gmm.fit(X)

# 3. データポイントのクラスタリング（所属するクラスタの推定）
predicted_labels = gmm.predict(X)
cluster_centers = gmm.means_  # 各クラスタの平均値（中心）

# 4. 結果のプロット（元のラベルとGMMによるクラスタリング結果）
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# 元のラベルによるプロット
ax1.scatter(X[:, 0], X[:, 1], c=true_labels, cmap='viridis', s=10)
ax1.scatter(mean1[0], mean1[1], c='red', s=200, marker='x', label='True Center 1')
ax1.scatter(mean2[0], mean2[1], c='blue', s=200, marker='x', label='True Center 2')
ax1.set_title("Original Data with True Labels")
ax1.legend()
ax1.grid(True)

# GMMによるクラスタリング結果のプロット
ax2.scatter(X[:, 0], X[:, 1], c=predicted_labels, cmap='viridis', s=10)
ax2.scatter(cluster_centers[0, 0], cluster_centers[0, 1], c='red', s=200, marker='x', label='GMM Center 1')
ax2.scatter(cluster_centers[1, 0], cluster_centers[1, 1], c='blue', s=200, marker='x', label='GMM Center 2')
ax2.set_title("Data Clustered by GMM")
ax2.legend()
ax2.grid(True)

plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
lambda_ = 5  # 平均到着率 (events per time unit)
time_horizon = 10  # シミュレーション時間の範囲 (time units)

# ポアソン過程のイベント発生時間を生成
inter_event_times = np.random.exponential(scale=1/lambda_, size=1000)  # 各イベント間の時間間隔を生成
arrival_times = np.cumsum(inter_event_times)  # 時間間隔を累積和で到着時間に変換

# 指定した時間範囲内のイベントに絞り込む
arrival_times = arrival_times[arrival_times <= time_horizon]

# 可視化
plt.figure(figsize=(10, 5))
plt.step(arrival_times, np.arange(1, len(arrival_times) + 1), where='post', color='blue')
plt.xlabel("Time")
plt.ylabel("Number of Events")
plt.title(f"Poisson Process (λ={lambda_})")
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

# データセットの読み込み
data = load_iris()
X = data.data  # 特徴量データ
y = data.target  # ラベル

# データの標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCAの実行 (2次元)
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

# プロット設定
plt.figure(figsize=(8, 6))
for target in np.unique(y):
    plt.scatter(X_pca[y == target, 0], X_pca[y == target, 1], label=data.target_names[target], s=50, alpha=0.7)

plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.title('PCA of Iris Dataset')
plt.legend()
plt.grid()
plt.show()