statsmodels 住宅価格予測

statsmodels 住宅価格予測

statsmodelsを使用して現実的な問題を解決する方法を示します。

ここでは、重回帰分析を例に取り上げます。

例えば、住宅価格を予測するためのモデルを構築します。

データセットとして、よく知られているボストン住宅価格データセットを使用します。

このデータセットには、住宅の価格とそれに関連する様々な特徴量が含まれています。

以下の手順で進めます：

  1. データの読み込み
  2. データの前処理
  3. モデルの構築
  4. 結果の解釈

手順

1. データの読み込み

まず、必要なライブラリをインポートし、データセットを読み込みます。

1
2
3
4
5
import statsmodels.api as sm
import pandas as pd

# ボストン住宅価格データセットの読み込み
boston = sm.datasets.get_rdataset("Boston", "MASS").data

2. データの前処理

データを見て、ターゲット変数（住宅価格）と説明変数（特徴量）を分けます。

1
2
3
4
5
6
# データの確認
print(boston.head())

# 説明変数とターゲット変数の分割
X = boston.drop(columns=['medv'])  # 説明変数
y = boston['medv']  # ターゲット変数

3. モデルの構築

statsmodelsを使って重回帰モデルを構築します。

1
2
3
4
5
# 定数項を追加
X = sm.add_constant(X)

# モデルの構築
model = sm.OLS(y, X).fit()

4. 結果の解釈

モデルの結果を表示し、解釈します。

1
2
# 結果の表示
print(model.summary())

完全なコード

以下に完全なコードを示します：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import statsmodels.api as sm
import pandas as pd

# ボストン住宅価格データセットの読み込み
boston = sm.datasets.get_rdataset("Boston", "MASS").data

# データの確認
print(boston.head())

# 説明変数とターゲット変数の分割
X = boston.drop(columns=['medv'])  # 説明変数
y = boston['medv']  # ターゲット変数

# 定数項を追加
X = sm.add_constant(X)

# モデルの構築
model = sm.OLS(y, X).fit()

# 結果の表示
print(model.summary())

[実行結果]

      crim    zn  indus  chas    nox     rm   age     dis  rad  tax  ptratio  \
0  0.00632  18.0   2.31     0  0.538  6.575  65.2  4.0900    1  296     15.3   
1  0.02731   0.0   7.07     0  0.469  6.421  78.9  4.9671    2  242     17.8   
2  0.02729   0.0   7.07     0  0.469  7.185  61.1  4.9671    2  242     17.8   
3  0.03237   0.0   2.18     0  0.458  6.998  45.8  6.0622    3  222     18.7   
4  0.06905   0.0   2.18     0  0.458  7.147  54.2  6.0622    3  222     18.7   

    black  lstat  medv  
0  396.90   4.98  24.0  
1  396.90   9.14  21.6  
2  392.83   4.03  34.7  
3  394.63   2.94  33.4  
4  396.90   5.33  36.2  
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   medv   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Mon, 08 Jul 2024   Prob (F-statistic):          6.72e-135
Time:                        03:02:00   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         36.4595      5.103      7.144      0.000      26.432      46.487
crim          -0.1080      0.033     -3.287      0.001      -0.173      -0.043
zn             0.0464      0.014      3.382      0.001       0.019       0.073
indus          0.0206      0.061      0.334      0.738      -0.100       0.141
chas           2.6867      0.862      3.118      0.002       0.994       4.380
nox          -17.7666      3.820     -4.651      0.000     -25.272     -10.262
rm             3.8099      0.418      9.116      0.000       2.989       4.631
age            0.0007      0.013      0.052      0.958      -0.025       0.027
dis           -1.4756      0.199     -7.398      0.000      -1.867      -1.084
rad            0.3060      0.066      4.613      0.000       0.176       0.436
tax           -0.0123      0.004     -3.280      0.001      -0.020      -0.005
ptratio       -0.9527      0.131     -7.283      0.000      -1.210      -0.696
black          0.0093      0.003      3.467      0.001       0.004       0.015
lstat         -0.5248      0.051    -10.347      0.000      -0.624      -0.425
==============================================================================
Omnibus:                      178.041   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              783.126
Skew:                           1.521   Prob(JB):                    8.84e-171
Kurtosis:                       8.281   Cond. No.                     1.51e+04
==============================================================================

結果の解釈

結果には以下の情報が含まれます：

  • R-squared: モデルがデータのどれだけを説明しているかを示す指標。
    $1$に近いほどモデルの説明力が高い。
  • Coefficients(coef): 各説明変数の回帰係数。
    これにより各特徴量が住宅価格にどのような影響を与えるかがわかります。
  • P-values: 各説明変数の有意性を示す値。
    一般的に、p値が$0.05$未満の変数は統計的に有意とされます。

このようにして、statsmodelsを使用して現実的な問題を解決することができます。

この例では、住宅価格を予測するための重回帰モデルを構築し、その結果を解釈しました。

statsmodels

statsmodels

statsmodelsは、統計モデルの推定や統計検定、データの可視化などを行うためのPythonライブラリです。

以下に、statsmodelsの便利な使い方をいくつか紹介します。

1. 線形回帰モデル

statsmodelsを使用して線形回帰モデルを構築し、データの関係性を解析できます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import statsmodels.api as sm
import pandas as pd

# データの準備
df = pd.DataFrame({
    'X': [1, 2, 3, 4, 5],
    'Y': [1, 2, 1.3, 3.75, 2.25]
})

X = df['X']
Y = df['Y']

# 定数項の追加
X = sm.add_constant(X)

# 回帰モデルの構築とフィッティング
model = sm.OLS(Y, X).fit()

# モデルの概要
print(model.summary())

[実行結果]

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.393
Model:                            OLS   Adj. R-squared:                  0.191
Method:                 Least Squares   F-statistic:                     1.942
Date:                Sun, 07 Jul 2024   Prob (F-statistic):              0.258
Time:                        02:26:29   Log-Likelihood:                -5.6369
No. Observations:                   5   AIC:                             15.27
Df Residuals:                       3   BIC:                             14.49
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.7850      1.012      0.776      0.494      -2.434       4.004
X              0.4250      0.305      1.393      0.258      -0.546       1.396
==============================================================================
Omnibus:                          nan   Durbin-Watson:                   3.369
Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.554
Skew:                           0.641   Prob(JB):                        0.758
Kurtosis:                       1.993   Cond. No.                         8.37
==============================================================================

2. ロジスティック回帰

カテゴリカルデータに対してロジスティック回帰を適用する例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import statsmodels.api as sm
import numpy as np

# データの準備
np.random.seed(42)
df = pd.DataFrame({
    'X1': np.random.rand(100),
    'X2': np.random.rand(100),
    'Y': np.random.randint(0, 2, 100)
})

X = df[['X1', 'X2']]
Y = df['Y']

# 定数項の追加
X = sm.add_constant(X)

# ロジスティック回帰モデルの構築とフィッティング
model = sm.Logit(Y, X).fit()

# モデルの概要
print(model.summary())

[実行結果]

Optimization terminated successfully.
         Current function value: 0.677050
         Iterations 4
                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      Y   No. Observations:                  100
