クラスタリング

クラスタリング問題の一例として、顧客セグメンテーションを考えてみましょう。

この問題では、顧客の購買行動や属性に基づいて、顧客を異なるグループに分けることを目指します。

ここでは、年齢と年間購入額の2つの特徴を持つ顧客データを考えてみます。

まず、必要なライブラリをインポートします。

1
2
3

import matplotlib.pyplot as pltfrom
 sklearn.cluster import KMeans
import numpy as np

次に、仮想的な顧客データを作成します。

1
2
3
4
5

# 仮想的な顧客データを作成
np.random.seed(0)
age = np.random.normal(loc=40, scale=10, size=100)
annual_amount = np.random.normal(loc=50000, scale=15000, size=100)
X = np.column_stack((age, annual_amount))

このデータをKMeansクラスタリングを用いて分析します。

1
2

# KMeansクラスタリング
kmeans = KMeans(n_clusters=3, random_state=0).fit(X)

最後に、結果をプロットします。

1
2
3
4
5
6

# クラスタリング結果をプロット
plt.scatter(X[:, 0], X[:, 1], c=kmeans.labels_)
plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], s=300, c='red')
plt.xlabel('Age')
plt.ylabel('Annual Amount')
plt.show()

このコードは、年齢と年間購入額に基づいて顧客を3つのクラスタに分けます。

そして、各クラスタの中心点を赤で表示します。これにより、顧客の購買行動のパターンを視覚的に理解することができます。

[実行結果]

このコードは仮想的なデータを使用していますので、実際の問題に適用する際には、適切なデータの前処理やパラメータの調整が必要になることに注意してください。

解説

ソースコードの内容を解説します。

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

• matplotlib.pyplot: データの可視化のためのライブラリ
• sklearn.cluster.KMeans: K-meansクラスタリングアルゴリズムの実装

2. 仮想的な顧客データの作成:

• age: 年齢データを平均40、標準偏差10の正規分布からランダムに生成
• annual_amount: 年間支出データを平均50000、標準偏差15000の正規分布からランダムに生成
• X: 年齢と年間支出を結合して特徴行列を作成

3. K-meansクラスタリングの実行:

• KMeans(n_clusters=3, random_state=0): クラスタ数を3に設定し、ランダムシードを0に設定してKMeansオブジェクトを作成
• fit(X): 特徴行列Xに対してクラスタリングを実行し、各サンプルのクラスタラベルを割り当てる

4. クラスタリング結果のプロット:

• plt.scatter(X[:, 0], X[:, 1], c=kmeans.labels_): 年齢をx軸、年間支出をy軸にし、各データポイントをクラスタラベルに基づいて散布図としてプロット
• plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], s=300, c='red'): クラスタの中心点を赤い色でプロット
• plt.xlabel('Age')とplt.ylabel('Annual Amount'): x軸とy軸のラベルを設定
• plt.show(): グラフの表示

このコードは、仮想的な顧客データをK-meansクラスタリングアルゴリズムを使ってクラスタリングし、その結果を散布図として可視化します。

クラスタ数は3と設定されており、年齢と年間支出の2つの特徴量をもとに顧客をクラスタリングします。

クラスタリングの結果、各データポイントが所属するクラスタと、クラスタの中心点がプロットされます。

クロスセル予測

クロスセル予測は、顧客がある商品を購入した後に他の商品を購入する可能性を予測するための手法です。

ここでは、scikit-learnを使用して、顧客がある商品を購入した後に他の商品を購入する可能性を予測する簡単な例を示します。

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

これらのライブラリは、データの操作、可視化、モデルの訓練と評価に必要です。

1
2
3
4
5

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import confusion_matrix

2.データの生成

ここでは、500人の顧客が商品Aを購入し、そのうち200人が商品Bも購入しているというデータを生成しています。

Xは顧客の特徴量を表し、yは商品Bを購入したかどうかを表します。

1
2
3
4
5
6
7
8
9
10
11
12

