効用最大化問題 CVXPY

効用最大化問題 CVXPY

経済学に関する最適化問題の一例として、消費者の効用最大化問題を解いてみましょう。

消費者の効用最大化は、限られた予算の中で効用(満足度)を最大化するためにどのように消費するかを決定する問題です。

問題設定

消費者は2つの財(商品)$ ( x ) $と$ ( y ) $を消費します。

それぞれの財の価格は$ ( p_x ) $と$ ( p_y ) $です。

消費者の予算は$ ( M ) $です。

消費者の効用関数は$ ( U(x, y) = x^{0.5} y^{0.5} ) $というコブ・ダグラス型の効用関数とします。

目標

予算制約を守りながら、消費者の効用を最大化するための財$ ( x ) $と$ ( y ) $の消費量を求めます。

解法

CVXPYを使ってこの最適化問題を解きます。

ステップ1: ライブラリのインストール

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

1
pip install cvxpy

ステップ2: データの準備

次に、価格と予算を設定します。

1
2
3
4
5
6
import cvxpy as cp

# データの設定
px = 20 # 財xの価格
py = 10 # 財yの価格
M = 100 # 予算

ステップ3: 最適化問題の定式化

次に、CVXPYを使って効用最大化問題を定式化します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 変数
x = cp.Variable(nonneg=True) # 財xの消費量
y = cp.Variable(nonneg=True) # 財yの消費量

# 効用関数
utility = cp.sqrt(x) * cp.sqrt(y)

# 最適化問題の定式化
objective = cp.Maximize(utility) # 効用を最大化
constraints = [
px * x + py * y <= M # 予算制約
]

# 問題の作成と解決
problem = cp.Problem(objective, constraints)
problem.solve(qcp=True)

# 最適な消費量
print("最適な財xの消費量:", x.value)
print("最適な財yの消費量:", y.value)
print("最大化された効用:", utility.value)

ステップ4: 結果の確認

最適な消費量最大化された効用を確認します。

1
2
3
print("最適な財xの消費量:", x.value)
print("最適な財yの消費量:", y.value)
print("最大化された効用:", (x.value**0.5) * (y.value**0.5))

まとめ

この例では、消費者の効用最大化問題を解きました。

CVXPYを使うことで、経済学に関するさまざまな最適化問題を解くことができます。

他にも以下のような問題に利用できます。

  • 生産者の利益最大化
  • 労働市場の均衡
  • 経済成長モデルの最適化
  • 政策シミュレーション

経済学の問題において最適化手法は非常に重要な役割を果たしており、CVXPYはそのための強力なツールとなります。

結果解説

[実行結果]

最適な財xの消費量: 2.4999980742405747
最適な財yの消費量: 5.000003845908935
最大化された効用: 3.5355339039482856
最適な財xの消費量: 2.4999980742405747
最適な財yの消費量: 5.000003845908935
最大化された効用: 3.5355339039482856

この結果は、予算制約のもとで消費者の効用を最大化するための最適な財$ ( x ) $と$ ( y ) $の消費量を示しています。

結果の詳細は次の通りです:

最適な消費量

  • 財 ( x ) の消費量: 2.4999980742405747
  • 財 ( y ) の消費量: 5.000003845908935

これらの数値は、消費者が最大限の満足度を得るために、それぞれの財をどの程度消費すべきかを示しています。

最大化された効用

  • 最大化された効用: 3.5355339039482856

この数値は、与えられた予算内で達成できる最大の効用(満足度)を表しています。

結果の解説

この問題では、消費者が以下の効用関数を最大化しようとしています:
$$
U(x, y) = x^{0.5} y^{0.5}
$$

また、予算制約は次のようになります:
$$
p_x \cdot x + p_y \cdot y \leq M
$$

具体的には、以下の条件のもとで解を求めました:

  • 財$ ( x ) $の価格$ ( p_x ) = 20$
  • 財$ ( y ) $の価格$ ( p_y ) = 10$
  • 予算$ ( M ) = 100$

結果の確認

結果から、予算制約の範囲内で以下のような消費を行うと、効用が最大化されることがわかります:

  • 財$ ( x ) $を約$2.5$単位消費する
  • 財$ ( y ) $を約$5.0$単位消費する

このとき、消費者の満足度(効用)は約$3.54$となります。
消費者が与えられた予算の中で、効用を最大化するためにどのような消費を行うべきかが明確になりました。

最適な消費量効用は、モデルの条件下で最適化されたものであり、経済学の理論に基づいた合理的な消費行動を示しています。

住宅価格予想 statsmodels

住宅価格予想 statsmodels

Pythonのstatsmodelsライブラリを使用して、現実的なデータ分析の問題を解決する例を提供します。

以下に、住宅価格のデータセットを使用して、価格を予測するための線形回帰モデルを構築する例を示します。

例題