Model:                          Logit   Df Residuals:                       97
Method:                           MLE   Df Model:                            2
Date:                Sun, 07 Jul 2024   Pseudo R-squ.:                 0.01611
Time:                        02:27:12   Log-Likelihood:                -67.705
converged:                       True   LL-Null:                       -68.814
Covariance Type:            nonrobust   LLR p-value:                    0.3299
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.4792      0.524      0.915      0.360      -0.547       1.506
X1             0.4080      0.688      0.593      0.553      -0.939       1.756
X2            -0.9362      0.701     -1.336      0.181      -2.309       0.437

3. 時系列解析

statsmodelsを使用してARIMAモデルなどの時系列解析を行うことができます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import statsmodels.api as sm
import pandas as pd

# 時系列データの準備
data = sm.datasets.co2.load_pandas().data
data = data.fillna(data.bfill())

# ARIMAモデルのフィッティング
model = sm.tsa.ARIMA(data['co2'], order=(1, 1, 1))
result = model.fit()

# モデルの概要
print(result.summary())

# 予測
forecast = result.forecast(steps=10)
print(forecast)

[実行結果]

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                    co2   No. Observations:                 2284
Model:                 ARIMA(1, 1, 1)   Log Likelihood               -1536.855
Date:                Sun, 07 Jul 2024   AIC                           3079.710
Time:                        02:27:49   BIC                           3096.910
Sample:                    03-29-1958   HQIC                          3085.984
                         - 12-29-2001                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.8959      0.022     40.360      0.000       0.852       0.939
ma.L1         -0.7587      0.034    -22.509      0.000      -0.825      -0.693
sigma2         0.2250      0.006     39.341      0.000       0.214       0.236
===================================================================================
Ljung-Box (L1) (Q):                  55.87   Jarque-Bera (JB):                74.67
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.13   Skew:                            -0.14
Prob(H) (two-sided):                  0.10   Kurtosis:                         3.84
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
2002-01-05    371.652800
2002-01-12    371.789687
2002-01-19    371.912318
2002-01-26    372.022178
2002-02-02    372.120597
2002-02-09    372.208766
2002-02-16    372.287753
2002-02-23    372.358514
2002-03-02    372.421906
2002-03-09    372.478695
Freq: W-SAT, Name: predicted_mean, dtype: float64

4. ANOVA（分散分析）

複数グループ間の平均値の差異を検定するためにANOVAを使用します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import statsmodels.api as sm
from statsmodels.formula.api import ols

# データの準備
df = sm.datasets.get_rdataset("iris").data

# 列名の変更
df.columns = [col.replace('.', '_') for col in df.columns]
print(df.columns)  # 変更後の列名の確認

# モデルの構築
model = ols('Sepal_Length ~ C(Species)', data=df).fit()

# ANOVAテーブル
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)

[実行結果]

Index(['Sepal_Length', 'Sepal_Width', 'Petal_Length', 'Petal_Width',
       'Species'],
      dtype='object')
               sum_sq     df           F        PR(>F)
C(Species)  63.212133    2.0  119.264502  1.669669e-31
Residual    38.956200  147.0         NaN           NaN

5. 残差プロット

モデルの適合度を視覚的に評価するために、残差プロットを作成できます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import statsmodels.api as sm
from statsmodels.formula.api import ols
import matplotlib.pyplot as plt

# データの準備
df = sm.datasets.get_rdataset("iris").data
df.columns = [col.replace('.', '_') for col in df.columns]

# モデルの構築（Sepal_LengthとPetal_Lengthの線形回帰）
model = ols('Sepal_Length ~ Petal_Length', data=df).fit()

# 残差のプロット
fig, ax = plt.subplots()
sm.graphics.plot_regress_exog(model, 'Petal_Length', fig=fig)
plt.show()

[実行結果]

6. 様々な統計検定

t検定やカイ二乗検定など、様々な統計検定を行うことができます。

t検定

1
2
3
4
5
6
7
8
9
10
from scipy import stats
import statsmodels.api as sm

# データの準備
df = sm.datasets.get_rdataset("iris").data
df.columns = [col.replace('.', '_') for col in df.columns]

# 2つの独立した標本のt検定
t_stat, p_value = stats.ttest_ind(df['Sepal_Length'], df['Petal_Length'])
print(f'T-statistic: {t_stat}, P-value: {p_value}')

[実行結果]

T-statistic: 13.098353108960858, P-value: 2.857104069581941e-31

カイ二乗検定

1
2
3
4
5
6
7
8
9
10
11
12
import pandas as pd
from scipy.stats import chi2_contingency
import statsmodels.api as sm

# データの準備
df = sm.datasets.get_rdataset("iris").data
df.columns = [col.replace('.', '_') for col in df.columns]

# カイ二乗検定
table = pd.crosstab(df['Species'], df['Sepal_Length'] > df['Sepal_Length'].median())
chi2, p, dof, ex = chi2_contingency(table)
print(f'Chi2: {chi2}, P-value: {p}')

[実行結果]

Chi2: 78.64285714285712, P-value: 8.37376079899524e-18

これらの例は、statsmodelsの基本的な使用方法を示していますが、このライブラリは他にも多くの機能を提供しており、統計解析やモデリングの幅広いニーズに応えることができます。

仮想通貨の価格分析

仮想通貨の価格分析

仮想通貨の分析には、価格データの取得、トレンドの分析、テクニカル指標の計算などが含まれます。

以下に、Pythonを使って簡単な仮想通貨の価格分析を行う方法を示します。

ここでは、Bitcoin（BTC）の価格データを取得し、移動平均を計算してプロットする例を紹介します。

必要なライブラリのインストール

まず、以下のライブラリをインストールしてください。

1
pip install pandas yfinance matplotlib

データの取得と分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import yfinance as yf
import pandas as pd
import matplotlib.pyplot as plt

# データの取得
btc = yf.Ticker("BTC-USD")
data = btc.history(period="1y")  # 過去1年間のデータを取得

# 移動平均の計算
data['MA50'] = data['Close'].rolling(window=50).mean()  # 50日移動平均
data['MA200'] = data['Close'].rolling(window=200).mean()  # 200日移動平均

# データのプロット
plt.figure(figsize=(14, 7))
plt.plot(data['Close'], label='BTC-USD Close Price')
plt.plot(data['MA50'], label='50-Day Moving Average')
plt.plot(data['MA200'], label='200-Day Moving Average')
plt.title('Bitcoin Price and Moving Averages')
plt.xlabel('Date')
plt.ylabel('Price (USD)')
plt.legend()
plt.show()

解説

  1. データの取得: yfinanceライブラリを使って、Bitcoinの過去1年間の価格データを取得します。
  2. 移動平均の計算: pandasを使って、50日移動平均と200日移動平均を計算します。
  3. データのプロット: matplotlibを使って、価格データと移動平均をプロットします。

この例では、Bitcoinの過去1年間の終値と、その50日および200日の移動平均を表示しています。
移動平均を使うことで、価格の長期的なトレンドを視覚的に把握することができます。

他の分析方法としては、相対力指数（RSI）、移動平均収束拡散指標（MACD）などのテクニカル指標の計算や、ボラティリティの分析なども可能です。