# 商品Aを購入した顧客の数
n_customers = 500

# 商品Aを購入した顧客のうち、商品Bも購入した顧客の数
n_cross_sell = 200

# 商品Aを購入した顧客の特徴量（ここではランダムに生成）
X = np.random.rand(n_customers, 10)

# 商品Bを購入したかどうかのラベル（最初のn_cross_sell人は1、残りは0）
y = np.zeros(n_customers)
y[:n_cross_sell] = 1

3.データの分割

データを訓練データとテストデータに分割します。

訓練データはモデルの学習に、テストデータはモデルの性能評価に使用します。

1
2

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

4.モデルの訓練

ランダムフォレスト分類器を使用してモデルを訓練します。

ランダムフォレストは、複数の決定木を組み合わせて予測を行うアルゴリズムです。

1
2
3
4
5

# ランダムフォレスト分類器のインスタンスを作成
clf = RandomForestClassifier(n_estimators=100, random_state=42)

# モデルを訓練
clf.fit(X_train, y_train)

5.予測と評価

訓練したモデルを使用してテストデータのクロスセルを予測し、混同行列を計算します。

混同行列は、モデルの性能を評価するための表です。

1
2
3
4
5

# テストデータのクロスセルを予測
y_pred = clf.predict(X_test)

# 混同行列を計算
cm = confusion_matrix(y_test, y_pred)

6.結果の可視化

混同行列をヒートマップとして表示します。

これにより、モデルの性能を視覚的に理解することができます。

1
2
3
4
5
6
7
8
9
10
11

# 混同行列を表示
plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)
plt.title('Confusion matrix')
plt.colorbar()
tick_marks = np.arange(2)
plt.xticks(tick_marks, ['No Cross-sell', 'Cross-sell'], rotation=45)
plt.yticks(tick_marks, ['No Cross-sell', 'Cross-sell'])
plt.tight_layout()
plt.ylabel('True label')
plt.xlabel('Predicted label')
plt.show()

[実行結果]

実際の問題では、顧客の特徴量をより詳細に把握するために、購入履歴や顧客のデモグラフィック情報など、さまざまなデータを使用することがあります。

また、モデルの性能を向上させるために、パラメータのチューニングや他の機械学習アルゴリズムの使用もあります。

混同行列（Confusion Matrix）

モデルの予測結果は、混同行列（Confusion Matrix）によって表現されます。

混同行列は、実際のクロスセルのラベルと予測されたクロスセルのラベルを比較し、それらの関係を可視化します。

混同行列は4つのセルで構成されます：

• 真陽性（True Positive; TP）:
  実際のラベルが「クロスセルあり」であり、予測も「クロスセルあり」の場合。
• 真陰性（True Negative; TN）:
  実際のラベルが「クロスセルなし」であり、予測も「クロスセルなし」の場合。
• 偽陽性（False Positive; FP）:
  実際のラベルが「クロスセルなし」であり、予測が「クロスセルあり」となった場合。
• 偽陰性（False Negative; FN）:
  実際のラベルが「クロスセルあり」であり、予測が「クロスセルなし」となった場合。

混同行列の各セルには、実際のデータポイントの数が表示されます。

モデルが正しく予測した場合は真陽性（TP）や真陰性（TN）の数が高くなり、誤った予測があった場合は偽陽性（FP）や偽陰性（FN）の数が増えます。

この例では、混同行列が表示されると、x軸が予測ラベル、y軸が実際のラベルであり、それぞれのセルに該当するデータポイントの数が表示されます。

これにより、モデルの予測結果の正確さやエラーの傾向を評価することができます。

また、混同行列から算出される指標として、精度（Accuracy）、再現率（Recall）、適合率（Precision）、F1スコア（F1 Score）などがあります。

これらの指標は、モデルの性能評価や比較に使用されます。

異常検知

異常検知の一つの方法として、Isolation Forestを使用した例を以下に示します。