問題: 住宅の面積、部屋数、築年数を使って、住宅価格を予測する線形回帰モデルを構築してください。

ステップ

  1. データを読み込む
  2. データの前処理を行う
  3. モデルを構築する
  4. モデルを評価する

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

まず、必要なライブラリをインストールします。これにはpandasstatsmodelsが必要です。

1
pip install pandas statsmodels

コード

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
import pandas as pd
import statsmodels.api as sm

# データの読み込み
# 例として、以下のようなCSVファイルを使うとします
# 面積 (平方フィート), 部屋数, 築年数, 価格 (USD)
data = {
'area': [1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400],
'rooms': [3, 3, 3, 4, 4, 4, 5, 5, 5, 5],
'age': [10, 15, 20, 25, 5, 10, 15, 20, 25, 5],
'price': [300000, 320000, 340000, 360000, 380000, 400000, 420000, 440000, 460000, 480000]
}

df = pd.DataFrame(data)

# 独立変数 (特徴量) と従属変数 (ターゲット) の設定
X = df[['area', 'rooms', 'age']]
y = df['price']

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

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

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

# 新しいデータで予測
new_data = pd.DataFrame({'const': 1, 'area': [2500], 'rooms': [4], 'age': [10]})
prediction = model.predict(new_data)

print(f"予測された価格: ${prediction[0]:,.2f}")

説明

  1. データの読み込み: データは辞書として定義され、pandasのデータフレームに変換します。
  2. 特徴量とターゲットの設定: 面積、部屋数、築年数を特徴量として、価格をターゲットとします。
  3. 定数項の追加: statsmodelsでは、定数項(切片)を明示的に追加する必要があります。
  4. モデルの構築: OLS (最小二乗法) を使用して回帰モデルを構築し、モデルをフィッティングします。
  5. モデルの評価: summary メソッドでモデルの概要を表示します。
  6. 予測: 新しいデータポイントで予測を行います。

この例では、住宅価格のデータを基にして線形回帰モデルを構築し、モデルの概要を出力します。また、新しいデータポイントに基づいて住宅価格を予測します。

これにより、statsmodelsを使った基本的な回帰分析の方法が理解できます。

ソースコード解説

以下のPythonコードは、statsmodelsライブラリを使用して線形回帰モデルを構築し、住宅価格を予測するプロセスを示しています。

コードの各章立てについて詳しく説明します。

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

1
2
import pandas as pd
import statsmodels.api as sm

ここでは、データ操作のためにpandasライブラリ、統計モデリングのためにstatsmodelsライブラリをインポートしています。

2. データの読み込み

1
2
3
4
5
6
7
data = {
'area': [1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400],
'rooms': [3, 3, 3, 4, 4, 4, 5, 5, 5, 5],
'age': [10, 15, 20, 25, 5, 10, 15, 20, 25, 5],
'price': [300000, 320000, 340000, 360000, 380000, 400000, 420000, 440000, 460000, 480000]
}
df = pd.DataFrame(data)

ここでは、住宅の面積、部屋数、築年数、および価格に関するデータを持つ辞書を定義し、それをpandasのデータフレームに変換しています。

3. 独立変数 (特徴量) と従属変数 (ターゲット) の設定

1
2
X = df[['area', 'rooms', 'age']]
y = df['price']

この部分では、住宅価格を予測するために使用する特徴量(面積、部屋数、築年数)をXとして、ターゲット変数(価格)をyとして設定しています。

4. 定数項を追加

1
X = sm.add_constant(X)

statsmodelsの線形回帰モデルでは、データに定数項(バイアス)を追加する必要があります。add_constant関数を使用して、特徴量データフレームに定数項を追加しています。

5. モデルの構築

1
model = sm.OLS(y, X).fit()

ここでは、OLS(Ordinary Least Squares: 最小二乗法)を使用して線形回帰モデルを構築し、データにフィットさせています。

6. モデルの概要を表示

1
print(model.summary())

model.summary()を使って、回帰モデルの詳細な結果を表示しています。この結果には、決定係数(R-squared)、回帰係数、p値などの統計情報が含まれています。

7. 新しいデータで予測

1
2
3
4
new_data = pd.DataFrame({'const': 1, 'area': [2500], 'rooms': [4], 'age': [10]})
prediction = model.predict(new_data)

print(f"予測された価格: ${prediction[0]:,.2f}")

最後に、新しいデータポイント(面積: 2500平方フィート、部屋数: 4、築年数: 10)を使って価格を予測しています。
この新しいデータには、前に追加した定数項も含める必要があります。
model.predictを使って価格を予測し、その結果をフォーマットして表示しています。

