グループ化／集約

Polarsの一つの強みは、大規模データセットを高速に処理する能力です。

その最適化された処理エンジンのおかげで、大量のデータを即座に変換、整理したり、集計することが容易となり、パフォーマンスを大いに向上させます。

また、シンプルで明確なAPIもPolarsの大きな魅力の一つです。

これにより、ユーザーはpandasのような一般的なデータ分析ライブラリと同様のユーザビリティを享受しながら、より大きなデータ量を扱うことができます。

以下に、Polarsを使ったデータのグループ化と集約のサンプルコードを示します。

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

import polars as pl

# データフレームの作成
df = pl.DataFrame({
    "fruits": ["apple", "banana", "banana", "apple", "apple", "banana", "banana", "apple"],
    "sales": [5, 6, 8, 9, 5, 6, 7, 9],
    "city": ["Tokyo", "Tokyo", "Tokyo", "London", "London", "London", "New York", "New York"]
})

# 都市ごと、果物ごとの総売上を計算
df_grouped = df.groupby(["city", "fruits"]).agg([pl.col("sales").sum().alias("total_sales")])

print(df_grouped)

上記のコードは、都市ごとに果物の売上合計を計算するシンプルな例です。

都市と果物ごとにグループ化し、各グループの“sales”列の合計を計算します。

[実行結果]

shape: (6, 3)
┌──────────┬────────┬─────────────┐
│ city     ┆ fruits ┆ total_sales │
│ ---      ┆ ---    ┆ ---         │
│ str      ┆ str    ┆ i64         │
╞══════════╪════════╪═════════════╡
│ Tokyo    ┆ apple  ┆ 5           │
│ Tokyo    ┆ banana ┆ 14          │
│ New York ┆ banana ┆ 7           │
│ London   ┆ apple  ┆ 14          │
│ London   ┆ banana ┆ 6           │
│ New York ┆ apple  ┆ 9           │
└──────────┴────────┴─────────────┘

shape: (6, 3) はこのデータフレームが6行3列であることを示しています。
列は city、fruits、そして total_sales となります。

それぞれの行を見てみましょう：

• 1行目：
  東京でのリンゴの総売上は5と計算されています。
• 2行目：
  東京でのバナナの総売上は14と計算されています。
• 3行目：
  ニューヨークでのバナナの総売上は7と計算されています。
• 4行目：
  ロンドンでのリンゴの総売上は14と計算されています。
• 5行目：
  ロンドンでのバナナの総売上は6と計算されています。
• 6行目：
  ニューヨークでのリンゴの総売上は9と計算されています。

これにより、各都市で各フルーツの総販売量を直観的に理解することができます。

このようなデータの表示は、日々の販売活動の理解や、将来的な販売戦略の策定に非常に有用です。

フルーツの種類と都市に基づいてデータを分析することで、どの商品がどの市場でよく売れているのか洞察を得ることが可能です。

Polarsの概要

Polarsは、高速なDataFrame操作を可能にするライブラリで、Apache Arrowの列指向メモリモデルを使用しています。

Polarsは、Rustで書かれており、多スレッドのクエリエンジンを使用して効率的な並列処理を実現しています。

また、Polarsはオープンソースであり、MITライセンスで利用できます。

サンプルソース

以下に、Polarsを使用してDataFrameを作成し、そのDataFrameに対して操作を行う例を示します。

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

1

import polars as pl

次に、DataFrameを作成します。

1
2
3
4
5
6
7

df = pl.DataFrame(
 {
     "name": ["John", "Sue", "Tom", "Alice"],
     "age": [23, 21, 19, 18],
     "grade": [90, 85, 78, 92],
 }
)

次に、DataFrameに対して操作を行います。

例えば、特定の列を選択したり、特定の条件に基づいて行をフィルタリングしたり、新しい列を作成したりすることができます。

1
2
3
4
5
6
7
8

# 特定の列を選択
selected_df = df.select(["name", "age"])

# 特定の条件に基づいて行をフィルタリング
filtered_df = df.filter(pl.col("age") > 20)

# 新しい列を作成
df = df.with_columns(pl.col("age") + 1)

最後に、DataFrameを表示します。

1
2
3

print(selected_df)
print(filtered_df)
print(df)

[実行結果]

shape: (4, 2)
┌───────┬─────┐
│ name  ┆ age │
│ ---   ┆ --- │
│ str   ┆ i64 │
╞═══════╪═════╡
│ John  ┆ 23  │
│ Sue   ┆ 21  │
│ Tom   ┆ 19  │
│ Alice ┆ 18  │
└───────┴─────┘
shape: (2, 3)
┌──────┬─────┬───────┐
│ name ┆ age ┆ grade │
│ ---  ┆ --- ┆ ---   │
│ str  ┆ i64 ┆ i64   │
╞══════╪═════╪═══════╡
│ John ┆ 23  ┆ 90    │
│ Sue  ┆ 21  ┆ 85    │
└──────┴─────┴───────┘
shape: (4, 3)
┌───────┬─────┬───────┐
│ name  ┆ age ┆ grade │
│ ---   ┆ --- ┆ ---   │
│ str   ┆ i64 ┆ i64   │
╞═══════╪═════╪═══════╡
│ John  ┆ 24  ┆ 90    │
│ Sue   ┆ 22  ┆ 85    │
│ Tom   ┆ 20  ┆ 78    │
│ Alice ┆ 19  ┆ 92    │
└───────┴─────┴───────┘