さらに高度な分析には、機械学習モデルを使った価格予測や、統計的手法を用いた異常検知などがあります。

[実行結果]

ソースコード解説

このPythonスクリプトは、Bitcoinの価格データを取得し、50日および200日移動平均を計算してプロットするものです。

以下、コードの各部分を詳しく説明します。

1. ライブラリのインポート

1
2
3
import yfinance as yf
import pandas as pd
import matplotlib.pyplot as plt

この章では、必要なライブラリをインポートしています。

  • yfinance: Yahoo Financeから金融データを取得するためのライブラリ。
  • pandas: データの操作や分析を行うためのライブラリ。
  • matplotlib.pyplot: データの可視化を行うためのライブラリ。

2. データの取得

1
2
btc = yf.Ticker("BTC-USD")
data = btc.history(period="1y")  # 過去1年間のデータを取得

この章では、Bitcoinの価格データを取得しています。

  • yf.Ticker("BTC-USD"): yfinanceライブラリを使用して、Bitcoin（ティッカーシンボル: BTC-USD）のデータを取得するためのオブジェクトを作成します。
  • btc.history(period="1y"): 過去1年間のBitcoinの履歴データを取得します。このデータには、日次の価格（始値、高値、安値、終値）、取引量などが含まれます。

3. 移動平均の計算

1
2
data['MA50'] = data['Close'].rolling(window=50).mean()  # 50日移動平均
data['MA200'] = data['Close'].rolling(window=200).mean()  # 200日移動平均

この章では、50日移動平均と200日移動平均を計算しています。

  • data['Close']: データフレームから終値の列を選択します。
  • rolling(window=50): 終値の列に対して、50日間の移動ウィンドウを作成します。
  • mean(): 移動ウィンドウ内の値の平均を計算します。この操作により、50日間の移動平均が計算されます。
  • data['MA50']: 50日移動平均を新しい列としてデータフレームに追加します。
  • 同様に、window=200として200日移動平均も計算し、data['MA200']列に追加します。

4. データのプロット

1
2
3
4
5
6
7
8
9
plt.figure(figsize=(14, 7))
plt.plot(data['Close'], label='BTC-USD Close Price')
plt.plot(data['MA50'], label='50-Day Moving Average')
plt.plot(data['MA200'], label='200-Day Moving Average')
plt.title('Bitcoin Price and Moving Averages')
plt.xlabel('Date')
plt.ylabel('Price (USD)')
plt.legend()
plt.show()

この章では、取得したデータと計算した移動平均をプロットしています。

  • plt.figure(figsize=(14, 7)): グラフのサイズを設定します。幅14インチ、高さ7インチの図を作成します。
  • plt.plot(data['Close'], label='BTC-USD Close Price'): Bitcoinの終値をプロットし、ラベルを「BTC-USD Close Price」に設定します。
  • plt.plot(data['MA50'], label='50-Day Moving Average'): 50日移動平均をプロットし、ラベルを「50-Day Moving Average」に設定します。
  • plt.plot(data['MA200'], label='200-Day Moving Average'): 200日移動平均をプロットし、ラベルを「200-Day Moving Average」に設定します。
  • plt.title('Bitcoin Price and Moving Averages'): グラフのタイトルを設定します。
  • plt.xlabel('Date'): x軸のラベルを「Date」に設定します。
  • plt.ylabel('Price (USD)'): y軸のラベルを「Price (USD)」に設定します。
  • plt.legend(): 凡例を表示します。
  • plt.show(): グラフを表示します。

これにより、Bitcoinの価格推移とその50日および200日移動平均が可視化されます。

グラフ解説

[実行結果]

このグラフは、Bitcoin（BTC-USD）の価格推移と移動平均を表示しています。

以下は、このグラフの詳細な説明です。

グラフの概要

  • 期間: 過去1年間（2023年7月から2024年7月）
  • データ: Bitcoinの終値（青線）
  • 移動平均:
    • 50日移動平均（オレンジ線）
    • 200日移動平均（緑線）

詳細な説明

  1. 価格の変動（青線）:

    • 2023年7月から2023年11月まで: Bitcoinの価格は$30,000 $USD以下で、比較的低い水準で推移しています。
    • 2023年11月から2024年3月にかけて: 大きな上昇トレンドが見られ、価格は$70,000 $USDに達しています。この期間は急激な価格上昇を示しています。
    • 2024年3月以降: 価格は乱高下を繰り返しながらも、全体的に下降傾向にあります。

  2. 50日移動平均（オレンジ線）:

    • 短期的な価格のトレンドを示しています。価格が50日移動平均を上回っているときは上昇トレンド、下回っているときは下降トレンドと解釈されます。
    • 2023年末から2024年3月にかけて上昇し、その後は価格の変動に応じて上下していますが、全体的には価格の動きに遅れて追従しています。

  3. 200日移動平均（緑線）:

    • 長期的な価格のトレンドを示しています。50日移動平均よりも滑らかな曲線で、価格の大きなトレンドを捉えています。
    • 2023年11月から2024年7月にかけて一貫して上昇しており、価格の全体的な上昇トレンドを示しています。

結論

  • 価格動向: Bitcoinの価格は2023年7月から2023年11月まで低迷していましたが、その後急激に上昇し、2024年3月にはピークに達しました。その後は価格が乱高下し、下降傾向にあります。
  • 移動平均の交差: 50日移動平均が200日移動平均を上回るとゴールデンクロスと呼ばれ、強気市場のサインとされます。
    一方、50日移動平均が200日移動平均を下回るとデッドクロスと呼ばれ、弱気市場のサインとされます。
    グラフからは、直近でデッドクロスが発生していることが確認でき、今後の価格動向に注意が必要です。

この分析は、投資の判断材料として役立つかもしれませんが、他の要因や指標も考慮することをお勧めします。

CVXPY ポートフォリオ最適化

CVXPY ポートフォリオ最適化

次のようなポートフォリオ最適化問題を考えてみましょう。

ポートフォリオ最適化問題

あなたは、以下の3つの資産に投資するポートフォリオを構築しようとしています：

  1. Stock_A
  2. Stock_B
  3. Bond_C

それぞれの資産の予想リターンとリスク（分散）は次の通りです：

  • Stock_A: リターン = $12%$、リスク（分散） = $0.04$
  • Stock_B: リターン = $10%$、リスク（分散） = $0.03$
  • Bond_C: リターン = $6%$、リスク（分散） = $0.01$

また、各資産間の共分散は以下のようになっています：

Stock_A Stock_B Bond_C
Stock_A 0.04 0.01 0.00
Stock_B 0.01 0.03 0.01
Bond_C 0.00 0.01 0.01

あなたの目標は、全体のリスクを最小限に抑えつつ、ポートフォリオの期待リターンが少なくとも9%となるように投資割合を決定することです。

制約条件

  1. 各資産への投資割合の合計は$1$（100%）でなければならない。
  2. 各資産への投資割合は$0$以上でなければならない（空売りはしない）。
  3. ポートフォリオの期待リターンは少なくとも9%でなければならない。

目的

上記の制約条件を満たしながら、ポートフォリオのリスク（分散）を最小化する投資割合を求めなさい。

ソースコード

2次計画法（Quadratic Programming, QP）を用いて上記の問題を解きます。

