システム同定

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.optimize import minimize

# 1. ランダムバイナリー入力信号（PRBS）の生成
def generate_prbs(length, rate):
    """擬似ランダムバイナリシーケンスを生成する関数"""
    prbs = np.random.choice([1, -1], size=length)
    prbs = np.repeat(prbs, rate)  # サンプルを繰り返して出力を増やす（サンプリング周波数を上げる）
    return prbs

# 2. サンプルデータの作成
fs = 100  # サンプリング周波数
duration = 10  # シミュレーション時間 (秒)
time = np.arange(0, duration, 1/fs)  # 時間ベクトル

# ランダムバイナリー入力信号を生成（長さを時間に合わせて設定）
prbs_signal = generate_prbs(length=int(duration/0.1), rate=int(fs*0.1))

# 3. 実際のシステムの伝達関数を設定し、応答を生成
# 例: システムの伝達関数 G(s) = 1 / (s^2 + s + 1)
numerator = [1]
denominator = [1, 1, 1]
system = signal.TransferFunction(numerator, denominator)

# 入力信号に対するシステム応答を計算
_, output_signal, _ = signal.lsim(system, prbs_signal, time)

# 4. システム同定用の伝達関数モデルを定義
def transfer_function_model(params, u, t):
    """伝達関数モデルの応答を計算する関数"""
    # パラメータを伝達関数の係数として使用
    num = [params[0]]  # 分子の係数
    den = [1, params[1], params[2]]  # 分母の係数
    # 伝達関数モデルを生成
    system = signal.TransferFunction(num, den)
    _, y_model, _ = signal.lsim(system, u, t)
    return y_model

# 5. 目的関数（最小化したい誤差関数）を定義
def objective_function(params, u, y, t):
    """目的関数: 実データとモデル応答の二乗誤差を最小化"""
    y_model = transfer_function_model(params, u, t)
    error = np.sum((y - y_model) ** 2)
    return error

# 6. 最適化によるパラメータ推定
initial_guess = [1, 1, 1]  # 初期パラメータ（分子、分母係数）
result = minimize(objective_function, initial_guess, args=(prbs_signal, output_signal, time), method='BFGS')
estimated_params = result.x
print(f"推定された伝達関数の係数: 分子 {estimated_params[0]}, 分母 [1, {estimated_params[1]}, {estimated_params[2]}]")

# 7. 推定された伝達関数を用いて応答をシミュレーション
estimated_system = signal.TransferFunction([estimated_params[0]], [1, estimated_params[1], estimated_params[2]])
_, estimated_output_signal, _ = signal.lsim(estimated_system, prbs_signal, time)

# 8. 結果のプロット
plt.figure(figsize=(12, 6))