この例では、2次元のデータを生成し、異常値を検出します。

また、matplotlibを使用して結果をグラフ化します。

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
50
51
52

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import IsolationForest

# データ生成
rng = np.random.RandomState(42)

# 正常なデータポイント
X_train = 0.3 * rng.randn(100, 2)
X_train = np.r_[X_train + 2, X_train - 2]

# 新しい正常なデータポイント
X_test = 0.3 * rng.randn(20, 2)
X_test = np.r_[X_test + 2, X_test - 2]

# 異ント
X_outliers = rng.uniform(low=-4, high=4, size=(20, 2))

# Isolation Forestの設定
clf = IsolationForest(max_samples=100, random_state=rng)

# モデルの学習
clf.fit(X_train)

# 予測
y_pred_train = clf.predict(X_train)
y_pred_test = clf.predict(X_test)
y_pred_outliers = clf.predict(X_outliers)

# プロット
xx, yy = np.meshgrid(np.linspace(-5, 5, 50), np.linspace(-5, 5, 50))
Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)

plt.title("IsolationForest")
plt.contourf(xx, yy, Z, cmap=plt.cm.Blues_r)

b1 = plt.scatter(X_train[:, 0], X_train[:, 1], c='white',
                 s=20, edgecolor='k')
b2 = plt.scatter(X_test[:, 0], X_test[:, 1], c='green',
                 s=20, edgecolor='k')
c = plt.scatter(X_outliers[:, 0], X_outliers[:, 1], c='red',
                s=20, edgecolor='k')

plt.axis('tight')
plt.xlim((-5, 5))
plt.ylim((-5, 5))
plt.legend([b1, b2, c],
           ["training observations",
            "new regular observations", "new abnormal observations"],
           loc="upper left")
plt.show()

このコードは、Isolation Forestを使用して異常検知を行い、結果をグラフ化するものです。

グラフでは、白色の点が訓練データ、緑色の点が新たな正常データ、赤色の点が異常データを表しています。

[実行結果]

解説

このコードは、異常検知のためのアルゴリズムであるIsolation Forestを使用しています。

以下に各部分の詳細な説明を記載します。

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

必要なライブラリをインポートしています。
numpyは数値計算、matplotlibはグラフ描画、sklearn.ensembleのIsolationForestは異常検知のためのライブラリです。

2. データ生成

正常なデータポイントと異常なデータポイントを生成しています。
正常なデータは、平均0、標準偏差0.3正規分布から生成され、その後2だけシフトされます。
これを2つのクラスターとして生成しています。
異常なデータポイントは、-4から4の一様分布から生成されます。

3. Isolation Forestの設定

IsolationForestのインスタンスを作成します。
ここでは、max_samplesパラメータを100に設定し、各ツリーの最大サンプル数を100に制限しています。

4. モデルの学習

正常なデータポイントを用いてIsolationForestモデルを学習させます。

5. 予測

学習したモデルを用いて、訓練データ、新たな正常データ、異常データの各データポイントが異常かどうかを予測します。

6. プロット

最後に、matplotlibを用いて結果をプロットします。
背景の色は各点の異常スコアを表し、白色の点は訓練データ、緑色の点は新たな正常データ、赤色の点は異常データを表しています。

このコードを実行すると、Isolation Forestがどのように異常データを検出するかを視覚的に理解することができます。

風力発電の出力予測

風力発電の出力予測に関する問題を考えてみましょう。

風力発電の出力は風速に大きく依存しますが、風速は時間とともに変化するため、出力の予測は重要な課題となります。

この問題を線形回帰モデルを用いて解くことを考えます。

風速を独立変数、風力発の出力を従属変数として、最小二乗によりモデルのパラメータを最適化します。

Pythonのライブラリであるscikit-learnを用いてこの問題を解くことができます。

以下に、風力発電の出力予測のための簡単なコードを示します。

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

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