解説

  • 定数項の追加: 線形回帰モデルでは、バイアス項(定数項)を追加することで、モデルの予測が正確になることが多いです。
    これは、入力特徴量がすべてゼロの場合でも一定の予測値を出力するために必要です。
  • OLSの使用: OLSは、線形回帰モデルをフィットさせるための一般的な方法です。
    ターゲット変数と特徴量の関係を線形に近似します。
  • モデル評価: summary()関数は、モデルのパフォーマンスを評価するための詳細な統計情報を提供します。
    これにより、モデルの有効性や各特徴量の影響を評価できます。
  • 予測: 新しいデータポイントを使って予測を行う際には、モデルに適用したのと同じ前処理(定数項の追加など)を行う必要があります。

このコードは、基本的な線形回帰モデルの構築と予測のプロセスを示しています。
データの前処理からモデルの評価、そして新しいデータに対する予測までの一連の流れを学ぶのに役立ちます。

結果解説

[実行結果]

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  price   R-squared:                       1.000
Model:                            OLS   Adj. R-squared:                  1.000
Method:                 Least Squares   F-statistic:                 2.164e+30
Date:                Sat, 13 Jul 2024   Prob (F-statistic):           1.73e-90
Time:                        01:29:56   Log-Likelihood:                 222.01
No. Observations:                  10   AIC:                            -436.0
Df Residuals:                       6   BIC:                            -434.8
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const      -2.328e-10   1.99e-10     -1.171      0.286   -7.19e-10    2.54e-10
area         200.0000   2.43e-13   8.22e+14      0.000     200.000     200.000
rooms       5.821e-11   8.44e-11      0.690      0.516   -1.48e-10    2.65e-10
age        -1.819e-12    3.3e-12     -0.552      0.601   -9.89e-12    6.25e-12
==============================================================================
Omnibus:                        0.133   Durbin-Watson:                   0.667
Prob(Omnibus):                  0.936   Jarque-Bera (JB):                0.294
Skew:                           0.198   Prob(JB):                        0.863
Kurtosis:                       2.258   Cond. No.                     1.80e+04
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.8e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
予測された価格: $500,000.00

この出力結果は、statsmodelsを用いて実行した線形回帰モデルの概要を示しています。

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

モデル全体の評価

  1. R-squared (決定係数): 1.000

    • モデルがターゲット変数(価格)の分散をどれだけ説明しているかを示します。

    1.000は、モデルがデータを完全に説明していることを示します。

  2. Adj. R-squared (自由度調整済み決定係数): 1.000

    • 自由度を考慮して調整された決定係数です。
      これも1.000であり、モデルが非常に良くフィットしていることを示しています。
  3. F-statistic: 2.164e+30

    • モデル全体の有意性を検定するための統計量です。
      非常に大きな値であり、モデルが有意であることを示しています。
  4. Prob (F-statistic): 1.73e-90

    • F検定のp値です。
      非常に小さい値であり、モデル全体が有意であることを示しています。

個々の係数の評価

  1. const (定数項): -2.328e-10 (標準誤差: 1.99e-10, p値: 0.286)

    • 定数項の値です。
      p値が0.286であり、統計的に有意ではないことを示しています。
  2. area (面積): 200.0000 (標準誤差: 2.43e-13, p値: 0.000)

    • 面積の係数です。
      係数が200であり、非常に低い標準誤差と0.000のp値から、面積が価格に非常に強い影響を持つことが分かります。
  3. rooms (部屋数): 5.821e-11 (標準誤差: 8.44e-11, p値: 0.516)

    • 部屋数の係数です。
      p値が0.516であり、統計的に有意ではないことを示しています。
  4. age (築年数): -1.819e-12 (標準誤差: 3.3e-12, p値: 0.601)

    • 築年数の係数です。
      p値が0.601であり、統計的に有意ではないことを示しています。

モデルの統計量

  1. Omnibus: 0.133 (Prob Omnibus: 0.936)

    • 残差の正規性を検定するための統計量です。
      p値が0.936であり、残差が正規分布に従っていることを示しています。
  2. Durbin-Watson: 0.667

    • 残差の自己相関を検定するための統計量です。
      0に近いと正の自己相関があることを示し、2に近いと自己相関がないことを示します。
      ここでは0.667であり、正の自己相関がある可能性があります。
  3. Jarque-Bera (JB): 0.294 (Prob JB: 0.863)

    • 残差の正規性を検定するための統計量です。
      p値が0.863であり、残差が正規分布に従っていることを示しています。
  4. Skew: 0.198

    • 残差の歪度を示します。
      0に近い値は対称性を示します。
  5. Kurtosis: 2.258

    • 残差の尖度を示します。
      3に近い値は正規分布に近いことを示します。

注意点

  1. Condition Number (条件数): 1.8e+04
    • 高い条件数は、多重共線性やその他の数値的な問題がある可能性を示唆します。
      ここでは非常に高い値であり、多重共線性が存在する可能性があります。

予測された価格