# 入力信号のプロット
plt.subplot(3, 1, 1)
plt.plot(time, prbs_signal, label='PRBS Input Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('PRBS Input Signal')
plt.grid(True)
plt.legend()

# 出力信号（実システム応答）のプロット
plt.subplot(3, 1, 2)
plt.plot(time, output_signal, label='Original System Output')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Original System Output')
plt.grid(True)
plt.legend()

# 推定されたモデルの出力をプロット
plt.subplot(3, 1, 3)
plt.plot(time, output_signal, label='Original System Output')
plt.plot(time, estimated_output_signal, '--', label='Estimated Model Output')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('System Identification: Original vs Estimated Model Output')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.optimize import minimize

# 1. ホワイトノイズ入力信号の生成
def generate_white_noise(length, amplitude=1.0):
    """ホワイトノイズを生成する関数"""
    return amplitude * np.random.normal(size=length)

# 2. サンプルデータの作成
fs = 100  # サンプリング周波数 (Hz)
duration = 10  # シミュレーション時間 (秒)
time = np.arange(0, duration, 1/fs)  # 時間ベクトル

# ホワイトノイズ入力信号を生成
white_noise_signal = generate_white_noise(length=len(time), amplitude=1.0)

# 3. 実際のシステムの伝達関数を設定し、応答を生成
# 例: システムの伝達関数 G(s) = 1 / (s^2 + s + 1)
numerator = [1]
denominator = [1, 1, 1]
system = signal.TransferFunction(numerator, denominator)

# 入力信号に対するシステム応答を計算
_, output_signal, _ = signal.lsim(system, white_noise_signal, time)

# 4. システム同定用の伝達関数モデルを定義
def transfer_function_model(params, u, t):
    """伝達関数モデルの応答を計算する関数"""
    # パラメータを伝達関数の係数として使用
    num = [params[0]]  # 分子の係数
    den = [1, params[1], params[2]]  # 分母の係数
    # 伝達関数モデルを生成
    system = signal.TransferFunction(num, den)
    _, y_model, _ = signal.lsim(system, u, t)
    return y_model

# 5. 目的関数（最小化したい誤差関数）を定義
def objective_function(params, u, y, t):
    """目的関数: 実データとモデル応答の二乗誤差を最小化"""
    y_model = transfer_function_model(params, u, t)
    error = np.sum((y - y_model) ** 2)
    return error

# 6. 最適化によるパラメータ推定
initial_guess = [1, 1, 1]  # 初期パラメータ（分子、分母係数）
result = minimize(objective_function, initial_guess, args=(white_noise_signal, output_signal, time), method='BFGS')
estimated_params = result.x
print(f"推定された伝達関数の係数: 分子 {estimated_params[0]}, 分母 [1, {estimated_params[1]}, {estimated_params[2]}]")

# 7. 推定された伝達関数を用いて応答をシミュレーション
estimated_system = signal.TransferFunction([estimated_params[0]], [1, estimated_params[1], estimated_params[2]])
_, estimated_output_signal, _ = signal.lsim(estimated_system, white_noise_signal, time)

# 8. 結果のプロット
plt.figure(figsize=(12, 8))

# 入力信号のプロット
plt.subplot(3, 1, 1)
plt.plot(time, white_noise_signal, label='White Noise Input Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('White Noise Input Signal')
plt.grid(True)
plt.legend()

# 出力信号（実システム応答）のプロット
plt.subplot(3, 1, 2)
plt.plot(time, output_signal, label='Original System Output')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Original System Output')
plt.grid(True)
plt.legend()

# 推定されたモデルの出力をプロット
plt.subplot(3, 1, 3)
plt.plot(time, output_signal, label='Original System Output')
plt.plot(time, estimated_output_signal, '--', label='Estimated Model Output')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('System Identification: Original vs Estimated Model Output')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as scipy_signal  # scipy.signal を明示的に別名でインポート
from scipy.optimize import minimize

# 1. ホワイトノイズ入力信号の生成
def generate_white_noise(length, amplitude=1.0):
    """ホワイトノイズを生成する関数"""
    return amplitude * np.random.normal(size=length)

# 2. データの前処理関数の定義
def preprocess_data(data_signal, fs, cutoff_freq=10, filter_order=4):
    """データの前処理を行う関数（正規化、フィルタリング、平均値の除去）"""
    # ローパスフィルタの設計
    nyquist_freq = fs / 2.0
    normal_cutoff = cutoff_freq / nyquist_freq
    b, a = scipy_signal.butter(filter_order, normal_cutoff, btype='low', analog=False)

    # フィルタリングを行う
    filtered_signal = scipy_signal.filtfilt(b, a, data_signal)

    # 平均値を除去してドリフトを除去
    detrended_signal = filtered_signal - np.mean(filtered_signal)

    # 正規化（0-1範囲）
    normalized_signal = (detrended_signal - np.min(detrended_signal)) / (np.max(detrended_signal) - np.min(detrended_signal))

    return normalized_signal

# 3. サンプルデータの作成
fs = 100  # サンプリング周波数 (Hz)
duration = 10  # シミュレーション時間 (秒)
time = np.arange(0, duration, 1/fs)  # 時間ベクトル

# ホワイトノイズ入力信号を生成
white_noise_signal = generate_white_noise(length=len(time), amplitude=1.0)

# 4. 実際のシステムの伝達関数を設定し、応答を生成
# 例: システムの伝達関数 G(s) = 1 / (s^2 + s + 1)
numerator = [1]
denominator = [1, 1, 1]
system = scipy_signal.TransferFunction(numerator, denominator)

# 入力信号に対するシステム応答を計算
_, output_signal, _ = scipy_signal.lsim(system, white_noise_signal, time)

# 5. データの前処理
preprocessed_input_signal = preprocess_data(white_noise_signal, fs)
preprocessed_output_signal = preprocess_data(output_signal, fs)

# 6. システム同定用の伝達関数モデルを定義
def transfer_function_model(params, u, t):
    """伝達関数モデルの応答を計算する関数"""
    # パラメータを伝達関数の係数として使用
    num = [params[0]]  # 分子の係数
    den = [1, params[1], params[2]]  # 分母の係数
    # 伝達関数モデルを生成
    system = scipy_signal.TransferFunction(num, den)
    _, y_model, _ = scipy_signal.lsim(system, u, t)
    return y_model

# 7. 目的関数（最小化したい誤差関数）を定義
def objective_function(params, u, y, t):
    """目的関数: 実データとモデル応答の二乗誤差を最小化"""
    y_model = transfer_function_model(params, u, t)
    error = np.sum((y - y_model) ** 2)
    return error

# 8. 最適化によるパラメータ推定
initial_guess = [1, 1, 1]  # 初期パラメータ（分子、分母係数）
result = minimize(objective_function, initial_guess, args=(preprocessed_input_signal, preprocessed_output_signal, time), method='BFGS')
estimated_params = result.x
print(f"推定された伝達関数の係数: 分子 {estimated_params[0]}, 分母 [1, {estimated_params[1]}, {estimated_params[2]}]")

# 9. 推定された伝達関数を用いて応答をシミュレーション
estimated_system = scipy_signal.TransferFunction([estimated_params[0]], [1, estimated_params[1], estimated_params[2]])
_, estimated_output_signal, _ = scipy_signal.lsim(estimated_system, preprocessed_input_signal, time)

# 10. 結果のプロット
plt.figure(figsize=(12, 12))

# 元の入力信号のプロット
plt.subplot(4, 1, 1)
plt.plot(time, white_noise_signal, label='Original White Noise Input Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Original White Noise Input Signal')
plt.grid(True)
plt.legend()

# 前処理後の入力信号のプロット
plt.subplot(4, 1, 2)
plt.plot(time, preprocessed_input_signal, label='Preprocessed Input Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Preprocessed Input Signal')
plt.grid(True)
plt.legend()

# 出力信号（実システム応答）のプロット
plt.subplot(4, 1, 3)
plt.plot(time, preprocessed_output_signal, label='Preprocessed Original System Output')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Preprocessed Original System Output')
plt.grid(True)
plt.legend()

# 推定されたモデルの出力をプロット
plt.subplot(4, 1, 4)
plt.plot(time, preprocessed_output_signal, label='Preprocessed Original System Output')
plt.plot(time, estimated_output_signal, '--', label='Estimated Model Output')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('System Identification: Preprocessed Original vs Estimated Model Output')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as scipy_signal  # scipy.signal を明示的に別名でインポート
from scipy.optimize import minimize

# 1. ホワイトノイズ入力信号の生成
def generate_white_noise(length, amplitude=1.0):
    """ホワイトノイズを生成する関数"""
    return amplitude * np.random.normal(size=length)

# 2. データの前処理関数の定義
def preprocess_data(data_signal, fs, cutoff_freq=10, filter_order=4):
    """データの前処理を行う関数（正規化、フィルタリング、平均値の除去）"""
    # ローパスフィルタの設計
    nyquist_freq = fs / 2.0
    normal_cutoff = cutoff_freq / nyquist_freq
    b, a = scipy_signal.butter(filter_order, normal_cutoff, btype='low', analog=False)

    # フィルタリングを行う
    filtered_signal = scipy_signal.filtfilt(b, a, data_signal)

    # 平均値を除去してドリフトを除去
    detrended_signal = filtered_signal - np.mean(filtered_signal)

    # 正規化（0-1範囲）
    normalized_signal = (detrended_signal - np.min(detrended_signal)) / (np.max(detrended_signal) - np.min(detrended_signal))

    return normalized_signal

# 3. サンプルデータの作成
fs = 100  # サンプリング周波数 (Hz)
duration = 10  # シミュレーション時間 (秒)
time = np.arange(0, duration, 1/fs)  # 時間ベクトル

# ホワイトノイズ入力信号を生成
white_noise_signal = generate_white_noise(length=len(time), amplitude=1.0)

# 4. 実際のシステムの伝達関数を設定し、応答を生成
# 例: システムの伝達関数 G(s) = 1 / (s^2 + s + 1)
numerator = [1]
denominator = [1, 1, 1]
system = scipy_signal.TransferFunction(numerator, denominator)

# 入力信号に対するシステム応答を計算
_, output_signal, _ = scipy_signal.lsim(system, white_noise_signal, time)

# 5. データの前処理
preprocessed_input_signal = preprocess_data(white_noise_signal, fs)
preprocessed_output_signal = preprocess_data(output_signal, fs)

# 6. システム同定用の伝達関数モデルを定義
def transfer_function_model(params, u, t):
    """伝達関数モデルの応答を計算する関数"""
    # パラメータを伝達関数の係数として使用
    num = [params[0]]  # 分子の係数
    den = [1, params[1], params[2]]  # 分母の係数
    # 伝達関数モデルを生成
    system = scipy_signal.TransferFunction(num, den)
    _, y_model, _ = scipy_signal.lsim(system, u, t)
    return y_model

# 7. 目的関数（最小化したい誤差関数）を定義
def objective_function(params, u, y, t):
    """目的関数: 実データとモデル応答の二乗誤差を最小化"""
    y_model = transfer_function_model(params, u, t)
    error = np.sum((y - y_model) ** 2)
    return error

# 8. 最適化によるパラメータ推定
initial_guess = [1, 1, 1]  # 初期パラメータ（分子、分母係数）
result = minimize(objective_function, initial_guess, args=(preprocessed_input_signal, preprocessed_output_signal, time), method='BFGS')
estimated_params = result.x
print(f"推定された伝達関数の係数: 分子 {estimated_params[0]}, 分母 [1, {estimated_params[1]}, {estimated_params[2]}]")

# 9. 推定された伝達関数を用いて応答をシミュレーション
estimated_system = scipy_signal.TransferFunction([estimated_params[0]], [1, estimated_params[1], estimated_params[2]])
_, estimated_output_signal, _ = scipy_signal.lsim(estimated_system, preprocessed_input_signal, time)

# 10. 結果のプロット
plt.figure(figsize=(12, 12))

# 元の入力信号のプロット
plt.subplot(4, 1, 1)
plt.plot(time, white_noise_signal, label='Original White Noise Input Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Original White Noise Input Signal')
plt.grid(True)
plt.legend()

# 前処理後の入力信号のプロット
plt.subplot(4, 1, 2)
plt.plot(time, preprocessed_input_signal, label='Preprocessed Input Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Preprocessed Input Signal')
plt.grid(True)
plt.legend()

# 出力信号（実システム応答）のプロット
plt.subplot(4, 1, 3)
plt.plot(time, preprocessed_output_signal, label='Preprocessed Original System Output')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Preprocessed Original System Output')
plt.grid(True)
plt.legend()

# 推定されたモデルの出力をプロット
plt.subplot(4, 1, 4)
plt.plot(time, preprocessed_output_signal, label='Preprocessed Original System Output')
plt.plot(time, estimated_output_signal, '--', label='Estimated Model Output')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('System Identification: Preprocessed Original vs Estimated Model Output')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from statsmodels.tsa.ar_model import AutoReg
from scipy.optimize import minimize

# 1. サンプルデータの生成
np.random.seed(0)  # 再現性のための乱数シード設定
fs = 100  # サンプリング周波数 (Hz)
duration = 10  # シミュレーション時間 (秒)
time = np.arange(0, duration, 1/fs)  # 時間ベクトル

# 入力信号（サイン波+ホワイトノイズ）
input_signal = 0.5 * np.sin(2 * np.pi * 1 * time) + 0.1 * np.random.normal(size=len(time))

# 実システムの伝達関数 G(s) = 1 / (s^2 + s + 1)
numerator = [1]
denominator = [1, 1, 1]
system = signal.TransferFunction(numerator, denominator)

# 入力信号に対するシステムの応答を計算
_, output_signal, _ = signal.lsim(system, input_signal, time)

# 出力信号にノイズを追加
output_signal += 0.05 * np.random.normal(size=len(time))

# 2. データの前処理（正規化）
def normalize(signal):
    return (signal - np.min(signal)) / (np.max(signal) - np.min(signal))

input_signal = normalize(input_signal)
output_signal = normalize(output_signal)

# 3. ARXモデルの定義とパラメータ推定
# ARXモデルの次数（ラグ項数）
na = 2  # 出力側のラグ
nb = 2  # 入力側のラグ
nk = 1  # 入力と出力の遅れ

# 入力と出力のラグ行列を作成
def create_lag_matrix(data, lag):
    """時系列データに基づくラグ行列を作成"""
    N = len(data)
    lag_matrix = np.zeros((N - lag, lag))
    for i in range(lag):
        lag_matrix[:, i] = data[lag - i - 1:N - i - 1]
    return lag_matrix

# ラグ行列の作成
y_lag_matrix = create_lag_matrix(output_signal, na)
u_lag_matrix = create_lag_matrix(input_signal, nb)

# データをそろえる（出力と入力の行数を合わせる）
lagged_output = output_signal[na:]
lagged_input = input_signal[na:]

# 回帰行列を作成（出力ラグ + 入力ラグ）
regression_matrix = np.hstack([y_lag_matrix, u_lag_matrix])

# 回帰係数の推定（最小二乗法）
theta = np.linalg.lstsq(regression_matrix, lagged_output, rcond=None)[0]

# 4. モデルの応答をシミュレーション
# 推定された係数を用いてモデル応答を計算
def simulate_arx_model(theta, u, na, nb, nk):
    """ARXモデルの係数を用いてシステム応答をシミュレーション"""
    N = len(u)
    y_model = np.zeros(N)
    for t in range(max(na, nb + nk), N):
        y_sum = 0
        # 出力側のラグの係数を計算
        for i in range(na):
            y_sum += -theta[i] * y_model[t - i - 1]

        # 入力側のラグの係数を計算
        for j in range(nb):
            y_sum += theta[na + j] * u[t - nk - j - 1]

        # モデルの出力を計算
        y_model[t] = y_sum
    return y_model

# ARXモデルの応答をシミュレーション
estimated_output = simulate_arx_model(theta, input_signal, na, nb, nk)

# 5. 結果のプロット
plt.figure(figsize=(12, 8))

# 入力信号のプロット
plt.subplot(3, 1, 1)
plt.plot(time, input_signal, label='Input Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Input Signal')
plt.grid(True)
plt.legend()

# 出力信号のプロット
plt.subplot(3, 1, 2)
plt.plot(time, output_signal, label='Original Output Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Original Output Signal')
plt.grid(True)
plt.legend()

# モデル応答のプロット
plt.subplot(3, 1, 3)
plt.plot(time, output_signal, label='Original Output Signal')
plt.plot(time, estimated_output, '--', label='Estimated Output (ARX Model)')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('System Identification: Original vs Estimated Output')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import statsmodels.api as sm

# サンプルデータの生成
np.random.seed(0)
n = 100
X = np.random.normal(size=(n, 1))  # 外部入力
e = np.random.normal(size=n)  # ノイズ
y = 2 * X.squeeze() + 0.5 * np.roll(X.squeeze(), 1) + e  # ARXモデル生成

# データの準備
X_lagged = sm.add_constant(np.column_stack([X, np.roll(X, 1)]))  # 定数項と1ラグの追加
X_lagged = X_lagged[1:]  # ラグの影響を考慮して最初のデータを削除
y = y[1:]

# ARXモデルのフィッティング
model_arx = sm.OLS(y, X_lagged).fit()
print(model_arx.summary())

import statsmodels.api as sm
from statsmodels.tsa.arima.model import ARIMA

# サンプルデータの生成
np.random.seed(0)
n = 100
X = np.random.normal(size=(n, 1))
e = np.random.normal(size=n)
y = 2 * X.squeeze() + 0.5 * np.roll(X.squeeze(), 1) + 0.3 * np.roll(e, 1) + e

# ARMAXモデルのフィッティング（ARIMAの拡張）
model_armax = ARIMA(endog=y, exog=X, order=(1, 0, 1)).fit()
print(model_armax.summary())

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX

# サンプルデータの生成
np.random.seed(0)
n = 100
X = np.random.normal(size=(n, 1))
y = 2 * X.squeeze() + 0.5 * np.roll(X.squeeze(), 1) + np.random.normal(size=n)

# 状態空間モデルの定義とフィッティング
model_ssm = SARIMAX(endog=y, exog=X, order=(1, 0, 1), trend='c').fit()
print(model_ssm.summary())

# 結果のプロット
model_ssm.plot_diagnostics(figsize=(10, 8))
plt.show()

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

# サンプルデータの生成
np.random.seed(0)
n = 100
y = np.cumsum(np.random.normal(size=n))  # ARIMAプロセスの生成

# ボックス・ジェンキンス（ARIMA）モデルのフィッティング
model_bj = sm.tsa.ARIMA(y, order=(1, 1, 1)).fit()
print(model_bj.summary())

# プロット
model_bj.plot_diagnostics(figsize=(10, 8))
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# データ生成（例: 入力 x とそれに対応する出力 y）
np.random.seed(0)  # 再現性のためのシード
x = np.linspace(0, 10, 20)  # 入力データ (0~10 の等間隔データ)
true_a, true_b = 2, 1  # 実際のパラメータ
y = true_a * x + true_b + np.random.normal(0, 1, len(x))  # 出力データ（ノイズを含む）

# 最小二乗法によるパラメータ推定
X = np.vstack([x, np.ones(len(x))]).T  # 入力データの行列 [x, 1]
a, b = np.linalg.lstsq(X, y, rcond=None)[0]  # 最小二乗法によるパラメータ推定

# 結果表示
print(f"推定されたパラメータ: a = {a:.2f}, b = {b:.2f}")
plt.scatter(x, y, label='Data (with noise)')  # 元データのプロット
plt.plot(x, a * x + b, color='red', label='Fitted Line')  # 推定されたモデルのプロット
plt.xlabel('Input x')
plt.ylabel('Output y')
plt.legend()
plt.title('Least Squares Method for System Identification')
plt.show()

import numpy as np
from scipy.linalg import hankel
from numpy.linalg import lstsq
import matplotlib.pyplot as plt

# サンプルデータの生成 (ARXモデル)
np.random.seed(0)
n = 100  # データ数
u = np.sin(np.linspace(0, 10, n)) + np.random.normal(0, 0.1, n)  # 外部入力 u(t)
y = np.zeros(n)
a1, a2 = -0.5, 0.25  # 出力の自己回帰係数
b1, b2 = 0.4, 0.2    # 外部入力の係数

# ARXモデルの生成
for t in range(2, n):
    y[t] = a1 * y[t-1] + a2 * y[t-2] + b1 * u[t-1] + b2 * u[t-2] + np.random.normal(0, 0.1)

# 同定用データ行列の構築
phi = np.column_stack([-y[1:n-1], -y[0:n-2], u[1:n-1], u[0:n-2]])
y_train = y[2:]

# 最小二乗法による係数推定
theta = lstsq(phi, y_train, rcond=None)[0]
print(f"推定された係数: {theta}")

# 結果のプロット
y_est = phi @ theta  # 推定出力
plt.plot(y_train, label='True Output')
plt.plot(y_est, '--', label='Estimated Output')
plt.legend()
plt.title('ARX Model Identification')
plt.show()

import statsmodels.api as sm

# サンプルデータ生成
np.random.seed(0)
n = 150
u = np.random.normal(0, 1, n)  # 外部入力
e = np.random.normal(0, 0.1, n)  # ノイズ
y = np.zeros(n)
a1, a2 = -0.5, 0.3
b1, b2 = 0.4, 0.2
c1 = 0.5  # 移動平均項

# ARMAXモデルの生成
for t in range(2, n):
    y[t] = a1 * y[t-1] + a2 * y[t-2] + b1 * u[t-1] + b2 * u[t-2] + c1 * e[t-1] + e[t]

# ARMAXモデルのフィッティング
model = sm.tsa.ARMA(y, order=(2, 1), exog=u)
model_fit = model.fit()
print(model_fit.summary())

# 結果プロット
plt.plot(y, label='True Output')
plt.plot(model_fit.fittedvalues, '--', label='Estimated Output')
plt.legend()
plt.title('ARMAX Model Identification')
plt.show()

from pykalman import KalmanFilter

# サンプルデータの生成
n_timesteps = 50
transition_matrix = [[1, 1], [0, 1]]  # 状態遷移行列
observation_matrix = [[1, 0]]         # 観測行列
initial_state_mean = [0, 0]           # 初期状態

# ノイズの生成
transition_covariance = np.eye(2) * 0.1  # 遷移ノイズの共分散行列
observation_covariance = np.eye(1) * 0.1  # 観測ノイズの共分散行列
initial_state_covariance = np.eye(2) * 0.1

kf = KalmanFilter(transition_matrices=transition_matrix,
                  observation_matrices=observation_matrix,
                  initial_state_mean=initial_state_mean,
                  initial_state_covariance=initial_state_covariance,
                  transition_covariance=transition_covariance,
                  observation_covariance=observation_covariance)

# 状態と観測データの生成
states, observations = kf.sample(n_timesteps)

# カルマンフィルタの適用
filtered_state_means, _ = kf.filter(observations)

# 結果表示
plt.plot(observations, 'o', label='Observations')
plt.plot(filtered_state_means[:, 0], label='Filtered State')
plt.legend()
plt.title('State-Space Model Identification using Kalman Filter')
plt.show()

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

# サンプルデータ生成
n = 100
x = np.linspace(-10, 10, n)
y = np.sin(x) + 0.1 * np.random.normal(size=n)

# ニューラルネットワークモデルの定義
class SimpleNN(nn.Module):
    def __init__(self):
        super(SimpleNN, self).__init__()
        self.fc1 = nn.Linear(1, 10)
        self.fc2 = nn.Linear(10, 1)
    
    def forward(self, x):
        x = torch.relu(self.fc1(x))
        x = self.fc2(x)
        return x

# モデルのインスタンス化と最適化の設定
model = SimpleNN()
criterion = nn.MSELoss()
optimizer = optim.SGD(model.parameters(), lr=0.01)

# データの準備
x_tensor = torch.tensor(x, dtype=torch.float32).unsqueeze(1)
y_tensor = torch.tensor(y, dtype=torch.float32).unsqueeze(1)

# 学習プロセス
for epoch in range(1000):
    model.train()
    optimizer.zero_grad()
    output = model(x_tensor)
    loss = criterion(output, y_tensor)
    loss.backward()
    optimizer.step()

print(f"最終損失: {loss.item():.4f}")

# 結果表示
y_pred = model(x_tensor).detach().numpy()
plt.plot(x, y, 'o', label='True Output')
plt.plot(x, y_pred, '--', label='NN Estimated Output')
plt.legend()
plt.title('Neural Network Based Nonlinear Model Identification')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# シミュレーションデータを生成
np.random.seed(0)
true_params = [2, 3]  # 真のパラメータ
x = np.linspace(0, 10, 100)
y = true_params[0] * x + true_params[1] + np.random.normal(0, 1, x.shape)  # ノイズを追加

# 最小二乗法を用いてパラメータ推定
X = np.vstack([x, np.ones(len(x))]).T  # デザイン行列
estimated_params, _, _, _ = np.linalg.lstsq(X, y, rcond=None)

# 結果の表示
print(f"真のパラメータ: {true_params}")
print(f"推定されたパラメータ: {estimated_params}")

# プロット
plt.scatter(x, y, label='Noisy data')
plt.plot(x, estimated_params[0] * x + estimated_params[1], color='red', label='Fitted line')
plt.legend()
plt.show()

from scipy.optimize import minimize

# 尤度関数（ガウス分布を仮定）
def negative_log_likelihood(params, x, y):
    a, b, sigma = params
    predicted = a * x + b
    likelihood = -np.sum(-0.5 * np.log(2 * np.pi * sigma**2) - ((y - predicted)**2) / (2 * sigma**2))
    return -likelihood

# 初期パラメータの設定
initial_params = [1, 1, 1]

# 尤度関数を最小化（負の対数尤度を最小化）
result = minimize(negative_log_likelihood, initial_params, args=(x, y), method='L-BFGS-B', bounds=[(-10, 10), (-10, 10), (1e-5, 10)])
mle_params = result.x

print(f"推定されたパラメータ（MLE）: {mle_params}")

# プロット
plt.scatter(x, y, label='Noisy data')
plt.plot(x, mle_params[0] * x + mle_params[1], color='green', label='MLE Fitted line')
plt.legend()
plt.show()

from scipy.stats import norm

# 事前分布（パラメータの先験情報）
prior_mean = [0, 0]  # パラメータの平均
prior_var = [1, 1]   # パラメータの分散

# 負の対数事後確率関数
def negative_log_posterior(params, x, y):
    a, b, sigma = params
    # 負の対数尤度
    predicted = a * x + b
    likelihood = -np.sum(-0.5 * np.log(2 * np.pi * sigma**2) - ((y - predicted)**2) / (2 * sigma**2))
    # 事前分布
    prior = -0.5 * np.sum(((params[:2] - prior_mean) ** 2) / prior_var)  # a, bの事前分布
    return -(likelihood + prior)

# 初期パラメータの設定
initial_params_map = [1, 1, 1]

# 負の対数事後確率を最小化
result_map = minimize(negative_log_posterior, initial_params_map, args=(x, y), method='L-BFGS-B', bounds=[(-10, 10), (-10, 10), (1e-5, 10)])
map_params = result_map.x

print(f"推定されたパラメータ（MAP）: {map_params}")

# プロット
plt.scatter(x, y, label='Noisy data')
plt.plot(x, map_params[0] * x + map_params[1], color='blue', label='MAP Fitted line')
plt.legend()
plt.show()

import numpy as np
from scipy.signal import welch, csd

# シミュレーションのためのデータを作成（サンプル数やサンプリング周波数の設定）
fs = 100  # サンプリング周波数
N = 1024  # データポイント数
t = np.linspace(0, N/fs, N, endpoint=False)

# 入力信号としてホワイトノイズを生成
np.random.seed(0)
u = np.random.randn(len(t))

# システムの伝達関数を用いて出力信号を生成
_, y, _ = lti(numerator, denominator).output(U=u, T=t)

# 入力信号と出力信号のパワースペクトル密度を計算
f, Puu = welch(u, fs, nperseg=256)
_, Pyy = welch(y, fs, nperseg=256)

# 入力信号と出力信号の相互パワースペクトル密度を計算
_, Puy = csd(u, y, fs, nperseg=256)

# 周波数応答を求める
G_hat = Puy / Puu  # 推定された周波数応答

# プロット
plt.figure(figsize=(12, 6))

# 周波数応答の振幅（Magnitude）をプロット
plt.subplot(2, 1, 1)
plt.semilogx(f, np.abs(G_hat))
plt.title('Estimated Frequency Response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude')
plt.grid(True)

# 周波数応答の位相（Phase）をプロット
plt.subplot(2, 1, 2)
plt.semilogx(f, np.angle(G_hat, deg=True))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [degrees]')
plt.grid(True)

plt.show()

import numpy as np
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt

# シミュレーション用に非線形システムを定義（例: 2次の非線形システム）
def true_system(input_signal):
    """非線形システム y = 2*u + 0.5*u^2 + 0.3*u^3 の例"""
    return 2 * input_signal + 0.5 * input_signal**2 + 0.3 * input_signal**3

# 入力データの生成（ホワイトノイズを使用）
np.random.seed(0)
N = 200
u = np.random.randn(N)

# システムに入力を加えたときの出力データを生成
y = true_system(u) + 0.1 * np.random.randn(N)  # ノイズを追加

# Volterra系列モデルの設計
def create_volterra_terms(input_signal, degree=3):
    """
    Volterra系列の項を作成する関数
    - input_signal: 入力信号
    - degree: Volterraモデルの次数（1次～指定した次数までの項を作成）
    """
    N = len(input_signal)
    terms = [np.ones(N)]  # 定数項

    # 1次～指定した次数までのVolterra項を作成
    for d in range(1, degree + 1):
        terms.append(input_signal**d)
    
    # (N, d)の行列として返す
    return np.vstack(terms).T

# Volterraモデルの次数を指定
degree = 3

# Volterra系列モデルの項を作成
X = create_volterra_terms(u, degree=degree)

# 線形回帰を用いてVolterraモデルの係数を推定
model = LinearRegression()
model.fit(X, y)

# 推定されたパラメータ
estimated_coefficients = model.coef_
print(f"推定されたVolterra系列モデルの係数: {estimated_coefficients}")

# モデルによる出力を計算
y_pred = model.predict(X)

# 結果のプロット
plt.figure(figsize=(10, 6))
plt.plot(y, label='True Output', color='blue')
plt.plot(y_pred, label='Volterra Model Output', color='red', linestyle='--')
plt.xlabel('Time step')
plt.ylabel('Output')
plt.title('Volterra Series Model Output vs True Output')
plt.legend()
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LinearRegression

# 1. システム同定用のサンプルデータを生成（例: 1次遅れ系のシステム）
np.random.seed(0)
N = 200  # データポイント数
u = np.random.randn(N)  # 入力信号（ホワイトノイズ）
y = np.zeros(N)  # 出力信号を初期化

# 出力信号: y(t) = 0.5 * y(t-1) + 0.3 * u(t-1) + 0.1 * u(t-2) + ノイズ
for t in range(2, N):
    y[t] = 0.5 * y[t-1] + 0.3 * u[t-1] + 0.1 * u[t-2] + 0.1 * np.random.randn()

# データフレームを作成
data = pd.DataFrame({'y': y, 'u': u})

# 2. ARXモデルの遅延データを作成する関数
def create_arx_dataframe(data, output_col, input_col, y_lags, u_lags):
    """
    ARXモデル用の遅延データフレームを作成する関数。
    - data: 入力と出力のデータを含むデータフレーム
    - output_col: 出力データのカラム名（例: 'y'）
    - input_col: 入力データのカラム名（例: 'u'）
    - y_lags: 出力の遅延次数（例: 1, 2, 3）
    - u_lags: 入力の遅延次数（例: 1, 2, 3）
    """
    df = pd.DataFrame()
    # 出力の遅延データを追加
    for i in range(1, y_lags + 1):
        df[f'{output_col}(t-{i})'] = data[output_col].shift(i)
    # 入力の遅延データを追加
    for i in range(1, u_lags + 1):
        df[f'{input_col}(t-{i})'] = data[input_col].shift(i)
    df[output_col] = data[output_col]  # 現在の出力データを追加
    return df.dropna()  # 欠損値（NaN）を削除

# 3. 遅延次数の設定
y_lags = 1  # 出力の遅延次数
u_lags = 2  # 入力の遅延次数

# ARXモデル用のデータフレームを作成
arx_data = create_arx_dataframe(data, 'y', 'u', y_lags, u_lags)

# 説明変数（X）と目的変数（y）に分離
X = arx_data.drop('y', axis=1)
y_arx = arx_data['y']

# 4. 線形回帰モデルを使用してARXモデルのパラメータを推定
model = LinearRegression()
model.fit(X, y_arx)

# 推定された係数を表示
coefficients = model.coef_
intercept = model.intercept_
print(f'推定されたARXモデルの係数: {coefficients}')
print(f'切片（定数項）: {intercept}')

# 5. 推定されたARXモデルによる予測
y_pred = model.predict(X)

# 6. 結果のプロット
plt.figure(figsize=(10, 6))
plt.plot(y_arx.values, label='True Output', color='blue')
plt.plot(y_pred, label='ARX Model Output', color='red', linestyle='--')
plt.xlabel('Time step')
plt.ylabel('Output')
plt.title('ARX Model Output vs True Output')
plt.legend()
plt.grid()
plt.show()