# 風速と風力発電の出力のデータ（仮）
wind_speed = np.array([2, 3, 5, 8, 10, 12, 15, 18, 20, 22, 25])
power_output = np0, 0, 10, 30, 50, 70, 90, 100, 100, 100, 100])

# 線形回帰モデルの作成と学習
model = LinearRegression()
model.fit(wind_speed.reshape(-1, 1), power_output)

# 風速と予測出力のグラフ化
plt.scatter(wind_speed, power_output, color='blue')
plt.plot(wind_speed, model.predict(wind_speed.reshape(-1, 1)), color='red')
plt.xlabel('Wind Speed')
plt.ylabel('Power Output')
plt.show()

このコードでは、風速と風力発電の出力の関係を学習し、その結果をグラフ化しています。

青色の点が実際のデータ、赤色の線が予測モデルによる出力を表しています。

ただし、実際の風力発電の出力は風速だけでなく、風向や気温、湿度などの他の要素にも影響を受けます。

また、風力発電の出力は風速が一定の値を超えると飽和するため、線形モデルよりも複雑なモデルを用いることが適切な場もあります。

[実行結果]

解説

このコードは、線形回帰を使用して風速と風力発電の出力の関係をモデリングし、その結果をグラフで表示するためのものです。

1. from sklearn.linear_model import LinearRegression:
   scikit-learnライブラリからLinearRegressionクラスをインポートします。これは、線形回帰モデルを作成するためのクラスです。

2. import matplotlib.pyplot as plt:
   matplotlibライブラリからpyplotモジュールをインポートします。
   これは、グラフを描画するための機能を提供します。

3. import numpy as np:
   numpyライブラリをインポートします。
   numpyは、数値計算や配列操作などの機能を提供します。

4. wind_speed = np.array([2, 3, 5, 8, 10, 12, 15, 18, 20, 22, 25]):
   風速のデータをnumpyの配列として定義します。

5. power_output = np.array([0, 0, 10, 30, 50, 70, 90, 100, 100, 100, 100]):
   風力発電の出力のデータをnumpyの配列として定義します。

6. model = LinearRegression():
   LinearRegressionクラスのインスタンスを作成し、modelという名前の変数に格納します。

7. model.fit(wind_speed.reshape(-1, 1), power_output):
   モデルをデータに適合させます。
   fitメソッドは、与えられた入力データ（風速）と出力データ（風力発電）を使用してモデルを学習します。

8. plt.scatter(wind_speed, power_output, color='blue'):
   scatter関数を使用して、散布図を描画します。
   x軸に風速、y軸に風力発電の出力をプロットします。

9. plt.plot(wind_speed, model.predict(wind_speed.reshape(-1, 1)), color='red'):
   plot関数を使用して、線形回帰モデルによる予測値を描画します。
   model.predictメソッドを使用して、与えられた風速に対する予測値を計算し、その結果をプロットします。

10. plt.xlabel('Wind Speed'):
    x軸のラベルを設定します。

11. plt.ylabel('Power Output'):
    y軸のラベルを設定します。

12. plt.show():
    グラフを表示します。

このコードの結果として、散布図と線形回帰モデルによる予測値を持つグラフが表示されます。

グラフを通じて、風速と風力発電の出力の関係を視覚化することができます。

散布図は青色でプロットされ、実際の風速と風力発電の出力のデータを表示します。

また、線形回帰モデルによる予測値は赤色の線で表示され、風速に対する予測された風力発電の出力を示します。

線形回帰モデルは、与えられたデータセットを基に学習し、風速と風力発電の出力の間の線形な関係を捉えようとします。

この例では、風速が増加すると風力発電の出力も増加する傾向があることが示されています。

このグラフを通じて、風速と風力発電の出力の関係性を直感的に理解することができます。

また、線形回帰モデルの予測を用いることで、未知の風速に対しても風力発電の出力を予測することが可能です。

Logitec LGB-8BNHEU3

選挙キャンペーン最適化

政治的な最適化問題として、選挙キャンペーンの最適化を考えてみましょう。