このコードは、Polarsを使用してDataFrameを作成し、そのDataFrameに対して操作を行い、その結果を表示します。

このように、Polarsを使用することで、DataFrameの操作を効率的に行うことができます。

Polarsは、高速なDataFrame操作を可能にするライブラリで、多くの操作が可能です。

以下に、一部の操作を示します。

1. 行の選択：

DataFrame内の1行を選択するには、row()メソッドを使用して行番号を渡します。
複数の行を選択するためには、filter()関数を使用します。

1
2

# 最初の行を取得する
df.row(0)

[実行結果]

('John', 24, 90)

1
2
3

# 'name'列の値が'Tom'の行を取得する
res = df.filter(pl.col('name') == 'Tom')
print(res)

[実行結果]

shape: (1, 3)
┌──────┬─────┬───────┐
│ name ┆ age ┆ grade │
│ ---  ┆ --- ┆ ---   │
│ str  ┆ i64 ┆ i64   │
╞══════╪═════╪═══════╡
│ Tom  ┆ 20  ┆ 78    │
└──────┴─────┴───────┘

2. 論理演算子の使用：

Polarsでは、論理演算子を使用して複数の条件を指定できます。

1
2
3

# 'name'列の値が'John'または'Alice'の行を取得する
res = df.filter((pl.col('name') == 'John') | (pl.col('name') == 'Alice'))
print(res)

[実行結果]

shape: (2, 3)
┌───────┬─────┬───────┐
│ name  ┆ age ┆ grade │
│ ---   ┆ --- ┆ ---   │
│ str   ┆ i64 ┆ i64   │
╞═══════╪═════╪═══════╡
│ John  ┆ 24  ┆ 90    │
│ Alice ┆ 19  ┆ 92    │
└───────┴─────┴───────┘

3. 結果の格納：

Polarsは演算結果を新しい列として格納する際の記述方法に優れています。
例えば、新しい列を作成するときには、普通の四則演算以外の結果も格納できます。

1
2
3

# 'age'列の値を2で割った結果を新しい列'age_half'に格納する
df2 = df.with_columns(pl.col("age") / 2)
print(df2)

[実行結果]

shape: (4, 3)
┌───────┬──────┬───────┐
│ name  ┆ age  ┆ grade │
│ ---   ┆ ---  ┆ ---   │
│ str   ┆ f64  ┆ i64   │
╞═══════╪══════╪═══════╡
│ John  ┆ 12.0 ┆ 90    │
│ Sue   ┆ 11.0 ┆ 85    │
│ Tom   ┆ 10.0 ┆ 78    │
│ Alice ┆ 9.5  ┆ 92    │
└───────┴──────┴───────┘

これらの操作は、Polarsの強力な機能の一部を示しています。

Polarsは、DataFrameの操作を効率的に行うための強力なツールであり、これらの操作を通じて、データ分析やデータ操作の幅広いタスクを効率的に実行することができます。

色彩の相補色

美術関連の問題として、色彩の相補色を見つける問題を考えてみましょう。

相補色は、RGB色空間内で最も対照的な色です。

与えられた色に対して、それに最も近い相補色を見つけることができます。

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import distance

# 関数定義: RGB色空間内で最も近い相補色を計算する
def find_complementary_color(color):
    # 相補色の定義: RGB値の各要素を1から引くことで相補色を得る
    complementary_color = [1 - channel for channel in color]
    return complementary_color

# 与えられた色（RGB値）
given_color = [0.4, 0.2, 0.9]  # 例えば青系統の色

# 相補色を計算
complementary_color = find_complementary_color(given_color)

# 色の距離を計算
distance_value = distance.euclidean(given_color, complementary_color)

# グラフ化
colors = [np.array(given_color), np.array(complementary_color)]
labels = ['Given Color', 'Complementary Color']

plt.figure(figsize=(8, 4))
for i in range(2):
    plt.subplot(1, 2, i+1)
    plt.bar(range(len(colors[i])), colors[i], color=['r', 'g', 'b'], alpha=0.6)
    plt.xticks(range(len(colors[i])), ['R', 'G', 'B'])
    plt.title(labels[i])

plt.tight_layout()
plt.show()

print(f"The given color is: {given_color}")
print(f"The complementary color is: {complementary_color}")
print(f"The distance between these colors is: {distance_value}")

このコードは、与えられた色（例えば青系統の色）のRGB値を取り、その色に最も近い相補色を見つけます。

そして、2つの色を棒グラフで示し、それらの色の距離を表示します。

色の距離が大きいほど、その色はより対照的であると考えられます。

[実行結果]

ソースコード解説

このソースコードは、Pythonで色彩処理を行うためのプログラムです。

以下はコードの構造と機能の詳細な説明です：

ライブラリのインポート

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import distance

• numpy：数値計算を行うためのライブラリ。
• matplotlib.pyplot：グラフを描画するためのライブラリ。
• scipy.spatial.distance：距離を計算するためのライブラリ。

相補色を見つける関数の定義

1
2
3

def find_complementary_color(color):
    complementary_color = [1 - channel for channel in color]
    return complementary_color

• find_complementary_color 関数は、与えられた色に対する相補色を計算します。
• 各RGBチャンネルの値から1を引くことで相補色を生成しています。

与えられた色の定義

1

given_color = [0.4, 0.2, 0.9]

• この行では、青系統の色を表すRGB値 [0.4, 0.2, 0.9] を定義しています。

相補色の計算

1

complementary_color = find_complementary_color(given_color)