以下は、cvxpyを用いたポートフォリオ最適化のコード例です：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
import cvxpy as cp
import numpy as np

# 投資可能な資産のリスト
assets = ['Stock_A', 'Stock_B', 'Bond_C']

# 予想リターンとリスク（分散）
returns = np.array([0.12, 0.10, 0.06])
risks = np.array([0.04, 0.03, 0.01])
covariance = np.array([
    [0.04, 0.01, 0.00],
    [0.01, 0.03, 0.01],
    [0.00, 0.01, 0.01]
])

# 変数の定義
weights = cp.Variable(len(assets))

# リスク（分散）の計算
risk = cp.quad_form(weights, covariance)

# 制約条件の定義
constraints = [
    cp.sum(weights) == 1,
    weights >= 0,
    cp.sum(weights * returns) >= 0.09
]

# 目的関数の定義（リスクの最小化）
prob = cp.Problem(cp.Minimize(risk), constraints)

# 問題の解決
prob.solve()

# 結果の表示
print("Status:", prob.status)
print("Optimal Weights:", weights.value)
print("Total Risk (Variance):", risk.value)

このコードは、cvxpyを用いて2次計画法を解き、ポートフォリオの最適な投資割合とリスクを求めます。

ソースコード解説

このソースコードは、ポートフォリオ最適化問題を解くために、Pythonの数理最適化ライブラリであるcvxpyを使用しています。

ソースコードを詳しく説明します。

1. 必要なライブラリのインポート

1
2
import cvxpy as cp
import numpy as np

最初に、数理最適化を行うためのcvxpyライブラリと、数値計算を行うためのnumpyライブラリをインポートします。

2. 投資可能な資産のリスト

1
assets = ['Stock_A', 'Stock_B', 'Bond_C']

ここでは、投資対象となる資産のリストを定義しています。

今回は、3つの資産（Stock_A, Stock_B, Bond_C）に投資します。

3. 予想リターンとリスク（分散）の定義

1
2
3
4
5
6
7
returns = np.array([0.12, 0.10, 0.06])
risks = np.array([0.04, 0.03, 0.01])
covariance = np.array([
    [0.04, 0.01, 0.00],
    [0.01, 0.03, 0.01],
    [0.00, 0.01, 0.01]
])
  • returnsは各資産の予想リターンを表します。
  • risksは各資産のリスク（分散）を表します。
  • covarianceは各資産間の共分散行列を表します。
    共分散行列は、資産間のリスクの相関関係を示します。

4. 変数の定義

1
weights = cp.Variable(len(assets))

ここでは、各資産に対する投資割合を表す変数weightsを定義しています。
weightsは、3つの要素を持つベクトルです（各資産に対する投資割合を示すため）。

5. リスク（分散）の計算

1
risk = cp.quad_form(weights, covariance)

ここでは、ポートフォリオ全体のリスクを計算しています。
cp.quad_formは、二次形式を計算するための関数で、ポートフォリオのリスク（分散）を求めるのに使用されます。

6. 制約条件の定義

1
2
3
4
5
constraints = [
    cp.sum(weights) == 1,
    weights >= 0,
    cp.sum(weights * returns) >= 0.09
]

ここでは、最適化問題の制約条件を定義しています。

  • cp.sum(weights) == 1：各資産の投資割合の合計は$1$（100%）でなければなりません。
  • weights >= 0：各資産の投資割合は$0$以上でなければなりません（空売りはしない）。
  • cp.sum(weights * returns) >= 0.09：ポートフォリオの期待リターンは少なくとも9%でなければなりません。

7. 目的関数の定義（リスクの最小化）

1
prob = cp.Problem(cp.Minimize(risk), constraints)

ここでは、最適化問題を定義しています。
目的関数はポートフォリオのリスク（分散）を最小化することです。
また、制約条件もここで設定します。

8. 問題の解決

1
prob.solve()

この行では、定義された最適化問題を解決します。
prob.solve()は、最適化問題を解くための関数です。

9. 結果の表示

1
2
3
print("Status:", prob.status)
print("Optimal Weights:", weights.value)
print("Total Risk (Variance):", risk.value)

最後に、最適化結果を表示します。

  • prob.status：最適化問題の解決ステータス（成功かどうか）を表示します。
  • weights.value：各資産に対する最適な投資割合を表示します。
  • risk.value：最適な投資割合で得られるポートフォリオの総リスク（分散）を表示します。

全体の流れ

このコードは、ポートフォリオ最適化問題を定義し、投資割合を決定することで、リスクを最小化しながら目標リターンを達成するためのものです。

数理最適化ライブラリcvxpyを使用して、制約条件を満たしつつ、ポートフォリオのリスクを最小限に抑える最適な投資割合を計算します。

結果解説

[実行結果]

Status: optimal
Optimal Weights: [0.38461538 0.17307692 0.44230769]
Total Risk (Variance): 0.011634615384615387

上記の実行結果は、ポートフォリオ最適化問題を解いた際に得られた最適な解答を示しています。

以下にそれぞれの項目について説明します。

1. Status: optimal

この「optimal」というステータスは、最適化問題が解決されたことを示しています。

具体的には、全ての制約条件を満たしながら、リスク（分散）を最小限に抑える投資割合が見つかったことを意味します。

2. Optimal Weights: [0.38461538 0.17307692 0.44230769]

これは、各資産に対する最適な投資割合を示しています。
具体的には、以下の通りです：

  • Stock_A: 38.46%
  • Stock_B: 17.31%
  • Bond_C: 44.23%

これらの割合に基づいて投資することで、リスクを最小限に抑えつつ、目標リターン（少なくとも9%）を達成することができます。

3. Total Risk (Variance): 0.011634615384615387

これは、最適な投資割合で得られるポートフォリオの総リスク（分散）を示しています。
具体的には、$0.0116$という値になっています。
この値が小さいほど、ポートフォリオ全体のリスクが低いことを意味します。

結論

この最適化結果から、以下の点が確認できます：

  • 最適な投資割合は、Stock_Aに約$38.46%$、Stock_Bに約$17.31%$、Bond_Cに約$44.23%$です。
  • この投資割合により、ポートフォリオの期待リターンが$9%$以上となり、かつリスク（分散）が最小化されます。

したがって、この結果を基に投資を行うことで、目標リターンを達成しつつ、リスクを最小限に抑えたポートフォリオを構築することができます。

PuLP 時間割作成問題

PuLP 時間割作成問題

学校や大学の時間割を作成する際には、特定の科目が特定の時間に配置されるようにする必要があります。

ここでは、時間割作成問題をPuLPで解決する方法を示します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import pulp as pl

# 教師、科目、時間スロットのリスト
teachers = ['T1', 'T2']
subjects = ['Math', 'Science', 'English']
time_slots = ['Monday 9am', 'Monday 11am', 'Wednesday 9am']

# 教師が教えられる科目
teacher_subjects = {
    'T1': ['Math', 'Science'],
    'T2': ['Science', 'English']
}

# 問題の定義
prob = pl.LpProblem("Timetable_Problem", pl.LpMinimize)

# 変数の定義
vars = pl.LpVariable.dicts("Schedule", (teachers, subjects, time_slots), 0, 1, pl.LpBinary)