候補者は、有限のリソース（例えば、時間、お金）を使って、可能な限り多くの票を得ることを目指します。

以下に、この問題を解くための簡単なPythonコードを示します。

このコードは、各選挙区でのキャンペーン活動による投票数の増加を最大化するように、キャンペーンリソースを配分します。

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

from pulp import LpMaximize, LpProblem, LpStatus, lpSum, LpVariable
import matplotlib.pyplot as plt

# 選挙区の数
districts = ['District_A', 'District_B', 'District_C', 'District_D', 'District_E']

# 各選挙区でのキャンペーンによる投票数の増加
votes_gain = {'District_A': 300, 'District_B': 250, 'District_C': 450, 'District_D': 400, 'District_E': 350}

# 各選挙区でのキャンペーンに必要なリソース
resources_needed = {'District_A': 20, 'District_B': 25, 'District_C': 30, 'District_D': 35, 'District_E': 40}

# 利用可能なリソースの総量
total_resources = 100

# 問題の定義
model = LpProblem(name="campaign-optimization", sense=LpMaximize)

# 変数の定義
x = {i: LpVariable(name=f"x{i}", lowBound=0) for i in districts}

# 目的関数の定義
model += lpSum(votes_gain[i] * x[i] for i in districts)

# 制約条件の定義
model += lpSum(resources_needed[i] * x[i] for i in districts) <= total_resources

# 問題の解
status = model.solve()

# 結果の表示
for var in x.values():
    print(f"{var.name}: {var.value()}")

# グラフの作成
plt.bar(x.keys(), [var.value() for var in x.values()])
plt.xlabel('Districts')
plt.ylabel('Resources allocated')
plt.show()

このコードを実行すると、各選挙区に割り当てられたリソースの量が表示され、それを基にバーグラフが作成されます。

これにより、どの選挙区にどれだけのリソースが割り当てられたかを視覚的に理解できます。

[実行結果]

解説

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

1. ライブラリのインポート:
   必要なライブラリをインポートします。
   pulpは線形計画問題を解くためのライブラリで、matplotlib.pyplotはデータを視覚化するためのライブラリです。

2. データの定義:
   選挙区のリスト、各選挙区でのキャンペーンによる投票数の増加、各選挙区でのキャンペーンに必要なリソース、利用可能なリソースの総量を定義します。

3. 問題の定義:
   LpProblemクラスを使用して問題を定義します。
   ここでは、問題の名前と最大化（LpMaximize）または最小化（LpMinimize）を指定します。

4. 変数の定義:
   LpVariableクラスを使用して問題の変数を定義します。
   ここでは、各選挙区に割り当てるリソースの量を表す変数を定義しています。

5. 目的関数の定義:
   lpSum関数を使用して目的関数を定義します。
   ここでは、各選挙区でのキャンペーンによる投票数の増加を最大化することを目指しています。

6. 制約条件の定義:
   同様に、lpSum関数を使用して制約条件を定義します。
   ここでは、利用可能なリソースの総量を超えないように、各選挙区に割り当てるリソースの量を制限しています。

7. 問題の解:
   solveメソッドを使用して問題を解きます。

8. 結果の表示:
   最後に、各選挙区に割り当てられたリソースの量を表示します。

9. グラフの作成:
   matplotlib.pyplotを使用して、各選挙区に割り当てられたリソースの量をバーグラフで表示します。

このコードは、選挙キャンペーンのリソース配分を最適化する一例です。
具体的な数値や制約は、実際の状況に応じて調整することができます。

無線通信ネットワーク最適化

無線通信ネットワーク最適化の一例として、ネットワークの帯域幅割り当て問題を考えてみましょう。

この問題では、各ユーザーが必要なデータレートを満たすために最小の帯域幅を割り当てることを目指します。

以下に、この問題を解くためのPythonとPuLPを使用したコードを示します。

この例では、3人のユーザーがいると仮定します。

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

from pulp import *
import numpy as np
import matplotlib.pyplot as plt