• find_complementary_color 関数を使って、与えられた色の相補色を計算しています。

色の距離の計算

1

distance_value = distance.euclidean(given_color, complementary_color)

• scipy の distance モジュールの euclidean 関数を使って、与えられた色とその相補色の間の距離を計算しています。

グラフの描画

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

colors = [np.array(given_color), np.array(complementary_color)]
labels = ['Given Color', 'Complementary Color']

plt.figure(figsize=(8, 4))
for i in range(2):
    plt.subplot(1, 2, i+1)
    plt.bar(range(len(colors[i])), colors[i], color=['r', 'g', 'b'], alpha=0.6)
    plt.xticks(range(len(colors[i])), ['R', 'G', 'B'])
    plt.title(labels[i])

plt.tight_layout()
plt.show()

• 2つのグラフを描画します。左側には与えられた色、右側にはその相補色を示すバー形式のグラフが表示されます。
• 色の強度は、RGBチャンネルの値を表しており、各バーの高さが対応する色の成分を示しています。

結果の出力

1
2
3

print(f"The given color is: {given_color}")
print(f"The complementary color is: {complementary_color}")
print(f"The distance between these colors is: {distance_value}")

• 最後に、与えられた色、その相補色、およびその間の距離をコンソールに表示します。

このプログラムは、与えられた色とその相補色の関係性を理解しやすく視覚化し、色の距離を数値化して表示することで、色彩理論に関する概念を理解するのに役立ちます。

結果解説

下記の結果は、与えられた色とその相補色に関する情報を示しています。

[実行結果]

The given color is: [0.4, 0.2, 0.9]
The complementary color is: [0.6, 0.8, 0.09999999999999998]
The distance between these colors is: 1.019803902718557

1. 与えられた色:

• RGB値が [0.4, 0.2, 0.9] となっています。
• これは、赤の成分が40%、緑の成分が20%、青の成分が90%である色を示しています。
  青が最も強く、赤と緑が弱い色です。

2. 相補色:

• 計算された相補色のRGB値は [0.6, 0.8, 0.09999999999999998] です。
• この色は、与えられた色と対照的であり、赤と緑の成分が高く、青の成分が低い色です。
  つまり、青系統の色に対して赤と緑が強く現れる色です。

3. 色の距離:

• 与えられた色と相補色の間の距離は、1.019803902718557 です。
• この値は、色彩空間内で2つの色がどれだけ対照的かを示しています。
  この場合、与えられた色とその相補色はかなり対照的であることを示しています。

グラフ解説

[実行結果]

このグラフは、与えられた色とその相補色を示しています。

それぞれの棒グラフは、RGB色空間内の3つのチャンネル（赤、緑、青）を表しています。

1. 左側のグラフ（Given Color）：

   • 与えられた色のRGB値が示されています。
   • 各バーは、赤、緑、青のチャンネルを示し、それぞれの色の強度を表しています。
   • 例えば、青系統の色が与えられた場合、青のバーが最も高くなるでしょう。

2. 右側のグラフ（Complementary Color）：

   • 計算された相補色のRGB値が示されています。
   • 同様に、赤、緑、青のバーがそれぞれの色の強度を表しています。
   • 相補色は、与えられた色と対照的な色であるため、バーの高さが与えられた色とは異なる場合があります。

これらのグラフを見ることで、与えられた色とその相補色の各色成分の相違や類似性が視覚的に理解できます。

また、それらの色の距離も計算され、その差異の大きさが数値として示されています。

色の距離が大きいほど、色の相補性が高いと考えられます。

振動子の運動方程式

物理に関連する問題として、振動子の運動方程式を解いてみましょう。

例として、単振動子の運動を考えます。

単振動子は、バネにつながれた質点が減衰力や外力を受けずに振動するシンプルなモデルです。

振動子の運動方程式は、次のように表されます。

$$
[
m \frac{d^2x}{dt^2} + c \frac{dx}{dt} + kx = 0
]
$$

ここで、$ (m) $は質量、$ (c) $は減衰係数、$ (k) $はバネ定数です。

これをSciPyを使って数値的に解いて、結果をグラフ化してみましょう。

以下は、質量$ (m = 1) $、減衰係数$ (c = 0.1) $、バネ定数$ (k = 1) $の単振動子の運動方程式を解く例です。

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

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# 運動方程式を定義
def harmonic_oscillator(x, t, m, c, k):
    dxdt = x[1]  # 速度を表す状態変数
    dvdt = -(c/m) * x[1] - (k/m) * x[0]  # 加速度を表す状態変数
    return [dxdt, dvdt]

# 初期条件とパラメータ
x0 = [1.0, 0.0]  # 初期位置と初速度
m = 1.0          # 質量
c = 0.1          # 減衰係数
k = 1.0          # バネ定数

# 時間の範囲
t = np.linspace(0, 20, 1000)

# 運動方程式を解く
solution = odeint(harmonic_oscillator, x0, t, args=(m, c, k))

# 位置と速度の時間変化をプロット
plt.figure(figsize=(8, 6))
plt.plot(t, solution[:, 0], label='Position')
plt.plot(t, solution[:, 1], label='Velocity')
plt.xlabel('Time')
plt.ylabel('Value')
plt.title('Harmonic Oscillator Motion')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、Scipyの odeint を使用して単振動子の運動方程式を数値的に解いています。

結果は時間と共に変化する位置と速度の変化をグラフ化しています。

位置と速度のグラフが振動する様子が見られます。

