株価分析 Prophet

実際の株価データを使ったProphetによる分析の具体例を示します。

ここでは、Yahoo Financeから取得したApple Inc.（AAPL）の過去の株価データを使って説明します。

手順

1. データの準備:

   • Yahoo Financeなどからデータを取得する（ここでは、pandas_datareaderを使用します）。
   • データの前処理を行う。

2. Prophetのインストールとインポート:

   1 2 3 4 5

   !pip install yfinance prophet import yfinance as yf from prophet import Prophet import pandas as pd import matplotlib.pyplot as plt

3. データの取得とフォーマット:

   • Yahoo FinanceからAppleの株価データを取得し、Prophetが理解できるフォーマット（ds列に日付、y列に値）に変換します。

     1 2 3 4 5

     # Appleの株価データを取得 df = yf.download('AAPL', start='2020-01-01', end='2023-01-01') df.reset_index(inplace=True) df = df[['Date', 'Close']] df.columns = ['ds', 'y']

4. モデルの作成とトレーニング:

   1 2	model = Prophet() model.fit(df)

5. 予測の作成:

   • 予測する期間のデータフレームを作成し、モデルに予測を行わせます。

     1 2	future = model.make_future_dataframe(periods=365) # 365日先まで予測 forecast = model.predict(future)

6. 予測結果の可視化:

   1 2	fig = model.plot(forecast) plt.show()

以下に、上記手順に基づくコードを示します。

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
from prophet import Prophet
import pandas as pd
import matplotlib.pyplot as plt

# Yahoo FinanceからAppleの株価データを取得
df = yf.download('AAPL', start='2020-01-01', end='2023-01-01')
df.reset_index(inplace=True)
df = df[['Date', 'Close']]
df.columns = ['ds', 'y']

# Prophetモデルの作成とトレーニング
model = Prophet()
model.fit(df)

# 未来のデータフレームを作成し、予測を行う
future = model.make_future_dataframe(periods=365)  # 365日先まで予測
forecast = model.predict(future)

# 予測結果の可視化
fig = model.plot(forecast)
plt.show()

このコードでは、$2020$年から$2023$年のAppleの株価データを取得し、それを基に$2024$年までの株価を予測しています。

実際のデータを用いることで、Prophetの予測性能や実用性を具体的に確認することができます。

結果解説

[実行結果]

グラフの説明をします。

このグラフはProphetを用いてAppleの株価データを予測した結果です。

以下に各要素の説明を行います。

グラフの要素

1. 黒い点（Observed Data）:

   • 黒い点は実際の観測データ（ここではAppleの過去の株価）を示しています。
   • 日付軸（$x$軸）は観測期間を示し、価格軸（$y$軸）は株価を示しています。

2. 青い線（Predicted Data）:

   • 青い線はProphetによって予測された株価を示しています。
   • 観測データの終了時点から未来のデータに対して予測が行われています。

3. 青いシェード（Uncertainty Interval）:

   • 青いシェードは予測の不確実性範囲を示しています。
   • これは$95%$信頼区間を示しており、予測値がこの範囲内にある確率が高いことを意味します。

グラフの解釈

• 過去の株価のトレンド:

  • 過去のデータ（黒い点）を見ると、株価は$2020$年から上昇トレンドを見せ、その後にいくつかのピークと谷が見られます。
  • $2021$年から$2022$年にかけて株価が特に高い値を示していますが、その後は下降トレンドに入っています。

• 未来の株価予測:

  • 青い線とシェードを見ると、$2023$年以降の株価が予測されています。
  • 予測によると、株価は$2023$年の終わりにかけてさらに変動しながらも緩やかに下がる傾向があることが示唆されています。
  • ただし、青いシェードが広がっていることから予測の不確実性が高いことも示されています。
    このため、実際の株価が予測値から大きく外れる可能性もあることを意味します。

このグラフを元に、予測の精度やモデルの調整が必要かを評価したり、他の外部要因（経済状況、企業のニュースなど）を考慮に入れて更に詳しい分析を行うことができます。

DEAPの使い方

DEAP（Distributed Evolutionary Algorithms in Python）は、進化的アルゴリズムを簡単に実装できる強力なPythonライブラリです。

DEAPは遺伝的アルゴリズム、遺伝的プログラミング、進化的戦略などの進化的計算技術をサポートします。

以下にDEAPの使い方を紹介します。