# ユーザー数
n_users = 3

# ユーザーごとの必要なデータレート (Mbps)
data_rates = np.array([5, 10, 8])

# ユーザーごとの信号対雑音比 (SNR) (dB)
snrs = np.array([10, 20, 15])

# 帯域幅 (Hz)
B = 2000000

# ノイズパワー (W)
N = 10**-13

# 問題の定義
prob = LpProblem("Bandwidth_Allocation", LpMinimize)

# 変数の定義
x = [LpVariable(f"x{i}", lowBound=0) for i in range(n_users)]

# 目的関数
prob += lpSum(x)

# 制約条件
for i in range(n_users):
    prob += x[i] * np.log2(1 + snrs[i]) >= data_rates[i]

# 最適化問題を解く
prob.solve()

# 結果を表示
print("最適な帯域幅割り当て: ", [value(xi) for xi in x])

# グラフ化
plt.bar(range(n_users), [value(xi) for xi in x])
plt.xlabel('User')
plt.ylabel('Bandwidth')
plt.show()

このコードは、各ユーザーに必要な最小の帯域幅を割り当てる問題を解きます。

目的関数は、割り当てられた帯域幅の合計を最小化することです。

制約条件は、各ユーザーが必要なデータレートを満たすことを保証します。

最後に、各ユーザーに割り当てられた最適な帯域幅を棒グラフで表示します。

これにより、各ユーザーがどの程度の帯域幅を使用しているかを視覚的に理解することができます。

[実行結果]

解説

この問題は、無線通信ネットワークにおける帯域幅割り当ての問題を解くものです。

各ユーザーは一定のデータレートを必要としており、それを満たすためには一定の帯域幅が必要です。

しかし、全体の帯域幅は限られているため、どのユーザーにどれだけの帯域幅を割り当てるかを最適化する必要があります。

目的関数は、割り当てられた帯域幅の合計を最小化することです。

これは、全体の帯域幅を最も効率的に使用するためのものです。

制約条件は、各ユーザーが必要なデータレートを満たすことを保証します。

これは、ユーザーごとの信号対雑音比（SNR）と割り当てられた帯域幅を用いて計算されます。

具体的には、シャノンの定理を用いて、割り当てられた帯域幅とSNRから達成可能なデータレートを計算し、それがユーザーの必要なデータレート以上であることを保証します。

この問題を解くと、各ユーザーに割り当てるべき最適な帯域幅が得られます。

そして、その結果を棒グラフで表示することで、各ユーザーがどの程度の帯域幅を使用しているかを視覚的に理解することができます。

ソフトマージンサポートベクターマシン

ソフトマージンサポートベクターマシン(SVM)は、データが完全に分離できない場合に使用されます。

以下に、2次元のデータセットに対するソフトマージンSVMの例を示します。

この例では、CVXPYを使用して最適化問題を解き、matplotlibを使用して結果をグラフ化します。

まず、必要なライブラリをインポートします。

1
2
3

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

次に、ランダムな2次元データセットを生成します。

1
2
3
4
5
6

np.random.seed(1)
n = 50
d = 2
X = np.random.randn(n, d)
delta = np.random.uniform(-0.5, 0.5, size=(n,1))
Y = 2*(X[:,0:1] + delta > X[:,1:2]) - 1

次に、ソフトマージンSVMの最適化問題を定義します。

1
2
3
4
5
6
7
8
9

w = cp.Variable((d, 1))
b = cp.Variable()
xi = cp.Variable((n, 1))
C = 1.0

objective = cp.Minimize(cp.norm(w, 2) + C*cp.sum(xi))
constraints = [cp.multiply(Y, X @ w + b) >= 1 - xi, xi >= 0]
prob = cp.Problem(objective, constraints)
prob.solve()

最後に、結果をグラフ化します。

1
2
3
4
5
6
7
8
9
10