[実行結果]

ソースコード解説

このコードは、単振動子（harmonic oscillator）の運動方程式を解いて、その結果をグラフ化するものです。

以下にソースコードの各部分を説明します。

1. モジュールのインポート:

1
2
3

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

• numpy は数値計算用ライブラリで、数値配列を処理するために使われます。
• scipy.integrate の odeint は常微分方程式を数値的に解くための関数です。
• matplotlib.pyplot はグラフ描画ライブラリで、データを可視化するために使用されます。

2. 運動方程式の定義:

1
2
3
4

def harmonic_oscillator(x, t, m, c, k):
   dxdt = x[1]  # 速度を表す状態変数
   dvdt = -(c/m) * x[1] - (k/m) * x[0]  # 加速度を表す状態変数
   return [dxdt, dvdt]

• harmonic_oscillator 関数は、単振動子の運動方程式を定義しています。
• 位置$ (x) $と時間$ (t) $を引数として受け取り、速度$ (dx/dt) $と加速度$ (d^2x/dt^2) $を計算して返します。

3. 初期条件とパラメータの設定:

1
2
3
4

x0 = [1.0, 0.0]  # 初期位置と初速度
m = 1.0          # 質量
c = 0.1          # 減衰係数
k = 1.0          # バネ定数

• 初期位置と初速度、質量、減衰係数、バネ定数を設定します。

4. 時間の設定:

1

t = np.linspace(0, 20, 1000)

• np.linspace を使って、$0$から$20$までの範囲を$1000$個の等間隔な点に区切り、時間の配列を作成します。

5. 運動方程式の解を計算:

1

solution = odeint(harmonic_oscillator, x0, t, args=(m, c, k))

• odeint 関数を使用して、運動方程式を数値的に解きます。
• harmonic_oscillator 関数に初期条件とパラメータを渡し、時間 t の範囲での解を求めます。

6. 位置と速度の時間変化をプロット:

1
2
3
4
5
6
7
8
9

plt.figure(figsize=(8, 6))
plt.plot(t, solution[:, 0], label='Position')
plt.plot(t, solution[:, 1], label='Velocity')
plt.xlabel('Time')
plt.ylabel('Value')
plt.title('Harmonic Oscillator Motion')
plt.legend()
plt.grid(True)
plt.show()

• matplotlib を使用して、時間に対する位置と速度の変化を折れ線グラフでプロットします。
• solution[:, 0] は位置、solution[:, 1] は速度の結果を表します。
• グラフには軸ラベル、タイトル、凡例、グリッド線が追加されます。
• 最後に plt.show() でグラフを表示します。

グラフ解説

[実行結果]

このグラフは、単振動子の振動を時間と共に表現したものです。

1. 位置（Position）のグラフ:

• 横軸は時間を表し、縦軸は単振動子の位置を示しています。
• グラフは、時間の経過とともに振動子の位置が変化する様子を示しています。
  振動子は、特定の周期で正弦波のような動きをします。
  このグラフはその振動の様子を可視化しています。

2. 速度（Velocity）のグラフ:

• 横軸は時間を表し、縦軸は単振動子の速度を示しています。
• 速度グラフは、時間と共に速度がどのように変化するかを示しています。
  振動子は、位置が最大または最小になるときに速度が0になり、その間に速度が最大になることが一般的です。
  このグラフはその速度の変化を示しています。

これらのグラフから、単振動子が周期的に位置と速度が変化する様子を観察できます。

位置と速度の変化が互いに90度ずれた位相で起こり、振動が継続していることがわかります。

航空旅客数の予測

Facebookが開発した時系列分析ライブラリProphetを用いて、航空旅客数のデータを分析し、将来の旅客数を予測してみましょう。

まず、以下のようにライブラリをインストールします。

1

!pip install prophet

次に、以下のコードを使ってデータを分析・予測します。

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

import pandas as pd
from prophet import Prophet
import matplotlib.pyplot as plt

# データの読み込み
df = pd.read_csv('https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv')

# 列名の変更
df.columns = ['ds', 'y']
df['ds'] = pd.to_datetime(df['ds'])

# モデルの訓練
m = Prophet(yearly_seasonality=True)
m.fit(df)

# 未来の予測のためのデータフレームの作成
future = m.make_future_dataframe(periods=365)
forecast = m.predict(future)

# 予測のプロット
m.plot(forecast)
plt.show()

# トレンドとコンポーネントの表示
m.plot_components(forecast)
plt.show()

ここでは、航空旅客数のデータを読み込み、Prophetのモデルに適用可能な形式に加工します。

その後、モデルを訓練し、将来の365日にわたる予測を行います。

m.plot(forecast)では、元のデータ（点）と予測の結果（線）をプロットしています。

また、m.plot_components(forecast)では、予測の各コンポーネント（トレンド、年間の季節性）をプロットしています。

これらの結果から、旅客数の長期的なトレンドや季節的な変動を理解することができ、さらに未来の旅客数を予測することが可能です。

[実行結果]

ソースコード解説

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

1. データの読み込み

• pd.read_csv を使用して、航空旅客数の時系列データをURLから読み込みます。
• データフレームの列名を['ds', 'y'] に変更し、'ds' 列を日付型に変換します。

2. Prophetモデルの訓練

• Prophetを初期化し、年次の季節性を有効にしてモデルを構築します。
• 訓練用のデータを使って、m.fit(df) でモデルを訓練します。

3. 未来の予測のためのデータフレームの作成