最後に、モデルを使って新しいデータポイントの価格を予測しました。
結果として、予測された価格は$500,000.00でした。

この結果は、面積が価格に強く影響していることを示していますが、部屋数や築年数の影響は有意ではない可能性があります。
また、高い条件数はデータの多重共線性の問題を示唆しているため、注意が必要です。

グラフィカル診断 Statsmodels

グラフィカル診断 Statsmodels

Statsmodelsは、Pythonの強力な統計モデリングツールです。

以下に、Statsmodelsグラフィカル診断を紹介します。

グラフィカル診断

モデルの診断には様々なプロットを使用できます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import statsmodels.api as sm
import numpy as np

# サンプルデータの作成
data = sm.datasets.get_rdataset('mtcars').data

# モデルの作成
model = sm.OLS(data['mpg'], sm.add_constant(data[['hp', 'wt']])).fit()

# 残差プロット
sm.graphics.plot_regress_exog(model, 'hp')
sm.graphics.plot_regress_exog(model, 'wt')

# QQプロット
sm.qqplot(model.resid, line='s')

[実行結果]


ソースコード解説

ソースコードの詳細な説明を行います。

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

1
2
import statsmodels.api as sm
import numpy as np
  • statsmodels.api: 回帰分析や統計モデルのライブラリで、データ分析とモデル評価に広く使われます。
  • numpy: 数値計算のためのライブラリで、配列操作や数値演算を簡単に行えます。

2. サンプルデータの作成

1
data = sm.datasets.get_rdataset('mtcars').data
  • sm.datasets.get_rdataset('mtcars'): Rのデータセットを取得するメソッド。ここでは有名なmtcarsデータセットを使用しています。
  • data: mtcarsデータセットをデータフレームとして取得。

3. モデルの作成

1
model = sm.OLS(data['mpg'], sm.add_constant(data[['hp', 'wt']])).fit()
  • sm.OLS(endog, exog): OLS(最小二乗法)回帰モデルを作成するための関数。
    • endog: 被説明変数(目的変数)、ここでは燃費(mpg)。
    • exog: 説明変数(特徴量)、ここでは馬力(hp)と車両重量(wt)。
    • sm.add_constant(data[['hp', 'wt']]): 説明変数に定数項(切片)を追加。
  • fit(): モデルをデータに適合させるメソッド。

4. 残差プロット

1
2
sm.graphics.plot_regress_exog(model, 'hp')
sm.graphics.plot_regress_exog(model, 'wt')
  • sm.graphics.plot_regress_exog(model, exog_idx): 指定した説明変数に対する残差プロットを作成。
    • model: 作成した回帰モデル。
    • exog_idx: 残差プロットを作成する説明変数(ここではhpwt)。

これにより、以下のような4つのグラフが生成されます。

  • Y and Fitted vs. X: 実測値と予測値のプロット。
  • Residuals versus X: 残差と説明変数のプロット。
  • Partial regression plot: 他の変数を固定した状態での説明変数と目的変数の関係を示すプロット。
  • CCPR Plot: 説明変数の影響を考慮した残差プロット。

5. QQプロット

1
sm.qqplot(model.resid, line='s')
  • sm.qqplot(data, line='s'): Q-Qプロットを作成するための関数。
    • data: プロットするデータ、ここでは回帰モデルの残差(model.resid)。
    • line='s': プロットに標準正規分布の参照線を追加。

まとめ