インストール

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

1

pip install deap

基本的な使い方

以下は、DEAPを使って簡単な遺伝的アルゴリズムを実装する例です。

この例では、1次元の関数$ ( f(x) = x^2 ) $の最小値を見つけることを目標とします。

ステップ1: 必要なモジュールのインポート

1
2

import random
from deap import base, creator, tools, algorithms

ステップ2: 適応度と個体の定義

1
2

creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)

ステップ3: 個体と集団の初期化関数の定義

1
2
3
4
5
6

def create_individual():
    return [random.uniform(-10, 10)]

toolbox = base.Toolbox()
toolbox.register("individual", tools.initIterate, creator.Individual, create_individual)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

ステップ4: 評価関数の定義

1
2
3
4

def evaluate(individual):
    return (individual[0]**2,)

toolbox.register("evaluate", evaluate)

ステップ5: 遺伝的操作（交叉、突然変異、選択）の定義

1
2
3

toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1)
toolbox.register("select", tools.selTournament, tournsize=3)

ステップ6: メインループの実行

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

def main():
    random.seed(42)
    
    # 初期集団の生成
    population = toolbox.population(n=300)
    
    # 統計情報の記録
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", lambda x: sum(f[0] for f in x) / len(x))
    stats.register("min", lambda x: min(f[0] for f in x))
    stats.register("max", lambda x: max(f[0] for f in x))
    
    # 遺伝的アルゴリズムの適用
    population, log = algorithms.eaSimple(population, toolbox, cxpb=0.5, mutpb=0.2, ngen=40, stats=stats, verbose=True)
    
    return population, log

if __name__ == "__main__":
    population, log = main()

詳細な説明

1. 個体と適応度の定義:

   • creator.create("FitnessMin", base.Fitness, weights=(-1.0,)) は最小化問題の適応度を定義します。
   • creator.create("Individual", list, fitness=creator.FitnessMin) は個体を定義します。

2. 個体と集団の初期化:

   • toolbox.register("individual", tools.initIterate, creator.Individual, create_individual) は個体の初期化方法を定義します。
   • toolbox.register("population", tools.initRepeat, list, toolbox.individual) は集団の初期化方法を定義します。

3. 評価関数の定義:

   • toolbox.register("evaluate", evaluate) は個体の適応度を評価する関数を登録します。

4. 遺伝的操作の定義:

   • toolbox.register("mate", tools.cxBlend, alpha=0.5) は交叉方法を定義します。
   • toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.1) は突然変異方法を定義します。
   • toolbox.register("select", tools.selTournament, tournsize=3) は選択方法を定義します。

5. メインループの実行:

   • 初期集団を生成し、遺伝的アルゴリズムを適用して進化させます。
   • algorithms.eaSimple はシンプルな進化アルゴリズムを実行します。

応用例

DEAPはこの他にも、複雑な最適化問題や機械学習モデルのハイパーパラメータチューニングなど、多岐にわたる応用が可能です。

進化的アルゴリズムの基本的な操作を組み合わせることで、さまざまな問題に対して効果的な解を見つけることができます。

結果解説

[実行結果]