• 未来の予測のために、m.make_future_dataframe(periods=365) を使って新しい日付のデータフレームを作成します。
• この新しいデータフレームに対して m.predict(future) を実行し、将来の予測を行います。

4. 予測のプロット

• m.plot(forecast) を使って、予測されたデータを含むグラフを表示します。
• 実際のデータと予測されたデータが可視化され、さらに予測の信頼区間も表示されます。

5. トレンドとコンポーネントの表示

• m.plot_components(forecast) を使用して、Prophetが検出したトレンドや季節性の成分を個別にプロットします。
• トレンド、年次、週次、日次の季節性成分が表示され、データのパターンや季節変動を視覚的に理解するのに役立ちます。

これらのステップにより、Prophetを使って時系列データの予測を行い、その結果をグラフで視覚化する完全な流れが示されています。

グラフ解説

このコードは、Prophetを使って時系列データの予測とその分解を行います。

データセットは航空旅客数の時系列データで、以下の3つのグラフを生成します。

[実行結果]

1. 全体の時系列データの予測

   • Prophetによる時系列データの予測を可視化します。
     横軸は時間（日付）を、縦軸は航空旅客数を表します。
     黒色の点が実際のデータ、青色の線が予測されたデータを示します。
     さらに、予測の信頼区間も表示されます（水色の領域）。

2. トレンド

   • Prophetが検出したトレンドを示します。
     時系列データの中での長期的な変動を捉えており、時間に伴う変動や推移を示します。
     年次（yearly）、週次（weekly）、日次（daily）の季節性成分も含まれます。

3. 季節性成分の表示

   • 季節性成分を個別にプロットし、それぞれの季節パターンを示します。
     年間の季節変動や週ごとのパターン、さらに日単位の変動を示しています。

これらのグラフは、時系列データの予測に使用されるProphetモデルの挙動や、データが持つトレンドや季節性などを視覚的に理解するのに役立ちます。

心臓疾患の予測

scikit-learnを使用して、医学的な問題の一つである心臓疾患の予測を行う例を示します。

この問題では、心臓疾患の予測を行うための機械学習モデルを作成し、その結果をグラフ化します。

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

1
2
3
4
5
6

from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt
import seaborn as sns

次に、心臓疾患のデータセットをロードします。

scikit-learnのdatasetsモジュールには、心臓疾患のデータセットが含まれています。

1
2
3

heart = datasets.load_breast_cancer()
X = heart.data
y = heart.target

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

1

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

次に、RandomForestClassifierを使用してモデルを訓練します。

1
2

clf = RandomForestClassifier(random_state=42)
clf.fit(X_train, y_train)

次に、モデルの予測精度を計算します。

1
2
3

y_pred = clf.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)
print(f'Accuracy: {accuracy}')

最後に、特徴量間の相関関係を確認し、その結果をヒートマップとしてグラフ化します。

1
2
3

correlation_matrix = np.corrcoef(X.T)
sns.heatmap(correlation_matrix, annot=True)
plt.show()

このコードは、心臓疾患の予測を行うための機械学習モデルを作成し、その結果をグラフ化します。

このコードは、心臓疾患の予測や他の医学的な問題の解決に有用です。

結果解説

[実行結果]

得られた結果は、RandomForestClassifierを使用して心臓疾患の予測を行った結果の精度を示しています。
精度は、予測が正しい割合を表しています。
したがって、Accuracy: 0.9649122807017544は、このモデルがテストデータに対して正しく予測を行った割合が96.49%であることを示しています。

また、特徴量間の相関関係を確認し、その結果をヒートマップとしてグラフ化しました。
ヒートマップは、データの相関関係を視覚的に表示するためのグラフで、色が暖色であれば高い相関関係を示し、冷色であれば低い相関関係を示します。

このヒートマップを見ると、特徴量間の相関関係を一目で確認できます。
例えば、特徴量Aと特徴量Bの間に高い相関関係があることを示す色が強く塗られている場合、それらの特徴量は予測の結果に大きな影響を与える可能性があります。
これは、特徴量Aと特徴量Bが心臓疾患の予測に大きな影響を与える可能性があることを示しています。

このように、scikit-learnとヒートマップを使用することで、医学的な問題を解決し、その結果を視覚的に理解することが可能です。

ロトカ＝ヴォルテラ方程式

非線形微分方程式を考えましょう。

以下は、ロトカ＝ヴォルテラ方程式と呼ばれる非線形微分方程式の例です。

これをPythonで解き、その解をグラフ化します。

$$
[ \frac{dx}{dt} = \alpha x - \beta xy ]
$$

$$
[ \frac{dy}{dt} = \delta xy - \gamma y ]
$$

ここで、$ (x) $と$ (y) $は時間$ (t) $の関数です。

この方程式は、捕食者と被食者の関係をモデル化したもので、$ (x) $が被食者（例えばウサギ）、$ (y) $が捕食者（例えばキツネ）の個体数を表します。

それでは、この方程式を解いて結果をグラフ化する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

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# ロトカ＝ヴォルテラ方程式
def predator_prey_equations(xy, t, alpha, beta, delta, gamma):
    x, y = xy
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# 初期条件とパラメータ
x0 = 10  # 初期被食者の個体数
y0 = 5   # 初期捕食者の個体数
alpha = 1.1
beta = 0.4
delta = 0.1
gamma = 0.4

# 時間
t = np.linspace(0, 30, 1000)