# 制約条件の定義
# 1. 各科目は1つの時間スロットに配置される
for s in subjects:
    prob += pl.lpSum([vars[t][s][ts] for t in teachers for ts in time_slots]) == 1

# 2. 各教師は1つの時間スロットで1つの科目のみ教える
for t in teachers:
    for ts in time_slots:
        prob += pl.lpSum([vars[t][s][ts] for s in subjects]) <= 1

# 3. 教師は教えられる科目のみ教える
for t in teachers:
    for s in subjects:
        if s not in teacher_subjects[t]:
            for ts in time_slots:
                prob += vars[t][s][ts] == 0

# 問題の解決
prob.solve()

# 結果の表示
print("Status:", pl.LpStatus[prob.status])
for t in teachers:
    for s in subjects:
        for ts in time_slots:
            if vars[t][s][ts].varValue == 1:
                print(f"{t} will teach {s} at {ts}")

ソースコード解説

このソースコードは、教師、科目、時間スロットに関する時間割問題を解決するためにPuLPを使用しています。

以下に、コードの各章を詳しく説明します。

1. インポートとデータの定義

1
2
3
4
5
6
7
8
9
10
11
12
import pulp as pl

# 教師、科目、時間スロットのリスト
teachers = ['T1', 'T2']
subjects = ['Math', 'Science', 'English']
time_slots = ['Monday 9am', 'Monday 11am', 'Wednesday 9am']

# 教師が教えられる科目
teacher_subjects = {
    'T1': ['Math', 'Science'],
    'T2': ['Science', 'English']
}

説明:

  • pulp ライブラリを pl としてインポートします。
  • 教師、科目、時間スロットのリストを定義します。
  • teacher_subjects 辞書を用いて、各教師が教えることができる科目を定義します。

2. 問題の定義

1
2
# 問題の定義
prob = pl.LpProblem("Timetable_Problem", pl.LpMinimize)

説明:

  • LpProblem クラスを使って、新しい最適化問題を定義します。
    この問題の目的は、時間割の問題を解決することです。

目的関数は LpMinimize で、最小化することを示していますが、実際の目的関数はこのコードでは特に定義されていません。

3. 変数の定義

1
2
# 変数の定義
vars = pl.LpVariable.dicts("Schedule", (teachers, subjects, time_slots), 0, 1, pl.LpBinary)

説明:

  • LpVariable.dicts を使って、3次元のバイナリ変数を定義します。

vars[t][s][ts] は、教師 t が科目 s を時間スロット ts に教えるかどうかを示すバイナリ変数です。

$0 $または$ 1 $の値を取ります。

4. 制約条件の定義

4.1 各科目は1つの時間スロットに配置される

1
2
3
# 1. 各科目は1つの時間スロットに配置される
for s in subjects:
    prob += pl.lpSum([vars[t][s][ts] for t in teachers for ts in time_slots]) == 1

説明:

  • 各科目 s は、全ての教師と時間スロットに対して、ちょうど1回だけ教えられることを保証する制約です。
    lpSum 関数を使って、全ての教師と時間スロットにわたる変数の合計が$1$になるようにします。

4.2 各教師は$1$つの時間スロットで$1$つの科目のみ教える

1
2
3
4
# 2. 各教師は1つの時間スロットで1つの科目のみ教える
for t in teachers:
    for ts in time_slots:
        prob += pl.lpSum([vars[t][s][ts] for s in subjects]) <= 1

説明:

  • 各教師 t は、特定の時間スロット ts で、最大1つの科目 s を教えることを保証する制約です。

lpSum 関数を使って、各時間スロットにおける変数の合計が$1$以下になるようにします。

4.3 教師は教えられる科目のみ教える

1
2
3
4
5
6
# 3. 教師は教えられる科目のみ教える
for t in teachers:
    for s in subjects:
        if s not in teacher_subjects[t]:
            for ts in time_slots:
                prob += vars[t][s][ts] == 0

説明:

  • 各教師 t は、自分が教えることができない科目 s を教えないことを保証する制約です。teacher_subjects 辞書を参照し、教師が教えられない科目に対する変数を$0$に設定します。

5. 問題の解決

1
2
# 問題の解決
prob.solve()

説明:

  • 定義された問題を解決します。

PuLP の solve メソッドを使って、最適な解を見つけます。

6. 結果の表示

1
2
3
4
5
6
7
# 結果の表示
print("Status:", pl.LpStatus[prob.status])
for t in teachers:
    for s in subjects:
        for ts in time_slots:
            if vars[t][s][ts].varValue == 1:
                print(f"{t} will teach {s} at {ts}")

説明:

  • 解のステータスを表示します。pl.LpStatus[prob.status] は解の状態を示します（例えば、最適解が見つかったかどうか）。
  • 各教師、科目、時間スロットの組み合わせに対して、変数の値が$1$の場合、その教師がその時間にその科目を教えることを示すメッセージを表示します。

まとめ

このコードは、PuLP を使用して時間割作成問題を解決するための例です。

教師、科目、時間スロットの組み合わせに対してバイナリ変数を定義し、制約条件を設定して、問題を解決します。

最終的に、最適な時間割を表示します。

結果解説

[実行結果]

Status: Optimal
T1 will teach Math at Monday 11am
T2 will teach Science at Wednesday 9am
T2 will teach English at Monday 11am

この結果は、PuLPを使って時間割作成問題を解決した際の出力です。

出力の内容を以下に説明します。

結果の解釈

  • Status: Optimal

    • このステータスは、最適解が見つかったことを示しています。
      つまり、与えられた制約条件の下で最良の時間割が作成されたことを意味します。

  • T1 will teach Math at Monday 11am

    • 教師T1が、月曜日の午前11時に数学を教えることが決定されたことを示しています。

  • T2 will teach Science at Wednesday 9am

    • 教師T2が、水曜日の午前9時に科学を教えることが決定されたことを示しています。

  • T2 will teach English at Monday 11am

    • 教師T2が、月曜日の午前11時に英語を教えることが決定されたことを示しています。

解の詳細

この結果は、以下のような制約条件と目標に基づいています。

  1. 各科目は1つの時間スロットに配置される

    • どの科目も$1$つの時間スロットに割り当てられています。

  2. 各教師は1つの時間スロットで1つの科目のみ教える

    • 各教師は$1$つの時間スロットで$1$つの科目しか教えていません。
      例えば、T1は月曜日11amにMathを教え、他の時間帯には何も教えていません。

  3. 教師は教えられる科目のみ教える

    • 各教師は指定された科目のみを教えています。
      T2はScienceとEnglishを教えており、T1はMathを教えています。

このスケジュールは、与えられた制約をすべて満たしつつ、可能な最適な配置となるように調整されています。

PuLPがこの問題を効率的に解決し、最適な時間割を提供していることが分かります。

NumPy

NumPy

NumPyは、Pythonで科学計算を行うための強力なライブラリです。

以下にNumPyの便利な使い方をいくつか紹介します。

1. 配列の作成

NumPyの基本は配列（ndarray）です。

以下の方法で配列を作成できます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np

# 配列の作成
a = np.array([1, 2, 3, 4, 5])
print(a)