gen    nevals    avg        min            max    
0      300       35.1213    0.000137411    99.7714
1      186       11.3106    0.000137411    142.071
2      189       2.06073    8.09465e-05    57.5294
3      182       0.311408    8.09465e-05    5.859  
4      201       0.0476175    1.43496e-07    0.537461
5      184       0.0264013    1.43496e-07    1.98748 
6      165       0.0187372    1.43496e-07    2.3111  
7      184       0.00584479    4.66901e-08    1.01736 
8      171       0.0143057     1.74459e-09    2.7338  
9      183       0.03766       3.42508e-11    3.13119 
10     196       0.0139445     2.70198e-13    1.71795 
11     174       0.00763203    9.98396e-12    2.2896  
12     183       0.0279807     2.19087e-12    3.16397 
13     172       0.0445861     5.71099e-14    4.10532 
14     197       0.0404401     9.19872e-16    4.16248 
15     193       0.00869942    1.0381e-16     1.68343 
16     168       0.017233      4.54004e-18    4.39632 
17     190       0.0282809     4.54004e-18    3.99426 
18     181       0.00581783    4.54004e-18    1.35411 
19     171       0.0169429     1.1955e-20     1.48987 
20     178       0.00411972    1.1955e-20     0.942579
21     168       1.49487e-05    1.1955e-20     0.00448462
22     162       0.0339901      2.64989e-24    3.10837   
23     190       0.0281804      2.64989e-24    2.458     
24     183       0.00136639     2.64989e-24    0.406396  
25     192       0.00143387     2.64989e-24    0.415259  
26     184       0.0123394      2.64989e-24    1.53045   
27     181       0.0065535      1.01205e-25    0.865262  
28     174       0.0104235      1.01205e-25    3.07322   
29     165       0.0185256      3.81818e-26    2.29195   
30     194       0.0189607      8.40103e-28    2.17665   
31     177       0.022405       5.31745e-28    4.24237   
32     168       0.00344535     1.27978e-30    0.688723  
33     198       0.0134132      1.27978e-30    1.96295   
34     196       0.0149899      2.52342e-33    2.83645   
35     174       0.0222204      8.61228e-34    3.91442   
36     199       0.0133543      2.52342e-33    2.43717   
37     186       0.0186542      1.59323e-33    2.38034   
38     185       0.0256567      9.92576e-36    4.04765   
39     182       0.0165402      9.92576e-36    3.32881   
40     171       0.00423898     8.06979e-36    0.47821   

この結果は、DEAPを用いて進化的アルゴリズムを実行した際の世代ごとの統計情報です。

各世代（gen）において、評価された個体数（nevals）、平均適応度（avg）、最小適応度（min）、最大適応度（max）が記録されています。

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

各項目の説明

1. gen: 世代数。進化のステップを示します。
2. nevals: 評価された個体の数。その世代で評価された個体の数を示します。
3. avg: 平均適応度。その世代の個体の平均適応度（この場合、xの2乗の平均）です。
4. min: 最小適応度。その世代の中で最も小さい適応度（最小値）を示します。
5. max: 最大適応度。その世代の中で最も大きい適応度（最大値）を示します。

結果の解説

初期世代（世代0）

• avg: $35.1213$
• min: $0.000137411$
• max: $99.7714$

初期集団では適応度の平均が$35.1213$で、最小値が$0.000137411$、最大値が$99.7714$です。
このことから、初期集団には非常に幅広い適応度の個体が含まれていることがわかります。

世代ごとの進化

進化が進むにつれて、適応度の平均値（avg）は減少していきます。
これは、進化が進むことで適応度が改善されている（より良い個体が選択され、交叉されている）ことを示しています。

例えば、世代$3$では、

• avg: $0.311408$
• min: $0.0000809465$
• max: $5.859$

適応度の平均が$0.311408$に減少し、最小値も$0.0000809465$と非常に小さくなっています。
進化の過程で良好な個体が選択され続けていることがわかります。

最終世代（世代$40$）

• avg: $0.00423898$
• min: $8.06979e-36$
• max: $0.47821$

最終世代では、平均適応度が$0.00423898$まで減少し、最小値は$8.06979e-36$とほぼゼロに近づいています。
これは、進化的アルゴリズムが目標とする適応度関数の最小値を非常に効果的に探索したことを示しています。

まとめ

この結果から、DEAPを用いた進化的アルゴリズムが適応度関数の最小化問題を効率的に解決できたことがわかります。

世代を重ねるごとに適応度が改善され、最終的には非常に小さな適応度の個体が得られました。

効用最大化問題 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

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

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

例題

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

ステップ

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

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

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

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は、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: 残差プロットを作成する説明変数（ここではhpとwt）。

これにより、以下のような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を目的変数、hpとwtを説明変数とする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 学校の時間割

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

問題

学校の時間割

ある学校では、次の科目を月曜日から金曜日の$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 ライブラリから Problem と AllDifferentConstraint をインポートしています。
  • 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 魔法陣

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

魔法陣とは、整数の$ ( 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 をインストールします。

1

pip install python-constraint

例題

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

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

1. 変数 x と y が $1$ から $10$ までの整数である。
2. 変数 x と y の和が $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}

このスクリプトは、x と y が $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を使用して現実的な問題を解決する方法を示します。

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

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

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

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

以下の手順で進めます：

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は、統計モデルの推定や統計検定、データの可視化などを行うための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の基本的な使用方法を示していますが、このライブラリは他にも多くの機能を提供しており、統計解析やモデリングの幅広いニーズに応えることができます。