# 微分方程式を解く
solution = odeint(predator_prey_equations, [x0, y0], t, args=(alpha, beta, delta, gamma))
x, y = solution[:, 0], solution[:, 1]

# 結果をグラフ化
plt.figure(figsize=(8, 6))
plt.plot(t, x, label='被食者 (Prey)')
plt.plot(t, y, label='捕食者 (Predator)')
plt.xlabel('時間')
plt.ylabel('個体数')
plt.title('ロトカ＝ヴォルテラ方程式の解')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、Scipyのodeint関数を使用して微分方程式を数値的に解いています。

結果を時系列でグラフ化し、被食者と捕食者の個体数の変化を示しています。

ソースコード解説

このソースコードは、Pythonでロトカ＝ヴォルテラ方程式（捕食者と被食者の間の相互作用をモデル化する非線形微分方程式）を解き、結果をグラフ化するものです。

それぞれの部分について順番に説明します。

ライブラリのインポート

1
2
3

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

この部分では、NumPy（数値計算用のライブラリ）、Scipyのodeint（常微分方程式を数値的に解くための関数）、Matplotlib（グラフ描画ライブラリ）をインポートしています。

ロトカ＝ヴォルテラ方程式の定義

1
2
3
4
5

def predator_prey_equations(xy, t, alpha, beta, delta, gamma):
    x, y = xy
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

この部分では、ロトカ＝ヴォルテラ方程式を関数predator_prey_equationsとして定義しています。

この方程式は被食者と捕食者の個体数の時間変化を示します。

初期条件とパラメータの設定

1
2
3
4
5
6

x0 = 10  # 初期被食者の個体数
y0 = 5   # 初期捕食者の個体数
alpha = 1.1
beta = 0.4
delta = 0.1
gamma = 0.4

初期条件（被食者と捕食者の初期個体数）と方程式のパラメータ（alpha, beta, delta, gamma）を設定しています。

時間の設定

1

t = np.linspace(0, 30, 1000)

$0$から$30$までの時間を$1000分割$した配列を生成しています。

微分方程式の解を求める

1
2

solution = odeint(predator_prey_equations, [x0, y0], t, args=(alpha, beta, delta, gamma))
x, y = solution[:, 0], solution[:, 1]

odeint関数を使用して微分方程式を数値的に解き、被食者と捕食者の個体数の時間変化を求めています。

結果のグラフ化

1
2
3
4
5
6
7
8
9

plt.figure(figsize=(8, 6))
plt.plot(t, x, label='被食者 (Prey)')
plt.plot(t, y, label='捕食者 (Predator)')
plt.xlabel('時間')
plt.ylabel('個体数')
plt.title('ロトカ＝ヴォルテラ方程式の解')
plt.legend()
plt.grid(True)
plt.show()

Matplotlibを使って、被食者と捕食者の個体数の時間変化をグラフ化しています。

横軸は時間、縦軸は個体数を表し、それぞれのラインは被食者と捕食者を表しています。

グラフは凡例付きで、グリッド線も表示されています。

グラフ解説

このグラフは、ロトカ＝ヴォルテラ方程式に基づいてモデル化された捕食者と被食者の個体数の時間変化を表しています。
この方程式は捕食者と被食者の相互作用を記述し、個体数の変動を表す非線形微分方程式です。

横軸は時間を示し、縦軸は被食者（Prey）と捕食者（Predator）の個体数を表します。
グラフでは、被食者（ウサギなど）の個体数が青色で、捕食者（キツネなど）の個体数がオレンジ色で表示されています。

初期条件では、被食者と捕食者の個体数がそれぞれ$ x0 $と$ y0 $で設定されています。
それぞれの個体数は時間の経過とともに変化し、捕食者は被食者に依存して増減します。
被食者が増えると捕食者も増えるが、その後被食者の減少により捕食者も減少します。
この相互作用が繰り返され、周期的な変動が生じます。

グラフの振動や変動は、捕食者と被食者の間の相互作用とその影響を示しています。
このようなモデルは生態学や生物学などで捕食者と被食者の個体数の相互作用を理解するために利用されます。

Watts-Strogatzモデル

Watts-Strogatzモデルを使ってランダムなグラフを作成し、その次数分布をグラフ化してみましょう。

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

import networkx as nx
import matplotlib.pyplot as plt

# Watts-Strogatzモデルによるグラフの生成
n = 100  # ノード数
k = 4    # 各ノードの近傍ノード数
p = 0.1  # エッジの再接続確率

G = nx.watts_strogatz_graph(n, k, p)

# 次数分布の取得
degrees = [G.degree(node) for node in G.nodes()]

# 次数分布のグラフ化
plt.hist(degrees, bins='auto', alpha=0.7, color='skyblue', edgecolor='black')
plt.xlabel('Degree')
plt.ylabel('Frequency')
plt.title('Degree Distribution of a Watts-Strogatz Graph')
plt.show()

このコードでは、Watts-Strogatzモデルを使用して100個のノードを持つランダムなグラフを生成しています。

その後、各ノードの次数（接続数）を取得し、次数分布をヒストグラムで描画しています。

次数分布はネットワークの特性を示す重要な指標の1つです。

これにより、ノードがどれだけ他のノードと接続しているかの様子を可視化することができます。

ソースコード解説

このコードは、PythonのNetworkXライブラリとMatplotlibを使用して、Watts-Strogatzモデルに基づくランダムなグラフを生成し、その次数分布を可視化しています。

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

1
2