このコードは、mtcarsデータセットを用いて、燃費(mpg馬力(hp車両重量(wtから予測する回帰モデルを作成し、モデルの評価を行っています。

具体的には、以下のことを行っています。

  1. mtcarsデータセットを取得。
  2. mpgを目的変数、hpwtを説明変数とするOLS回帰モデルを作成し、適合。
  3. 説明変数に対する残差プロットを作成し、モデルの適合性を評価
  4. 残差のQ-Qプロットを作成し、残差が正規分布に従っているかを確認。

これにより、モデルの性能適合性残差の特性を視覚的に評価することができます。

グラフ解説

表示されるグラフを解説します。


[実行結果1]

Q-Q プロット

  • 説明: Q-Qプロットは、データの分布が理論的な分布(通常は正規分布)にどれだけ近いかを視覚化するためのプロットです。
  • 見方: 点が赤い直線に近いほど、データが正規分布に従っていることを示しています。
  • 解釈: 多くの点が直線上にあるため、データは正規分布に近いと考えられます。ただし、端のほうでは少し外れている点もあり、これがデータの外れ値や分布の尾部の特性を示しています。

[実行結果2]

回帰プロット(hpに対する)

  • Y and Fitted vs. X: 実測値 (mpg) と予測値 (fitted) のプロット。予測値が実測値にどれだけ近いかを示しています。
  • Residuals versus hp: 残差(実測値と予測値の差)とhpのプロット。残差がランダムに散らばっている場合、モデルは適切です。
  • Partial regression plot: 他の変数を一定とした場合のhpとmpgの関係を示すプロット。
  • CCPR Plot: 残差プロットに相関を持たせたもの。各観測点の影響を視覚化しています。

解釈

  • Y and Fitted vs. X: 多くの点が一致しているため、モデルの予測は概ね正確です。
  • Residuals versus hp: 残差はランダムに分布しており、特定のパターンがないため、モデルは適切です。
  • Partial regression plot: hpが増加するにつれてmpgが減少していることを示しています。
  • CCPR Plot: hpの影響が残差に与える影響を示しており、負の関係が見られます。

[実行結果3]

回帰プロット(wtに対する)

  • Y and Fitted vs. X: 実測値 (mpg) と予測値 (fitted) のプロット。予測値が実測値にどれだけ近いかを示しています。
  • Residuals versus wt: 残差とwtのプロット。残差がランダムに散らばっている場合、モデルは適切です。
  • Partial regression plot: 他の変数を一定とした場合のwtとmpgの関係を示すプロット。
  • CCPR Plot: 残差プロットに相関を持たせたもの。各観測点の影響を視覚化しています。

解釈

  • Y and Fitted vs. X: 実測値と予測値が近いため、モデルの予測は正確です。
  • Residuals versus wt: 残差がランダムに分布しており、特定のパターンがないため、モデルは適切です。
  • Partial regression plot: wtが増加するにつれてmpgが減少していることを示しています。
  • CCPR Plot: wtの影響が残差に与える影響を示しており、負の関係が見られます。

これらのプロットは、回帰モデルがどの程度データに適合しているか、残差がランダムに分布しているか、特定の変数の影響を視覚化するための有用なツールです。

python-constraint 学校の時間割

python-constraint 学校の時間割

以下の問題を解いてみましょう。

問題

学校の時間割

ある学校では、次の科目を月曜日から金曜日の$5$日間の授業時間に割り当てる必要があります。
各科目は$1$回ずつ行われ、各科目の授業は$1$日1時間だけです。

  • 数学(Math)
  • 英語(English)
  • 理科(Science)
  • 社会(SocialStudies)
  • 体育(PhysicalEducation)

条件:

  • 数学は月曜日か火曜日に行われる必要があります。
  • 英語は水曜日に行われる必要があります。
  • 体育は金曜日に行われる必要があります。

この条件に従って、時間割を作成してください。

解法

Python-Constraintライブラリを使って、条件に合った時間割を求めます。

コード

以下のコードを実行して、解を求めます。

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
from constraint import Problem, AllDifferentConstraint

# 科目のリスト
subjects = ['Math', 'English', 'Science', 'SocialStudies', 'PhysicalEducation']

# 曜日のリスト
days = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday']

# 問題を作成
problem = Problem()

# 各曜日に科目を割り当てる
for day in days:
problem.addVariable(day, subjects)

# 各曜日には異なる科目が割り当てられるように制約を追加
problem.addConstraint(AllDifferentConstraint(), days)

# 数学は月曜日か火曜日に行われる必要がある
def math_on_monday_or_tuesday(monday, tuesday):
return monday == 'Math' or tuesday == 'Math'

# 英語は水曜日に行われる必要がある
def english_on_wednesday(wednesday):
return wednesday == 'English'

# 体育は金曜日に行われる必要がある
def pe_on_friday(friday):
return friday == 'PhysicalEducation'

# 制約を追加
problem.addConstraint(math_on_monday_or_tuesday, ['Monday', 'Tuesday'])
problem.addConstraint(english_on_wednesday, ['Wednesday'])
problem.addConstraint(pe_on_friday, ['Friday'])

# 解を求める
solutions = problem.getSolutions()

# 結果を表示
if solutions:
for solution in solutions:
print(solution)
else:
print("解が見つかりませんでした。")

このコードでは、指定された条件に従って$5$日間の学校の時間割を作成します。

結果として、各曜日にどの科目を割り当てるかが表示されます。

ソースコード解説

以下は、ソースコードの詳細な説明です。

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

1
from constraint import Problem, AllDifferentConstraint
  • constraint ライブラリから ProblemAllDifferentConstraint をインポートしています。
    • Problem は制約充足問題を定義するためのクラスです。
    • AllDifferentConstraint はすべての変数に異なる値を割り当てる制約を定義するためのクラスです。

2. 科目と曜日のリスト定義

1
2
3
4
5
# 科目のリスト
subjects = ['Math', 'English', 'Science', 'SocialStudies', 'PhysicalEducation']

# 曜日のリスト
days = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday']
  • subjects は授業科目のリストです。
  • days は曜日のリストです。

3. 問題の作成

1
2
# 問題を作成
problem = Problem()
  • Problem クラスのインスタンスを作成し、問題を定義します。

4. 変数の追加

1
2
3
# 各曜日に科目を割り当てる
for day in days:
problem.addVariable(day, subjects)
  • 各曜日に対して、可能な科目を変数として追加します。
  • addVariable メソッドを使用して、各曜日(変数)が subjects のいずれかの科目(値)を取ることを指定します。

5. 制約の追加

1
2
# 各曜日には異なる科目が割り当てられるように制約を追加
problem.addConstraint(AllDifferentConstraint(), days)
  • AllDifferentConstraint を使用して、すべての曜日に異なる科目が割り当てられるように制約を追加します。

6. 特定の制約を定義

1
2
3
4
5
6
7
8
9
10
11
# 数学は月曜日か火曜日に行われる必要がある
def math_on_monday_or_tuesday(monday, tuesday):
return monday == 'Math' or tuesday == 'Math'

# 英語は水曜日に行われる必要がある
def english_on_wednesday(wednesday):
return wednesday == 'English'

# 体育は金曜日に行われる必要がある
def pe_on_friday(friday):
return friday == 'PhysicalEducation'
  • 特定の制約条件を関数として定義します。
    • math_on_monday_or_tuesday は、数学が月曜日火曜日に行われる制約です。
    • english_on_wednesday は、英語が水曜日に行われる制約です。
    • pe_on_friday は、体育が金曜日に行われる制約です。

7. 制約の追加

1
2
3
4
# 制約を追加
problem.addConstraint(math_on_monday_or_tuesday, ['Monday', 'Tuesday'])
problem.addConstraint(english_on_wednesday, ['Wednesday'])
problem.addConstraint(pe_on_friday, ['Friday'])
  • 定義した制約関数を問題に追加します。
    • math_on_monday_or_tuesday月曜日火曜日に対して追加します。
    • english_on_wednesday水曜日に対して追加します。
    • pe_on_friday金曜日に対して追加します。

8. 解の求解

1
2
# 解を求める
solutions = problem.getSolutions()
  • getSolutions メソッドを使用して、問題のすべての解を求めます。

9. 結果の表示

1
2
3
4
5
6
# 結果を表示
if solutions:
for solution in solutions:
print(solution)
else:
print("解が見つかりませんでした。")
  • 解が見つかった場合、すべての解を表示します。
  • 解が見つからなかった場合、その旨を表示します。

まとめ

このコードは、制約充足問題を用いて学校の時間割を作成するものです。

与えられた制約(曜日ごとに異なる科目を割り当て、特定の曜日に特定の科目を割り当てる)を満たすすべての解を求め、表示します。

結果解説

[実行結果]

{'Monday': 'SocialStudies', 'Tuesday': 'Math', 'Friday': 'PhysicalEducation', 'Wednesday': 'English', 'Thursday': 'Science'}
{'Monday': 'Science', 'Tuesday': 'Math', 'Friday': 'PhysicalEducation', 'Wednesday': 'English', 'Thursday': 'SocialStudies'}
{'Monday': 'Math', 'Tuesday': 'SocialStudies', 'Friday': 'PhysicalEducation', 'Wednesday': 'English', 'Thursday': 'Science'}
{'Monday': 'Math', 'Tuesday': 'Science', 'Friday': 'PhysicalEducation', 'Wednesday': 'English', 'Thursday': 'SocialStudies'}

これは、学校の時間割問題に対する解がいくつか見つかったことを示しています。

各解は、月曜日から金曜日までの各曜日にどの科目が割り当てられているかを示しています。

各解は異なる時間割を示していますが、いずれも指定された制約を満たしています。

それぞれの解について説明します:

  1. 解1:

    • 月曜日: 社会(SocialStudies)
    • 火曜日: 数学(Math)
    • 水曜日: 英語(English)
    • 木曜日: 理科(Science)
    • 金曜日: 体育(PhysicalEducation)
  2. 解2:

    • 月曜日: 理科(Science)
    • 火曜日: 数学(Math)
    • 水曜日: 英語(English)
    • 木曜日: 社会(SocialStudies)
    • 金曜日: 体育(PhysicalEducation)
  3. 解3:

    • 月曜日: 数学(Math)
    • 火曜日: 社会(SocialStudies)
    • 水曜日: 英語(English)
    • 木曜日: 理科(Science)
    • 金曜日: 体育(PhysicalEducation)
  4. 解4:

    • 月曜日: 数学(Math)
    • 火曜日: 理科(Science)
    • 水曜日: 英語(English)
    • 木曜日: 社会(SocialStudies)
    • 金曜日: 体育(PhysicalEducation)

これらの解は、以下の制約をすべて満たしています:

  • 数学は月曜日か火曜日に行われる。
  • 英語は水曜日に行われる。
  • 体育は金曜日に行われる。

このように、制約を満たしつつ複数の異なる時間割が存在することがわかります。

python-constraint 魔法陣

python-constraint 魔法陣

魔法陣を解く問題を考えます。

魔法陣とは、整数の$ ( n \times n ) $の正方行列で、行・列および対角線の各要素の和がすべて同じになるものです。

魔法陣問題の解決

以下の例では、$3\times3$の魔法陣を解きます。

$3\times3$の魔法陣の各行、列、および対角線の和は$15$である必要があります。

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
47
48
49
import constraint

# 魔法陣のサイズ
n = 3

# 魔法陣の各セルに変数を割り当てる
problem = constraint.Problem()

# 各セルは1から9の整数を持つ
problem.addVariables(range(n * n), range(1, n * n + 1))

# 各セルの値は一意である必要がある
problem.addConstraint(constraint.AllDifferentConstraint())

# 行、列、対角線の和がすべて等しくなる制約を追加
magic_sum = n * (n * n + 1) // 2

# 行の制約
for row in range(n):
problem.addConstraint(
lambda a, b, c: a + b + c == magic_sum,
[row * n + i for i in range(n)]
)

# 列の制約
for col in range(n):
problem.addConstraint(
lambda a, b, c: a + b + c == magic_sum,
[col + n * i for i in range(n)]
)

# 対角線の制約
problem.addConstraint(
lambda a, b, c: a + b + c == magic_sum,
[i * (n + 1) for i in range(n)]
)
problem.addConstraint(
lambda a, b, c: a + b + c == magic_sum,
[(i + 1) * (n - 1) for i in range(n)]
)

# 解を取得
solutions = problem.getSolutions()

# 結果を表示
for solution in solutions:
for row in range(n):
print([solution[row * n + col] for col in range(n)])
print()

説明

  1. 問題定義:

    • 魔法陣のサイズ n = 3 を定義します。
    • 各セルに変数を割り当て、各変数は$1$から$9$までの値を取ることができます。
  2. 制約の追加:

    • constraint.AllDifferentConstraint() を使用して、各セルの値が一意であることを保証します。
    • 各行、列、および対角線の和が magic_sum になるように制約を追加します。
  3. 解の取得と表示:

    • problem.getSolutions() を使用してすべての解を取得し、各解を表示します。

このコードは$3\times3$の魔法陣を解くもので、すべての解を表示します。

他のサイズの魔法陣に対しても、同様の方法で解くことができます。

ただし、サイズが大きくなると計算時間が増加するため注意が必要です。

結果解説

[実行結果]

[6, 1, 8]
[7, 5, 3]
[2, 9, 4]

[6, 7, 2]
[1, 5, 9]
[8, 3, 4]

[8, 1, 6]
[3, 5, 7]
[4, 9, 2]

[8, 3, 4]
[1, 5, 9]
[6, 7, 2]

[4, 3, 8]
[9, 5, 1]
[2, 7, 6]

[4, 9, 2]
[3, 5, 7]
[8, 1, 6]

[2, 9, 4]
[7, 5, 3]
[6, 1, 8]

[2, 7, 6]
[9, 5, 1]
[4, 3, 8]

これは、$3\times3$の魔法陣のすべての可能な解を表示した結果です。

$3\times3$の魔法陣は、各行、各列、そして$2$つの対角線の数の合計がすべて$15$になるような配置です。

以下に各魔法陣の配置を解説します。

  1. 最初の魔法陣:

    1
    2
    3
    [6, 1, 8]
    [7, 5, 3]
    [2, 9, 4]
    • 各行の和: $6+1+8=15$, $7+5+3=15$, $2+9+4=15$
    • 各列の和: $6+7+2=15$, $1+5+9=15$, $8+3+4=15$
    • 対角線の和: $6+5+4=15$, $8+5+2=15$
  2. 二番目の魔法陣:

    1
    2
    3
    [6, 7, 2]
    [1, 5, 9]
    [8, 3, 4]
    • 各行の和: $6+7+2=15$, $1+5+9=15$, $8+3+4=15$
    • 各列の和: $6+1+8=15$, $7+5+3=15$, $2+9+4=15$
    • 対角線の和: $6+5+4=15$, $2+5+8=15$
  3. 三番目の魔法陣:

    1
    2
    3
    [8, 1, 6]
    [3, 5, 7]
    [4, 9, 2]
    • 各行の和: $8+1+6=15$, $3+5+7=15$, $4+9+2=15$
    • 各列の和: $8+3+4=15$, $1+5+9=15$, $6+7+2=15$
    • 対角線の和: $8+5+2=15$, $6+5+4=15$

これらの結果は、すべて$3\times3$の魔法陣であり、どの行、列、および対角線の和もすべて$15$になります。

このように、魔法陣の問題は数学的な整合性を満たすように数値を配置するパズルです。

python-constraintの使い方

python-constraint

python-constraint は、制約プログラミングのライブラリで、制約を満たす解を探索するために使用されます。

ここでは、python-constraint を使用した基本的な使い方を説明します。

基本的な使い方

インストール

まず、python-constraint をインストールします。

1
pip install python-constraint

例題

以下は、基本的な制約プログラムの例です。

この例では、変数 xy を設定し、それらに特定の制約を課して解を探索します。

  1. 変数 xy が $1$ から $10$ までの整数である。
  2. 変数 xy の和が $15$ である。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from constraint import *

# 問題を作成
problem = Problem()

# 変数を追加
problem.addVariable("x", range(1, 11))
problem.addVariable("y", range(1, 11))

# 制約を追加(和が15になるように)
def constraint_sum(x, y):
return x + y == 15

problem.addConstraint(constraint_sum, ["x", "y"])

# 解を取得
solutions = problem.getSolutions()

# 結果を表示
for solution in solutions:
print(solution)

[実行結果]

{'x': 10, 'y': 5}
{'x': 9, 'y': 6}
{'x': 8, 'y': 7}
{'x': 7, 'y': 8}
{'x': 6, 'y': 9}
{'x': 5, 'y': 10}

このスクリプトは、xy が $1$ から $10$ の範囲にあり、それらの和が $15$ であるすべての解を探します。

より複雑な例

次に、より複雑な例として、異なる色でグラフの頂点を塗り分ける問題を解きます。

この例では、4つの頂点$ (A, B, C, D) $があり、隣接する頂点は異なる色で塗る必要があります。

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

# 問題を作成
problem = Problem()

# 変数を追加(各頂点に色を割り当てる)
colors = ["red", "green", "blue"]
vertices = ["A", "B", "C", "D"]
for vertex in vertices:
problem.addVariable(vertex, colors)

# 隣接する頂点の制約を追加(異なる色で塗る)
problem.addConstraint(AllDifferentConstraint(), ["A", "B"])
problem.addConstraint(AllDifferentConstraint(), ["A", "C"])
problem.addConstraint(AllDifferentConstraint(), ["B", "D"])
problem.addConstraint(AllDifferentConstraint(), ["C", "D"])

# 解を取得
solutions = problem.getSolutions()

# 結果を表示
for solution in solutions:
print(solution)

[実行結果]

{'A': 'blue', 'B': 'green', 'C': 'green', 'D': 'blue'}
{'A': 'blue', 'B': 'green', 'C': 'green', 'D': 'red'}
{'A': 'blue', 'B': 'green', 'C': 'red', 'D': 'blue'}
{'A': 'blue', 'B': 'red', 'C': 'green', 'D': 'blue'}
{'A': 'blue', 'B': 'red', 'C': 'red', 'D': 'green'}
{'A': 'blue', 'B': 'red', 'C': 'red', 'D': 'blue'}
{'A': 'green', 'B': 'blue', 'C': 'blue', 'D': 'red'}
{'A': 'green', 'B': 'blue', 'C': 'blue', 'D': 'green'}
{'A': 'green', 'B': 'blue', 'C': 'red', 'D': 'green'}
{'A': 'green', 'B': 'red', 'C': 'blue', 'D': 'green'}
{'A': 'green', 'B': 'red', 'C': 'red', 'D': 'blue'}
{'A': 'green', 'B': 'red', 'C': 'red', 'D': 'green'}
{'A': 'red', 'B': 'green', 'C': 'green', 'D': 'red'}
{'A': 'red', 'B': 'green', 'C': 'green', 'D': 'blue'}
{'A': 'red', 'B': 'green', 'C': 'blue', 'D': 'red'}
{'A': 'red', 'B': 'blue', 'C': 'green', 'D': 'red'}
{'A': 'red', 'B': 'blue', 'C': 'blue', 'D': 'green'}
{'A': 'red', 'B': 'blue', 'C': 'blue', 'D': 'red'}

このスクリプトは、4つの頂点があり、それぞれが3つの色(赤、緑、青)のいずれかに塗られるようにし、隣接する頂点が異なる色で塗られるすべての解を探します。

その他の制約

python-constraint には他にもさまざまな制約を追加する方法があります。

例えば、以下のような制約があります。

  • AllDifferentConstraint(): すべての変数が異なる値を持つようにする。
  • ExactSumConstraint(total): 変数の和が特定の値に等しくなるようにする。
  • MaxSumConstraint(max_sum): 変数の和が特定の最大値を超えないようにする。

これらの制約を適用することで、より複雑な問題を解くことができます。

まとめ

python-constraint は、制約プログラミングを簡単に行うための強力なツールです。

変数と制約を定義し、getSolutions() メソッドを使用して解を探索します。

簡単な例から始めて、徐々に複雑な問題に取り組むことで、python-constraint の使い方に習熟できます。

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%$以上となり、かつリスク(分散)が最小化されます。

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