plt.figure(figsize=(8,8))
plt.scatter(X[:,0], X[:,1], c=Y.flatten())
slope = -w.value[0] / w.value[1]
intercept = -b.value / w.value[1]
x_vals = np.linspace(-2, 2, 100)
plt.plot(x_vals, slope * x_vals + intercept, 'k-')
plt.fill_between(x_vals, slope * x_vals + intercept - 1 / np.linalg.norm(w.value), slope * x_vals + intercept + 1 / np.linalg.norm(w.value), color='k', alpha=0.1)
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.show()

このコードは、2次元空間にランダムに生成されたデータポイントをプロットし、ソフトマージンSVMによって計算された分離超平面（黒線）とマージン（灰色の領域）を表示します。

[実行結果]

解説

この例では、ソフトマージンサポートベクターマシン（SVM）を使用して、2次元空間にランダムに生成されたデータポイントを分類しています。

この超平面は、異なるクラスのデータポイント間のマージン（つまり、最も近いデータポイントまでの距離）を最大化します。

これは、スラック変数ξ（xi）を導入することで実現され、これにより一部のデータポイントがマージンの「間違った側」に存在することを許容します。

この例の結果をグラフ化すると、黒線がSVMによって計算された分離超平面を表し、灰色の領域がマージンを示します。

データポイントは色でクラスが示され、一部のデータポイントがマージンを侵害していることがわかります。

これはソフトマージンSVMが許容する「間違い」を示しています。

最適化問題の目的関数は、マージンを最大化（つまり、wのノルムを最小化）し、同時にマージン侵害を最小化（つまり、スラック変数ξの合計を最小化）することです。

ここで、Cはこれら二つの目標の間のトレードオフを制御するパラメータです。

Cが大きいほど、マージン侵害をより厳しく罰することになります。

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

CVXPYを使用してポートフォリオ最適化問題を解く基本的なPythonコードを示します。

この例では、3つの異なる資産に対する投資の最適な割合を見つけることを目指しています。

まず、必要なライブラリをインポートします。

1
2
3

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

次に、各資産のリターンとリスク（ここでは標準偏差）を表すデータを作成します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 資産のリターン
returns = np.array([0.1, 0.2, 0.15])

# 資産のリスク（標準偏差）
risks = np.array([0.05, 0.1, 0.08])

# 資産間の相関係数
correlation = np.array([
    [1, 0.2, 0.5],
    [0.2, 1, 0.3],
    [0.5, 0.3, 1]
])

# 共分散行列を計算
cov_matrix = np.outer(risks, risks) * correlation

次に、CVXPYを使用して最適化問題を定義します。

ここでは、リターンを最大化し、リスク（ここではポートフォリオの標準偏差）を制約としています。

1
2
3
4
5
6
7
8
9
10
11
12
13

# 投資の割合を表す変数
weights = cp.Variable(len(returns))

# リターンを最大化する問題を定義
objective = cp.Maximize(weights @ returns)

# 制約条件を定義
constraints = [cp.sum(weights) == 1,
               weights >= 0,
               cp.quad_form(weights, cov_matrix) <= 0.05**2]

# 問題を定義
problem = cp.Problem(objective, constraints)

最後に、問題を解き、最適な投資の割合を表示します。

1
2
3
4
5

# 問題を解く
problem.solve()

# 最適な投資の割合を表示
print(weights.value)

このコードは、与えられたリターンとリスクに基づいて、3つの資産に対する最適な投資の割合を計算します。

制約条件は、全ての投資の割合の合計が1（つまり、全ての資産に対する投資の割合）であることと、ポートフォリオのリスクが特定の値以下であることを保証します。

最後に、各資産の投資の割合をグラフ化します。

1
2
3
4

plt.bar(range(len(weights.value)), weights.value)
plt.xlabel('Assets')
plt.ylabel('Investment proportion')
plt.show()

このグラフは、各資産に対する最適な投資の割合を視覚的に示しています。

[実行結果]

栄養最適化