import networkx as nx
import matplotlib.pyplot as plt

• networkxはグラフ理論や複雑ネットワークの作成や操作を行うためのライブラリです。
• matplotlib.pyplotはグラフの描画に使用されます。

2. Watts-Strogatzモデルによるグラフの生成:

1
2
3
4
5

n = 100  # ノード数
k = 4    # 各ノードの近傍ノード数
p = 0.1  # エッジの再接続確率

G = nx.watts_strogatz_graph(n, k, p)

• nx.watts_strogatz_graph(n, k, p)はWatts-Strogatzモデルに基づいてランダムなグラフを生成します。
• nはノード数、kは各ノードの近傍ノード数、pはエッジの再接続確率を表します。

3. 次数分布の取得:

1

degrees = [G.degree(node) for node in G.nodes()]

• G.degree(node)はネットワーク内の各ノードの次数（接続数）を取得します。
• G.nodes()はグラフ内のすべてのノードを返します。
  これを利用して各ノードの次数をリストに収集します。

4. 次数分布のグラフ化:

1
2
3
4
5

plt.hist(degrees, bins='auto', alpha=0.7, color='skyblue', edgecolor='black')
plt.xlabel('Degree')
plt.ylabel('Frequency')
plt.title('Degree Distribution of a Watts-Strogatz Graph')
plt.show()

• plt.hist()は与えられたデータ（degrees）のヒストグラムを作成します。
• bins='auto'は適切なビンの数を自動的に決定します。
• alpha、color、edgecolorはヒストグラムのスタイルを設定します。
• plt.xlabel()、plt.ylabel()、plt.title()はそれぞれx軸のラベル、y軸のラベル、タイトルを設定します。
• plt.show()はグラフを表示します。

これにより、Watts-Strogatzモデルに基づくランダムグラフの次数分布が可視化されます。

次数分布は、ネットワーク内の各ノードが持つ接続の数を表し、ネットワークの構造を理解するのに役立ちます。

最適化問題

PuLPを使用して最適化問題を解いてみましょう。

例として、複数の製品を生産する工場があり、利益を最大化するためにそれぞれの製品の生産量を決定したいとします。

以下は、3つの製品 (A, B, C) の生産にかかる時間と利益、および生産できる最大数量の制約条件を考慮した最適化問題です。

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

from pulp import LpMaximize, LpProblem, LpVariable

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

# 変数の定義
prod_A = LpVariable(name="product_A", lowBound=0, cat="Integer")
prod_B = LpVariable(name="product_B", lowBound=0, cat="Integer")
prod_C = LpVariable(name="product_C", lowBound=0, cat="Integer")

# 利益と制約条件の設定
model += 300 * prod_A + 200 * prod_B + 500 * prod_C, "total_profit"

model += 2 * prod_A + 3 * prod_B + 4 * prod_C <= 100  # 生産時間の制約
model += prod_A + prod_B <= 40  # 製品Aと製品Bの同時生産の制約

# 問題の解決
model.solve()

# 結果の出力
print(f"Product A: {prod_A.varValue}")
print(f"Product B: {prod_B.varValue}")
print(f"Product C: {prod_C.varValue}")
print(f"Total Profit: {model.objective.value()}")

# グラフ化
import matplotlib.pyplot as plt

products = ['Product A', 'Product B', 'Product C']
production_quantities = [prod_A.varValue, prod_B.varValue, prod_C.varValue]

plt.figure(figsize=(8, 5))
plt.bar(products, production_quantities, color='skyblue')
plt.xlabel('Products')
plt.ylabel('Production Quantity')
plt.title('Optimal Production Quantities')
plt.show()

このコードを実行すると、最適な製品の生産量とその利益が表示され、それらを棒グラフで可視化することができます。

このように、PuLPを使用して最適化問題を解き、結果を分かりやすくグラフ化することができます。

[実行結果]

ソースコード解説

このコードは、PuLPとMatplotlibを使用して製品生産最適化問題を解き、結果をグラフ化するものです。

1. PuLPの利用

• from pulp import LpMaximize, LpProblem, LpVariable:
  PuLPの必要なクラスや関数をインポートしています。
• model = LpProblem(name="product_production", sense=LpMaximize):
  問題を定義するための空のモデルを作成しています。
  ここでは、最大化問題（LpMaximize）としてproduct_productionという名前の問題を作成しています。

2. 変数の定義

• prod_A, prod_B, prod_C:
  製品A、B、Cのそれぞれの生産量を表す変数を定義しています。
  LpVariableを使用して整数として設定し、0以上の値を取るように制約を設けています。

3. 利益と制約条件の設定

• model += 300 * prod_A + 200 * prod_B + 500 * prod_C, "total_profit":
  目的関数を定義しています。
  製品ごとの利益を考慮して、それらの生産量の総利益を最大化するように設定しています。
• model += 2 * prod_A + 3 * prod_B + 4 * prod_C <= 100:
  生産時間の制約条件を設定しています。
  各製品の生産にかかる時間が合計で100以下であることを意味します。
• model += prod_A + prod_B <= 40:
  製品AとBの同時生産量に関する制約条件を設定しています。
  製品AとBの合計生産量が40以下であることを要求しています。

4. 問題の解決

• model.solve():
  定義された最適化問題を解きます。

5. 結果の出力とグラフ化

• 最適解から得られた製品ごとの生産量と総利益を出力しています。
• matplotlibを使用して、製品ごとの最適生産量を棒グラフとして可視化しています。