# ゼロ配列
b = np.zeros((3, 3))
print(b)

# 1配列
c = np.ones((2, 2))
print(c)

# 連続した数値の配列
d = np.arange(10)
print(d)

# 等間隔の数値の配列
e = np.linspace(0, 1, 5)
print(e)

[実行結果]

[1 2 3 4 5]
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
[[1. 1.]
 [1. 1.]]
[0 1 2 3 4 5 6 7 8 9]
[0.   0.25 0.5  0.75 1.  ]

2. 配列の形状変更

配列の形状を変更する方法です。

1
2
3
4
5
6
7
8
# 1次元配列を2次元に変更
a = np.array([1, 2, 3, 4, 5, 6])
b = a.reshape((2, 3))
print(b)

# 転置
c = b.T
print(c)

[実行結果]

[[1 2 3]
 [4 5 6]]
[[1 4]
 [2 5]
 [3 6]]

3. 配列の演算

配列同士の演算や、ブロードキャストを用いた演算が簡単にできます。

1
2
3
4
5
6
7
8
9
10
11
12
13
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

# 要素ごとの演算
c = a + b
print(c)

d = a * b
print(d)

# ブロードキャスト
e = a + 10
print(e)

[実行結果]

[5 7 9]
[ 4 10 18]
[11 12 13]

4. 統計関数

NumPyには多くの統計関数があります。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
a = np.array([1, 2, 3, 4, 5])

# 平均
mean = np.mean(a)
print(mean)

# 中央値
median = np.median(a)
print(median)

# 標準偏差
std = np.std(a)
print(std)

# 和
sum_ = np.sum(a)
print(sum_)

[実行結果]

3.0
3.0
1.4142135623730951
15

5. インデックスとスライシング

NumPy配列のインデックスとスライシングは非常に柔軟です。

1
2
3
4
5
6
7
8
9
10
a = np.array([1, 2, 3, 4, 5])

# インデックス
print(a[0])
print(a[-1])

# スライシング
print(a[1:3])
print(a[:2])
print(a[::2])

[実行結果]

1
5
[2 3]
[1 2]
[1 3 5]

6. 条件を使ったフィルタリング

条件を使って配列の要素をフィルタリングできます。

1
2
3
4
5
6
7
8
9
10
a = np.array([1, 2, 3, 4, 5])

# 条件を使ったフィルタリング
b = a[a > 2]
print(b)

# ブールインデックス
c = np.array([True, False, True, False, True])
d = a[c]
print(d)

[実行結果]

[3 4 5]
[1 3 5]

7. 線形代数

NumPyには線形代数のための関数も豊富に用意されています。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 行列の作成
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# 行列積
C = np.dot(A, B)
print(C)

# 逆行列
D = np.linalg.inv(A)
print(D)

# 固有値と固有ベクトル
eigvals, eigvecs = np.linalg.eig(A)
print(eigvals)
print(eigvecs)

[実行結果]

[[19 22]
 [43 50]]
[[-2.   1. ]
 [ 1.5 -0.5]]
[-0.37228132  5.37228132]
[[-0.82456484 -0.41597356]
 [ 0.56576746 -0.90937671]]

これらはNumPyの基本的な使い方の一部に過ぎません。

NumPyは非常に強力で多機能なライブラリなので、公式ドキュメントやチュートリアルを参考にしながら、さらに深く学んでいくことをお勧めします。

Pyomo

Pyomo

PyomoはPythonで線形計画法や整数計画法などの数理最適化問題をモデル化、解決するための強力なツールです。

以下はPyomoの便利な使い方をいくつか紹介します。

1. 簡単な線形計画問題のモデル化

まず、Pyomoを使って簡単な線形計画問題をモデル化してみましょう。

問題

次の線形計画問題を解いてみます。
最大化：
$$
\text{maximize} \quad 2x + 3y
$$
制約：
$$
x + 2y \leq 20
$$
$$
3x + y \leq 15
$$
$$
x \geq 0, \quad y \geq 0
$$

解法

まず、Pyomoをインストールします。

1
2
!pip install pyomo
!apt-get install -y -qq glpk-utils

次に、Pyomoのコードを書きます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from pyomo.environ import *

# モデルの作成
model = ConcreteModel()

# 変数の定義
model.x = Var(within=NonNegativeReals)
model.y = Var(within=NonNegativeReals)

# 目的関数の定義
model.objective = Objective(expr=2*model.x + 3*model.y, sense=maximize)

# 制約条件の定義
model.constraint1 = Constraint(expr=model.x + 2*model.y <= 20)
model.constraint2 = Constraint(expr=3*model.x + model.y <= 15)

# 求解器の呼び出し
solver = SolverFactory('glpk')
solver.solve(model)

# 結果の表示
print("Optimal value for x:", model.x())
print("Optimal value for y:", model.y())
print("Optimal objective value:", model.objective())

[実行結果]

Optimal value for x: 2.0
Optimal value for y: 9.0
Optimal objective value: 31.0

2. 制約付き最適化問題のモデル化

Pyomoは制約の数が多くても簡単にモデル化できます。

例えば、多くの変数や複雑な制約がある場合でも、Pythonのリストや辞書を使って効率的に定義できます。

問題

多変数の制約付き最適化問題を解く例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
from pyomo.environ import *

# モデルの作成
model = ConcreteModel()

# データの定義
N = 5  # 変数の数
coefficients = [1, 2, 3, 4, 5]
constraints = [
    [2, 1, 2, 1, 1, 10],
    [1, 3, 2, 2, 1, 15]
]

# 変数の定義
model.x = Var(range(N), within=NonNegativeReals)

# 目的関数の定義
model.objective = Objective(expr=sum(coefficients[i] * model.x[i] for i in range(N)), sense=maximize)

# 制約条件の定義
model.constraints = ConstraintList()
for con in constraints:
    model.constraints.add(sum(con[i] * model.x[i] for i in range(N)) <= con[-1])

# 求解器の呼び出し
solver = SolverFactory('glpk')
solver.solve(model)

# 結果の表示
for i in range(N):
    print(f"Optimal value for x[{i}]:", model.x[i]())
print("Optimal objective value:", model.objective())

[実行結果]

Optimal value for x[0]: 0.0
Optimal value for x[1]: 0.0
Optimal value for x[2]: 0.0
Optimal value for x[3]: 0.0
Optimal value for x[4]: 10.0
Optimal objective value: 50.0

3. 整数計画問題のモデル化

Pyomoは整数計画問題も扱えます。

例えば、整数変数を使ったモデルの作成です。

問題

整数計画問題の例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from pyomo.environ import *

# モデルの作成
model = ConcreteModel()

# 変数の定義
model.x = Var(within=NonNegativeIntegers)
model.y = Var(within=NonNegativeIntegers)

# 目的関数の定義
model.objective = Objective(expr=2*model.x + 3*model.y, sense=maximize)

# 制約条件の定義
model.constraint1 = Constraint(expr=model.x + 2*model.y <= 20)
model.constraint2 = Constraint(expr=3*model.x + model.y <= 15)

# 求解器の呼び出し
solver = SolverFactory('glpk')
solver.solve(model)