食事に関する最適化問題の一つとして、特定の栄養素を必要量摂取しつつ、食事のコストを最小化する問題があります。

以下に、この問題をCVXPYを用いて解く例を示します。

まず、必要なライブラリをインポートします。

1
2
3

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

次に、問題のパラメータを定義します。ここでは、3つの食品（リンゴ、バナナ、オレンジ）と2つの栄養素（ビタミンCとカリウム）を考えます。

1
2
3
4
5
6
7
8

# 各食品のコスト（円）
costs = np.array([100, 50, 80])

# 各食品のビタミンCとカリウムの含有量
nutrients = np.array([[0.5, 0.2], [0.1, 0.3], [0.2, 0.2]])

# 必要なビタミンCとカリウムの量
required_nutrients = np.array([1.0, 0.5])

CVXPYを用いて最適化問題を定義します。

1
2
3
4
5
6
7

# 食品の購入量を表す変数
x = cp.Variable(3, nonneg=True)

# コスト最小化問題を定義
objective = cp.Minimize(costs @ x)
constraints = [nutrients.T @ x >= required_nutrients, x >= 0]
problem = cp.Problem(objective, constraints)

最適化問題を解きます。

1

problem.solve()

最適な購入量を表示します。

1
2
3

print(f"Optimal purchase quantity for Apple: {x.value[0]}")
print(f"Optimal purchase quantity for Banana: {x.value[1]}")
print(f"Optimal purchase quantity for Orange: {x.value[2]}")

最後に、この問題とその解をグラフ化します。

ただし、3つ以上の食品を扱う場合、それら全てを一つのグラフに描画することは難しいため、各食品の最適な購入量を棒グラフで表示します。

1
2
3
4

plt.bar(['Apple', 'Banana', 'Orange'], x.value)
plt.xlabel('Food')
plt.ylabel('Optimal purchase quantity')
plt.show()

このコードは、リンゴ、バナナ、オレンジの最適な購入量を示す棒グラフをプロットします。

[実行結果]

値引き最適化

値引き最適化問題とは、特定の制約条件下で最大の利益を得るための最適な値引き率を求める問題です。

ここでは、ある店が3つの商品を販売していて、それぞれの商品に対する最大の値引き率と、それぞれの商品がもたらす利益が異なるという状況を考えてみましょう。

まずは、必要なライブラリをインポートします。

1
2
3

from pulp import *
import matplotlib.pyplot as plt
import numpy as np

次に、問題を定義します。

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

# 問題の定義
prob = LpProblem("Discount Optimization", LpMaximize)

# 変数の定義
discounts = ['product_1', 'product_2', 'product_3']
discount_vars = LpVariable.dicts("Discount", discounts, 0, 1)

# 利益の定義
profits = {'product_1': 100, 'product_2': 200, 'product_3': 300}

# 目的関数の定義
prob += lpSum([profits[i]*discount_vars[i] for i in discounts])

# 制約条件の定義
prob += lpSum([discount_vars[i] for i in discounts]) <= 1.5, "Total Discount Constraint"

# 問題の解決
prob.solve()

最後に、結果を表示し、グラフ化します。

1
2
3
4
5
6
7
8
9
10
11
12
13

# 結果の表示
for v in prob.variables():
    print(v.name, "=", v.varValue)

# グラフ化
plt.figure(figsize=(8,6))
plt.bar(range(len(discounts)), [value(var) for var in prob.variables()])

plt.xticks(range(len(discounts)), [var.name for var in prob.variables()])
plt.xlabel('Products')
plt.ylabel('Discount Rate')
plt.title('Optimal Discount Rates for Products')
plt.show()

このコードは、3つの商品に対する最適な値引き率を求め、それを棒グラフで表示します。

制約条件は、全商品の値引き率の合計が1.5以下であることを指定しています。

また、目的関数は各商品の利益と値引き率の積の合計を最大化することを目指しています。

これにより、最大の利益を得るための最適な値引き率が求められます。

[実行結果]