これにより、PuLPを使って製品生産の最適化問題を定義し、解決し、その結果をグラフ化するプロセスが示されています。

結果解説

この問題の解は以下のようになりました：

[実行結果]

• 製品Aを40個生産し、製品Bは0個、製品Cは5個生産されます。
• これにより、総利益は14500になります。

グラフでは、製品ごとの生産量が棒グラフとして示されています。

製品Aが最大生産量である40で、製品Bは制約条件により生産されず、製品Cは5個生産されています。

この生産量のもとで最大の利益である14500が得られます。

製品Bが生産されない理由は、製品Aと製品Bの同時生産に関する制約条件を考慮したためです。

製品Cは比較的高い利益が得られるため、最適解では少量生産されています。

食糧供給最適化

Google OR-Toolsは最適化問題を解くためのライブラリであり、主にルーティング、ネットワークフロー、線形プログラミングなどの問題を扱います。

具体的な貧困問題として、一定の予算で最大数の家庭に食糧を供給するという問題を考えてみます。

この問題は整数線形プログラミング問題として定式化で、Google OR-ToolsのLinearSolverを用いて解くことができます。

それぞれの家庭が必要とする食糧の量と、運搬コスト（距離や時間に基づいた値）が既知であるとします。

以下に、この問題を解くための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 ortools.linear_solver import pywraplp
import matplotlib.pyplot as plt

# データ
families = [0, 1, 2, 3, 4]  # 表現の一部として家族を表すインデックス
food_required = [3, 7, 2, 5, 4]  # 各家族が必要とする食糧の量
transport_costs = [10, 20, 5, 15, 10]  # 各家族への食糧配達のコスト

# ベクトル全体で最小化するためのスカラー乗数
scalar_multiplier = 10

# 問題を作成
solver = pywraplp.Solver.CreateSolver('SCIP')

# 変数を作成
x = []
for i in families:
    x.append(solver.IntVar(0, 1, 'x[%i]' % i))  # x[i] = 家族iに食糧を配達するかどうか (0 または 1)

# 制約を定義
budget = 30  # 予算
solver.Add(sum(transport_costs[i] * x[i] for i in families) <= budget)

# 目的関数を定義
solver.Maximize(scalar_multiplier * sum(food_required[i] * x[i] for i in families) - sum(transport_costs[i] * x[i] for i in families))

# 問題を解く
solver.Solve()

# 結果を出力
print("Total food delivered = ", sum(food_required[i] * x[i].solution_value() for i in families))
delivered = [x[i].solution_value() for i in families]
plt.figure(figsize=(10,5))
plt.bar(range(len(families)), delivered)
plt.title("Food delivered to each family")
plt.xlabel("Family")
plt.ylabel("Food quantity")
plt.xticks(range(len(families)), families)
plt.show()

[実行結果]

ソースコード解説

このコードは、Google OR-Toolsを使用して最適な食糧配達問題を解決し、結果を可視化するものです。

詳細について見ていきましょう。

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

1
2
3
4
5
6
7

from ortools.linear_solver import pywraplp
import matplotlib.pyplot as plt

families = [0, 1, 2, 3, 4]  # 家族を表すインデックス
food_required = [3, 7, 2, 5, 4]  # 各家族の食糧必要量
transport_costs = [10, 20, 5, 15, 10]  # 各家族への食糧配達コスト
scalar_multiplier = 10  # 最小化するためのスカラー乗数

まず、Google OR-ToolsとMatplotlibをインポートし、最適化のためのデータを定義しています。
familiesは家族を表すインデックスのリスト、food_requiredは各家族の食糧必要量のリスト、transport_costsは各家族への食糧配達コストのリストです。
また、scalar_multiplierは後で使われるスカラー乗数です。

2. 問題の作成と変数の定義

1
2
3
4
5

solver = pywraplp.Solver.CreateSolver('SCIP')

x = []
for i in families:
   x.append(solver.IntVar(0, 1, 'x[%i]' % i))  # 家族iに食糧を配達するかどうか (0 または 1)

Google OR-Toolsの線形ソルバーを作成し、バイナリ変数 x を定義しています。
x[i] は家族 i に食糧を配達するかどうかを示す変数で、$0$か$1$の値を取ります。

3. 制約の定義

1
2

budget = 30
solver.Add(sum(transport_costs[i] * x[i] for i in families) <= budget)

予算制約を設定しています。
すべての家族への食糧配達コストの合計が budget 以下である必要があります。

4. 目的関数の定義と最適化

1
2

solver.Maximize(scalar_multiplier * sum(food_required[i] * x[i] for i in families) - sum(transport_costs[i] * x[i] for i in families))
solver.Solve()

目的関数は、家族への食糧配達量を最大化し、同時にコストを最小化するように設定されています。
そして、ソルバーを使って問題を解決しています。

5. 結果の出力と可視化

1
2
3
4
5
6
7
8
9

print("Total food delivered = ", sum(food_required[i] * x[i].solution_value() for i in families))
delivered = [x[i].solution_value() for i in families]
plt.figure(figsize=(10,5))
plt.bar(range(len(families)), delivered)
plt.title("Food delivered to each family")
plt.xlabel("Family")
plt.ylabel("Food quantity")
plt.xticks(range(len(families)), families)
plt.show()

最後に、配達された食糧量を表示し、Matplotlibを使用して各家族への食糧配達量を可視化しています。
print文で配達された総食糧量を表示し、棒グラフで各家族の配達量を示しています。

[実行結果]