# 結果の表示
print("Optimal value for x:", model.x())
print("Optimal value for y:", model.y())
print("Optimal objective value:", model.objective())

[実行結果]

Optimal value for x: 2.0
Optimal value for y: 9.0
Optimal objective value: 31.0

4. 複雑な制約条件のモデル化

Pyomoを使えば、複雑な非線形制約条件もモデル化できます。

問題

非線形制約条件を含む最適化問題の例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
!apt-get update
!apt-get install -y coinor-cbc
!pip install -q pyomo

from pyomo.environ import *

# モデルの作成
model = ConcreteModel()

# データの定義
N = 5  # 変数の数
coefficients = [1, 2, 3, 4, 5]
constraints = [
    [2, 1, 2, 1, 1, 10],
    [1, 3, 2, 2, 1, 15]
]

# 変数の定義
model.x = Var(range(N), within=NonNegativeReals)

# 目的関数の定義
model.objective = Objective(expr=sum(coefficients[i] * model.x[i] for i in range(N)), sense=maximize)

# 制約条件の定義
model.constraints = ConstraintList()
for con in constraints:
    model.constraints.add(sum(con[i] * model.x[i] for i in range(N)) <= con[-1])

# CBCソルバーを使用する場合
solver = SolverFactory('cbc', executable='/usr/bin/cbc')
result = solver.solve(model)

# 結果の表示
print("CBC solver results:")
for i in range(N):
    print(f"Optimal value for x[{i}]:", model.x[i]())
print("Optimal objective value:", model.objective())

[実行結果]

CBC solver results:
Optimal value for x[0]: 0.0
Optimal value for x[1]: 0.0
Optimal value for x[2]: 0.0
Optimal value for x[3]: 0.0
Optimal value for x[4]: 10.0
Optimal objective value: 50.0

Pyomoは非常に柔軟で強力なツールであり、さまざまな最適化問題に対応できます。

これらの例を参考に、実際の問題に応用してみてください。

scikit-learn

scikit-learn

scikit-learnは、Pythonで機械学習モデルを構築、評価、実行するためのライブラリです。

以下に、scikit-learnの便利な使い方をいくつか紹介します。

1. 基本的なデータの読み込みと前処理

1
2
3
4
5
6
7
8
9
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split

# データの読み込み
iris = load_iris()
X, y = iris.data, iris.target

# トレーニングセットとテストセットに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

2. モデルの構築とトレーニング

1
2
3
4
5
6
7
from sklearn.ensemble import RandomForestClassifier

# モデルの構築
model = RandomForestClassifier(n_estimators=100, random_state=42)

# モデルのトレーニング
model.fit(X_train, y_train)

3. モデルの評価

1
2
3
4
5
6
7
8
9
10
11
12
from sklearn.metrics import accuracy_score, classification_report

# 予測
y_pred = model.predict(X_test)

# 精度の評価
accuracy = accuracy_score(y_test, y_pred)
print(f'Accuracy: {accuracy}')

# 詳細な評価
report = classification_report(y_test, y_pred)
print(report)

[実行結果]

Accuracy: 1.0
              precision    recall  f1-score   support

           0       1.00      1.00      1.00        10
           1       1.00      1.00      1.00         9
           2       1.00      1.00      1.00        11

    accuracy                           1.00        30
   macro avg       1.00      1.00      1.00        30
weighted avg       1.00      1.00      1.00        30

4. クロスバリデーション

1
2
3
4
5
6
7
8
9
from sklearn.model_selection import cross_val_score

# クロスバリデーションの実行
scores = cross_val_score(model, X, y, cv=5)

# 平均スコアと標準偏差
print(f'Cross-validation scores: {scores}')
print(f'Mean score: {scores.mean()}')
print(f'Standard deviation: {scores.std()}')

[実行結果]

Cross-validation scores: [0.96666667 0.96666667 0.93333333 0.96666667 1.        ]
Mean score: 0.9666666666666668
Standard deviation: 0.02108185106778919

5. ハイパーパラメータチューニング

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from sklearn.model_selection import GridSearchCV

# ハイパーパラメータの候補を設定
param_grid = {
    'n_estimators': [50, 100, 200],
    'max_depth': [None, 10, 20, 30]
}

# グリッドサーチの実行
grid_search = GridSearchCV(model, param_grid, cv=5)
grid_search.fit(X_train, y_train)

# 最良のパラメータとスコア
best_params = grid_search.best_params_
best_score = grid_search.best_score_

print(f'Best parameters: {best_params}')
print(f'Best cross-validation score: {best_score}')

[実行結果]

Best parameters: {'max_depth': None, 'n_estimators': 200}
Best cross-validation score: 0.95

6. パイプラインの作成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# パイプラインの構築
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('classifier', RandomForestClassifier(n_estimators=100, random_state=42))
])

# パイプラインでのトレーニングと評価
pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
print(f'Accuracy: {accuracy}')

[実行結果]

Accuracy: 1.0

7. 主成分分析 (PCA)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from sklearn.decomposition import PCA

# PCAの適用
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

# トレーニングセットとテストセットに分割
X_train, X_test, y_train, y_test = train_test_split(X_pca, y, test_size=0.2, random_state=42)

# モデルのトレーニングと評価
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
print(f'Accuracy with PCA: {accuracy}')

[実行結果]

Accuracy with PCA: 1.0

scikit-learnは、シンプルなAPIで多くの機械学習タスクを簡単に実行できるように設計されています。

これらの基本的な使い方を理解することで、さらに高度な機械学習の実装や実験に役立てることができます。

SciPy

SciPy

SciPyは科学計算やエンジニアリングの分野で広く利用されているPythonライブラリで、便利な機能が多数揃っています。

ここでは、SciPyの便利な使い方をいくつか紹介します。

1. 数値積分

SciPyのintegrateモジュールを使って数値積分を行うことができます。

以下は、関数の定積分を計算する例です。

1
2
3
4
5
6
7
8
9
import scipy.integrate as integrate

# 定義する関数
def f(x):
    return x ** 2

# 数値積分
result, error = integrate.quad(f, 0, 1)
print("積分の結果:", result)

[実行結果]

積分の結果: 0.33333333333333337

2. 最適化

SciPyのoptimizeモジュールを使って関数の最小化や方程式の解を求めることができます。

関数の最小化

1
2
3
4
5
6
7
8
9
10
import numpy as np
import scipy.optimize as optimize

# 定義する関数
def f(x):
    return x ** 2 + 10 * np.sin(x)

# 最小化
result = optimize.minimize(f, x0=0)
print("最適化の結果:", result.x)

[実行結果]

最適化の結果: [-1.30644012]

方程式の解を求める

1
2
3
4
5
6
7
# 定義する方程式
def equation(x):
    return x ** 2 - 2

# 解を求める
solution = optimize.root(equation, x0=1)
print("方程式の解:", solution.x)

[実行結果]

方程式の解: [1.41421356]

3. 線形代数

SciPyのlinalgモジュールを使って行列演算を行うことができます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
import scipy.linalg as linalg
import numpy as np

# 行列の定義
A = np.array([[3, 2], [1, 4]])
B = np.array([1, 2])

# 線形方程式の解を求める
x = linalg.solve(A, B)
print("方程式の解:", x)

# 行列の逆行列を求める
A_inv = linalg.inv(A)
print("逆行列:", A_inv)

[実行結果]

方程式の解: [0.  0.5]
逆行列: [[ 0.4 -0.2]
 [-0.1  0.3]]

4. 信号処理

SciPyのsignalモジュールを使ってフィルタリングや周波数解析を行うことができます。

信号のフィルタリング

1
2
3
4
5
6
7
8
9
10
11
12
from scipy.signal import butter, filtfilt
import numpy as np

# バターワースフィルタの設計
b, a = butter(N=4, Wn=0.2)

# 信号の定義（繰り返し）
signal = np.tile([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 3)  # 繰り返して長さを増やす

# フィルタリング
filtered_signal = filtfilt(b, a, signal)
print("フィルタリング後の信号:", filtered_signal)

[実行結果]

フィルタリング後の信号: [ 0.97688963  2.12759398  3.32584065  4.55768608  5.73180471  6.68699207
  7.23854322  7.25256462  6.72119364  5.80208259  4.79233683  4.03518382
  3.79495244  4.15630622  4.99146245  6.00594272  6.84255021  7.20634692
  6.96867975  6.2130705   5.20271618  4.28015452  3.74274467  3.74971619
  4.29610177  5.25159104  6.43442239  7.6839036   8.90514189 10.07593652]

5. 統計

SciPyのstatsモジュールを使って統計解析を行うことができます。

基本的な統計量の計算

1
2
3
4
5
6
7
8
9
10
11
12
13
from scipy import stats

# データの定義
data = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

# 平均と標準偏差
mean = np.mean(data)
std_dev = np.std(data)
print("平均:", mean, "標準偏差:", std_dev)

# 正規分布のフィッティング
param = stats.norm.fit(data)
print("正規分布のパラメータ:", param)

[実行結果]

平均: 5.5 標準偏差: 2.8722813232690143
正規分布のパラメータ: (5.5, 2.8722813232690143)

6. 補間

SciPyのinterpolateモジュールを使って補間を行うことができます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

# データの定義
x = np.linspace(0, 10, 10)
y = np.sin(x)

# 線形補間
f = interp1d(x, y, kind='linear')

# 新しいデータ点の生成
x_new = np.linspace(0, 10, 100)
y_new = f(x_new)

# プロット
plt.plot(x, y, 'o', label='Original data')
plt.plot(x_new, y_new, '-', label='Interpolated data')
plt.legend()
plt.show()

[実行結果]

これらはSciPyの便利な使い方のほんの一部です。

SciPyは多機能で多岐にわたる分野で使用できるため、目的に応じて必要なモジュールや関数を活用してみてください。

PuLP

PuLP

PuLPはPython用の線形プログラミングライブラリで、最適化問題を簡単に定義し解くことができます。

以下にPuLPの基本的な使い方と便利な機能を紹介します。

基本的な使い方

1. インストール

まず、PuLPをインストールします。

1
pip install pulp

2. 問題の定義

次に、線形プログラミング問題を定義します。

例えば、以下は最大化問題の定義です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import pulp

# 問題の定義
model = pulp.LpProblem("Maximize_Profit", pulp.LpMaximize)

# 変数の定義
x = pulp.LpVariable("x", lowBound=0, cat='Continuous')
y = pulp.LpVariable("y", lowBound=0, cat='Continuous')

# 目的関数の定義
model += 40 * x + 30 * y, "Total Profit"

# 制約条件の定義
model += 2 * x + y <= 20, "Constraint 1"
model += 4 * x + 3 * y <= 45, "Constraint 2"

# 問題を解く
model.solve()

# 結果の表示
print(f"Status: {pulp.LpStatus[model.status]}")
print(f"x: {pulp.value(x)}")
print(f"y: {pulp.value(y)}")
print(f"Total Profit: {pulp.value(model.objective)}")

[実行結果]

Status: Optimal
x: 0.0
y: 15.0
Total Profit: 450.0

応用的な使い方

3. 整数計画問題

PuLPでは整数計画問題も扱うことができます。

以下は整数計画問題の例です。

1
2
3
4
5
# 変数の定義 (整数)
x = pulp.LpVariable("x", lowBound=0, cat='Integer')
y = pulp.LpVariable("y", lowBound=0, cat='Integer')

# 問題の定義、目的関数、制約条件は前述の通り

4. 複数の目的関数

PuLPでは複数の目的関数を扱うことができますが、ここでは重み付きの合成目的関数として実装します。

1
2
# 目的関数の定義 (例: 50% の重みを持つ2つの目的関数)
model += 0.5 * (40 * x + 30 * y) + 0.5 * (30 * x + 20 * y), "Combined Objective"

5. 制約条件の名前付きリスト

制約条件をリストで管理し、名前を付けて定義することができます。

1
2
3
4
5
6
7
constraints = [
    (2 * x + y <= 20, "Raw Material Constraint"),
    (4 * x + 3 * y <= 45, "Labor Constraint")
]

for constraint, name in constraints:
    model += constraint, name

6. コマンドラインソルバーの利用

PuLPは多くのコマンドラインソルバー（例: CBC, Gurobi, CPLEX）と連携できます。

1
2
3
# CBCソルバーを使用
solver = pulp.PULP_CBC_CMD()
model.solve(solver)

実践的な例：供給と需要の最適化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
import pulp

# 問題の定義
model = pulp.LpProblem("Minimize_Transportation_Cost", pulp.LpMinimize)

# サプライチェーンの定義
warehouses = ["A", "B"]
stores = ["X", "Y", "Z"]
supply = {"A": 300, "B": 400}
demand = {"X": 250, "Y": 350, "Z": 100}

# 輸送コストの定義
costs = {"A": {"X": 2, "Y": 4, "Z": 5}, "B": {"X": 3, "Y": 1, "Z": 3}}

# 変数の定義
vars = pulp.LpVariable.dicts("Route", (warehouses, stores), lowBound=0, cat='Continuous')

# 目的関数の定義
model += pulp.lpSum([vars[w][s] * costs[w][s] for w in warehouses for s in stores]), "Total Cost"

# 制約条件の定義
for w in warehouses:
    model += pulp.lpSum([vars[w][s] for s in stores]) <= supply[w], f"Supply_{w}"

for s in stores:
    model += pulp.lpSum([vars[w][s] for w in warehouses]) >= demand[s], f"Demand_{s}"

# 問題を解く
model.solve()

# 結果の表示
print(f"Status: {pulp.LpStatus[model.status]}")
for w in warehouses:
    for s in stores:
        print(f"Route {w} to {s}: {pulp.value(vars[w][s])} units")
print(f"Total Cost: {pulp.value(model.objective)}")

[実行結果]

Status: Optimal
Route A to X: 250.0 units
Route A to Y: 0.0 units
Route A to Z: 50.0 units
Route B to X: 0.0 units
Route B to Y: 350.0 units
Route B to Z: 50.0 units
Total Cost: 1250.0

まとめ

PuLPを使うと、線形計画問題や整数計画問題を簡単に定義し解くことができます。

複雑な制約条件や複数の目的関数も扱えるため、実際のビジネスや研究の最適化問題に広く応用可能です。