Python を用いて情報工学基礎を実行する

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

本演習群ではnumpy・scipy・matplotlib・scikit-learnを使用する。いずれも標準ライブラリではないためインストールが必要である。本手順ではpipに--no-userオプションを付け、ユーザ専用ディレクトリではなく標準のsite-packagesへインストールする（--userの対義の指定）。外部のデータファイルは本演習群では使用しない（実在データはscikit-learn内蔵のものを使う）。

管理者権限でコマンドプロンプトを起動する （手順：Windowsキーまたはスタートメニュー → cmd と入力 → 右クリック → 「管理者として実行」）。 そして、コマンド全体をコマンドプロンプトにコピー＆ペーストする。

pip install -U --no-user numpy scipy matplotlib scikit-learn
python -c "import numpy, scipy, matplotlib, sklearn; print(numpy.__version__, scipy.__version__, matplotlib.__version__, sklearn.__version__)"

実行結果の例（バージョン番号は環境によって異なる）

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。2つの数値の列（ベクトル）x,yに対し、dot(x,y)=Σxᵢyᵢ（各要素どうしを掛けて全部足す）を自分で実装し、numpy.dotの結果と比較する。

# coding: utf-8
import numpy as np


def dot(x, y):
    s = 0.0
    for xi, yi in zip(x, y):
        s += xi * yi
    return s


x = [1, 2, 3, 4, 5]
y = [2, 0, -1, 3, 1]

manual = dot(x, y)
by_numpy = float(np.dot(x, y))

print("自作のdot:", manual)
print("numpy.dot:", by_numpy)
print("一致しているか:", manual == by_numpy)

実行結果の例

ヒント：dot(x,y)は、各要素どうしを掛けてから全部足す操作である。1×2+2×0+3×(-1)+4×3+5×1=2+0-3+12+5=16。この1つの演算が、この先の多くの演習で「平均」「分散」「共分散」「距離」「畳み込み」「スペクトル」という異なる名前で繰り返し登場する。

考察ポイント：内積は「2つの数値の列がどれだけ同じ方向を向いているか（似ているか）」を測る演算である。平均は「データと“全部1”という列」の内積、標本分散は「データ自身との内積」、標本共分散は「2つのデータの内積」、距離は「差との内積」、畳み込みは「核（カーネル）との内積」、スペクトルは「基準の波との内積」である。本演習群を通じて、その都度この対応を確認していく。

3. 演習2：平均

テーマ：平均が「データベクトルと、全要素が1のベクトルとの内積を、要素数で割ったもの」であることを、内積の計算として確認する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。ランダムに8個の数値を生成し、その平均を(1)sum/n、(2)内積を使った計算、(3)numpy.meanの3つの方法で求め、一致することを確認する。

# coding: utf-8
import random
import numpy as np


def dot(x, y):
    s = 0.0
    for xi, yi in zip(x, y):
        s += xi * yi
    return s


random.seed(101)
x = [round(random.uniform(0, 10), 3) for _ in range(8)]
n = len(x)
ones = [1.0] * n

mean_by_sum = sum(x) / n
mean_by_dot = dot(x, ones) / n
mean_by_numpy = float(np.mean(x))

print("データ:", x)
print("sum/nで求めた平均:", mean_by_sum)
print("内積で求めた平均  :", mean_by_dot)
print("numpy.meanの平均  :", mean_by_numpy)

実行結果の例

ヒント：dot(x, ones)は、x₁×1+x₂×1+…+xₙ×1=Σxᵢ、つまり単純な合計と同じ値になる。それをnで割れば平均になる。つまり平均は「全部1のベクトルとの内積」をnで割ったものである。

考察ポイント：平均という、誰もが知っている計算が、実は内積の特別な場合（相手のベクトルが「全部1」であるだけ）だったことを確認した。演習3・演習4では、相手のベクトルを変えることで、同じ内積から標本分散・標本共分散が得られることを見る。

4. 演習3：標本分散

テーマ：標本分散が「データを平均で中心化したベクトルと、その自分自身との内積を、要素数で割ったもの」であることを確認する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。8個の数値を生成し、各要素から平均を引いた「中心化したベクトル」を作り、それと自分自身の内積をnで割って分散を求め、numpy.varと比較する。

# coding: utf-8
import random
import numpy as np


def dot(x, y):
    s = 0.0
    for xi, yi in zip(x, y):
        s += xi * yi
    return s


random.seed(101)
x = [round(random.uniform(0, 10), 3) for _ in range(8)]
n = len(x)
mean_x = sum(x) / n

centered = [xi - mean_x for xi in x]
var_by_dot = dot(centered, centered) / n
var_by_numpy = float(np.var(x))  # ddof=0（標本分散。n-1で割る不偏分散とは異なる）

print("平均:", mean_x)
print("中心化したデータ:", [round(c, 3) for c in centered])
print("内積で求めた標本分散:", var_by_dot)
print("numpy.varの標本分散 :", var_by_numpy)

実行結果の例

ヒント：中心化したベクトルと自分自身の内積はΣ(xᵢ-平均)²、つまり「各データが平均からどれだけ離れているか」の2乗を全部足したものである。これをnで割ったのが標本分散である（numpy.varの既定値はこのnで割る定義であり、n-1で割る「不偏分散」とは異なるので注意する）。

考察ポイント：分散は「自分自身との内積」であり、内積は「似ている度合い」を測る演算だったから、分散は「データが平均からどれだけ広がっているか」を測っていることになる。演習2との違いは、相手のベクトルが「全部1」から「中心化した自分自身」に変わっただけである。

5. 演習4：標本共分散

テーマ：標本共分散が「2つの変数をそれぞれ中心化したベクトルどうしの内積を、要素数で割ったもの」であることを確認し、無関係なデータと相関のあるデータで値がどう変わるかを比較する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。2組のデータ対を用意する。1組目は独立な乱数どうし、2組目は一方からほぼ決まるように作った組である。それぞれの共分散を内積で計算し、numpy.covと比較する。

# coding: utf-8
import random
import numpy as np


def dot(x, y):
    s = 0.0
    for xi, yi in zip(x, y):
        s += xi * yi
    return s


def mean_(x):
    return sum(x) / len(x)


def cov_(x, y):
    n = len(x)
    mx, my = mean_(x), mean_(y)
    xc = [xi - mx for xi in x]
    yc = [yi - my for yi in y]
    return dot(xc, yc) / n


random.seed(102)
a = [round(random.uniform(0, 10), 3) for _ in range(12)]
b_indep = [round(random.uniform(0, 10), 3) for _ in range(12)]
b_corr = [round(2 * ai + random.uniform(-1.5, 1.5), 3) for ai in a]

cov_indep = cov_(a, b_indep)
cov_corr = cov_(a, b_corr)
cov_indep_np = float(np.cov(a, b_indep, ddof=0)[0, 1])
cov_corr_np = float(np.cov(a, b_corr, ddof=0)[0, 1])

print("a              :", a)
print("b(無関係)       :", b_indep)
print("b(相関有り,2a+雑音):", b_corr)
print("共分散(無関係)   内積:", cov_indep, " numpy:", cov_indep_np)
print("共分散(相関有り) 内積:", cov_corr, " numpy:", cov_corr_np)

実行結果の例

ヒント：共分散は標本分散（演習3）の拡張で、「自分自身との内積」を「もう一方のデータとの内積」に変えただけである。一方が増えると他方も増える傾向（正の相関）があれば共分散は正の大きな値になり、無関係ならゼロに近い値になる。

考察ポイント：b_indepはaと無関係に生成したので共分散はゼロに近い値（ここでは-1.42）になり、b_corrは2×aを基準に作ったので共分散は正の大きな値（9.31）になった。平均・標本分散・標本共分散がすべて内積で表せることを、ここまでで確認した。

6. 演習5：日時（西暦）

テーマ：西暦による日時を「年内日数」「時刻」という2つの数値に変換する方法（符号化）を確認し、指定した期間内でランダムな日時を生成する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。date_set関数は、指定した期間内で一様乱数により日付を選び、さらに時・分・秒も一様乱数で生成して、ランダムな日時を1つ返す。symbol_set関数は、指定した比率（重み）に従ってラベル（ここでは「平日」「土曜」「日曜」）をランダムに1つ選んで返す。生成した日時をto_numeric_datetimeで「年内日数 t」「時刻 h」という2つの数値に変換する。

# coding: utf-8
import random
import datetime


def date_set(yy1, mm1, dd1, yy2, mm2, dd2):
    date1 = datetime.date(yy1, mm1, dd1)
    date2 = datetime.date(yy2, mm2, dd2)
    length = (date2 - date1).days + 1
    d = date1 + datetime.timedelta(days=int(length * random.random()))
    result = datetime.datetime(
        d.year, d.month, d.day,
        int(24 * random.random()),
        int(60 * random.random()),
        int(60 * random.random())
    )
    return result


def symbol_set(s, d):
    c = [0] * (len(d) + 1)
    for i in range(1, len(d) + 1):
        c[i] = c[i - 1] + d[i - 1]
    m = c[len(d)]
    t = [None] * (m + 1)
    t[0] = -1
    for i in range(1, len(d) + 1):
        for j in range(c[i - 1], c[i]):
            t[j + 1] = i - 1
    idx = int(m * random.random()) + 1
    return s[t[idx]]


def to_numeric_datetime(dt):
    # 西暦による日時を「年内日数 t」「時刻 h」の2つの数値に変換する（=符号化）
    day_of_year = (dt.date() - datetime.date(dt.year, 1, 1)).days
    hour = dt.hour + dt.minute / 60 + dt.second / 3600
    return day_of_year, hour


def cov_(x, y):
    n = len(x)
    mx, my = sum(x) / n, sum(y) / n
    return sum((xi - mx) * (yi - my) for xi, yi in zip(x, y)) / n


def corr_(x, y):
    return cov_(x, y) / ((cov_(x, x) ** 0.5) * (cov_(y, y) ** 0.5))


random.seed(5)
print("--- 5回分の生成例 ---")
for _ in range(5):
    dt = date_set(2012, 4, 1, 2013, 3, 31)
    label = symbol_set(["平日", "土曜", "日曜"], [5, 1, 1])
    t, h = to_numeric_datetime(dt)
    print(dt, label, "-> t =", t, ", h =", round(h, 3))

# tとhはdate_setの中で別々の乱数で決めているので、無関係(無相関)のはずである。
# サンプル数を増やすと相関係数がどう変化するかを確認する。
for seed, N in [(55, 300), (55, 3000)]:
    random.seed(seed)
    ts, hs = [], []
    for _ in range(N):
        dt = date_set(2012, 4, 1, 2013, 3, 31)
        t, h = to_numeric_datetime(dt)
        ts.append(t)
        hs.append(h)
    print(f"N={N} のときの相関係数(t,h):", round(corr_(ts, hs), 4))

実行結果の例

ヒント：コンピュータは「日時」を直接理解しているわけではなく、内部的にはすべて数値として扱っている。「年内日数」「時刻」という2つの数値に変換する方法はその一例であり、これ以外にも「ある基準日から何日目か（演習6で扱う）」のように1つの整数に変換する方法もある。どちらも西暦という日時の表現を、別の数値表現に変換し直しているだけである。

考察ポイント：t(年内日数)とh(時刻)は、date_setの中でそれぞれ独立な乱数で決めているため、理論上は無関係（無相関＝相関係数0）であるはずである。実際にN=300個のサンプルで相関係数を計算すると-0.125とゼロから少しずれるが、これは標本の数が少ないことによる誤差（ばらつき）である。N=3000に増やすと-0.020まで0に近づく。これは演習4で扱った「標本共分散」（から計算される相関係数）を、サンプル数を増やしながら確認した例である。

7. 演習6：4バイト整数のコード化

テーマ：日時を表す1つの整数（基準日からの日数）を、コンピュータが実際にメモリ上で使う「4バイト（32ビット）整数」という形に変換し、そのビット列を直接確認する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。期間内でランダムな日時を1つ生成し、date.toordinal()で「西暦1年1月1日を1日目とする日数（シリアル値）」という1つの整数に変換する。これをstruct.pack("<i", 値)で4バイトのバイト列に変換し、16進数で表示する。負の数を4バイト整数にした場合の見た目（2の補数表現）も確認する。

# coding: utf-8
import random
import datetime
import struct


def date_set(yy1, mm1, dd1, yy2, mm2, dd2):
    date1 = datetime.date(yy1, mm1, dd1)
    date2 = datetime.date(yy2, mm2, dd2)
    length = (date2 - date1).days + 1
    d = date1 + datetime.timedelta(days=int(length * random.random()))
    return datetime.datetime(
        d.year, d.month, d.day,
        int(24 * random.random()), int(60 * random.random()), int(60 * random.random())
    )


random.seed(5)
dt = date_set(2012, 4, 1, 2013, 3, 31)
serial = dt.date().toordinal()  # 西暦0001-01-01を1とする日数

bytes_le = struct.pack('<i', serial)   # リトルエンディアン（下位バイトが先）
bytes_be = struct.pack('>i', serial)   # ビッグエンディアン（上位バイトが先）
back = struct.unpack('<i', bytes_le)[0]

print("日時          :", dt)
print("シリアル値(整数):", serial)
print("4バイト(リトルエンディアン):", bytes_le.hex())
print("4バイト(ビッグエンディアン):", bytes_be.hex())
print("バイト列から復元した値    :", back, "（元の値と一致:", back == serial, "）")
print("-1 を4バイト整数にした場合:", struct.pack('<i', -1).hex(), "（2の補数表現）")

実行結果の例

ヒント：4バイト＝32ビットの符号付き整数は、最大で約21億（2の31乗-1）まで表せる。日付のシリアル値（西暦1年から数えた日数）は数十万程度なので十分に収まる。-1がffffffff（すべてのビットが1）になるのは、負の数を「2の補数」という方式で表しているためである。リトルエンディアンとビッグエンディアンの違いは、バイトを並べる順序（下位バイトを先に置くか、上位バイトを先に置くか）の違いであり、CPUの種類によって標準が異なる。

考察ポイント：演習5では日時を「年内日数」「時刻」という2つの実数として表現したが、ここでは「基準日からの日数」という1つの整数として表現し、それを実際にコンピュータが使う4バイトのビット列に変換した。同じ「日時」というデータでも、目的に応じて複数の数値表現（符号化方式）があり得ることが分かる。

8. 演習7：文字コード

テーマ：文字（テキスト）データが、コンピュータの中では数値（バイト列）として符号化されていることを、日本語と英語の文字数・バイト数の違いから確認する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。日本語の文字列をUTF-8という方式でバイト列に変換し、文字数とバイト数を比較する。英語の文字列とも比較する。

# coding: utf-8
label = "土曜"
en_text = "Weekday"

label_bytes = label.encode('utf-8')
en_bytes = en_text.encode('ascii')

print("日本語:", label, "/ 文字数:", len(label), "/ UTF-8のバイト数:", len(label_bytes))
print("バイト列(16進):", label_bytes.hex())
print("英語  :", en_text, "/ 文字数:", len(en_text), "/ ASCIIのバイト数:", len(en_bytes))

実行結果の例

ヒント：ASCII（英数字や記号）は1文字を1バイトで表せるが、UTF-8という符号化方式では、日本語の文字（漢字・かな）の多くは1文字を3バイトで表す。そのため、「土曜」(2文字)は2×3=6バイトになる。英語は1文字につき1バイトなので、文字数とバイト数が一致する。

考察ポイント：演習6では数値（整数）の符号化を、ここでは文字の符号化を見た。どちらも「人間にとって意味のある情報を、有限長のバイト列に変換する」という同じ目的を持つが、変換の方式（符号化方式）はデータの種類ごとに異なる。同じ日本語の文章でも、文字数だけでは実際の保存に必要なバイト数（データ量）は分からない、という点は実用上重要である。

9. 演習8：緯度・経度

テーマ：緯度・経度という2つの数値の組で「場所」を表す方法を確認し、ランダムに複数の点を生成して散布図に描く。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。緯度・経度はそれぞれ単独の数値であり、その組(経度,緯度)が地球上の1点を表す。次のコードでは、ある中心点のまわりに、ある方向に偏って広がるように80個の点をランダムに生成し、散布図としてex8_scatter.pngに保存・表示する。

# coding: utf-8
import random
import math
import os
import matplotlib.pyplot as plt

# Windows環境での日本語フォントの文字化け（警告）対策
plt.rcParams['font.family'] = 'Meiryo'

random.seed(8)
N = 80
theta = math.radians(35)  # 点が広がる方向（基準点からの角度）

points = []
for _ in range(N):
    s = random.uniform(-1, 1) * 0.25   # 主な広がりの方向に沿った位置
    u = random.uniform(-1, 1) * 0.05   # それと垂直な方向の小さなばらつき
    dx = s * math.cos(theta) - u * math.sin(theta)
    dy = s * math.sin(theta) + u * math.cos(theta)
    lon = 137.00 + dx   # 経度
    lat = 35.00 + dy    # 緯度
    points.append((lon, lat))

lons = [p[0] for p in points]
lats = [p[1] for p in points]

print("生成した点の数:", N)
for lon, lat in points[:5]:
    print("経度=%.5f, 緯度=%.5f" % (lon, lat))

plt.figure()
plt.scatter(lons, lats, s=20)
plt.xlabel("経度 (lon)")
plt.ylabel("緯度 (lat)")
plt.title("生成した点(経度,緯度)の散布図")
plt.gca().set_aspect("equal")

# プログラムファイルと同じディレクトリのパスを取得して保存
script_dir = os.path.dirname(os.path.abspath(__file__))
save_path = os.path.join(script_dir, "ex8_scatter.png")
plt.savefig(save_path)

plt.show()

実行結果の例

ヒント：緯度・経度はどちらも「数値」であり、特別な型ではない。地球上の場所は、結局「2つの実数の組」として表現されている。今回はある方向（35度方向）に偏って点を広げて生成している。

考察ポイント：散布図を見ると、点が斜め方向に偏って分布していることが分かる。この「偏り」を目で見て確認することはできるが、これを数値として表す方法（共分散・主成分分析）は演習10で扱う。

10. 演習9：浮動小数点数（8バイト）

テーマ：緯度・経度のような小数を含む数値が、コンピュータの中で「8バイト（64ビット）の浮動小数点数（IEEE754形式）」としてどのように符号化されているかを確認する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。経度を表す小数の値を1つ取り、8バイトのバイト列に変換する。さらに、そのビット列を符号・指数部・仮数部の3つに分けて表示する。

# coding: utf-8
import struct

v = 136.86154911909023   # 経度を表す小数の値の例

b = struct.pack("<d", v)
back = struct.unpack("<d", b)[0]

bits = struct.unpack(">Q", struct.pack(">d", v))[0]
binstr = format(bits, "064b")
sign, exponent, mantissa = binstr[0], binstr[1:12], binstr[12:]

print("値          :", v)
print("8バイト(16進):", b.hex())
print("復元した値  :", back, "（元の値と一致:", back == v, "）")
print("符号部(1bit) :", sign)
print("指数部(11bit):", exponent)
print("仮数部(52bit):", mantissa)

実行結果の例

ヒント：IEEE754という規格では、64ビット（8バイト）を「符号1ビット」「指数部11ビット」「仮数部52ビット」に分けて実数を表す。整数（演習6）が「そのままの値」をビットに変換するのに対し、浮動小数点数は「仮数部 × 2の指数部乗」という形に分解してから変換する点が異なる。

考察ポイント：整数の符号化（演習6、4バイト）と浮動小数点数の符号化（本演習、8バイト）は、同じ「数値をバイト列にする」という目的でも仕組みが異なる。緯度・経度のような小数を含む量は浮動小数点数として、日付のシリアル値のような整数量は整数として、それぞれ別の方式で符号化されている。

11. 演習10：２次元

テーマ：実在のデータ（scikit-learn内蔵のiris）から2つの特徴量を取り出し、クラスごとに色分けした散布図を描く。そのうえで、2次元データに特徴的な分析として、2×2の共分散行列・相関係数と、主成分分析（PCA）による「最もデータが広がっている方向（主軸）」を求め、その主軸を散布図に重ねて描く。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。irisの150個の花について、特徴量のうち「花びらの長さ（petal length）」と「花びらの幅（petal width）」の2つを取り出す。(1)3つの品種（クラス）ごとに色を変えて散布図を描く、(2)2×2の共分散行列と相関係数を計算する、(3)主成分分析で第1主成分（最もばらつきが大きい方向）を求め、その方向の直線を散布図の上に重ねてex10_iris_2d.pngに保存・表示する。

# coding: utf-8
# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import os
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

# Windows環境での日本語フォントの文字化け対策
plt.rcParams['font.family'] = 'Meiryo'

iris = load_iris()
# petal length (列2), petal width (列3) の2特徴量を使う
X = iris.data[:, [2, 3]]
y = iris.target
names = iris.target_names

# (1) クラスごとに色分けした散布図
plt.figure()
for k in range(len(names)):
    plt.scatter(X[y == k, 0], X[y == k, 1], s=20, label=names[k])

# (2) 共分散行列と相関係数
cov = np.cov(X.T, ddof=0)
corr = cov[0, 1] / np.sqrt(cov[0, 0] * cov[1, 1])
print("2x2 共分散行列:")
print(cov)
print("相関係数(petal length, petal width):", corr)

# (3) 主成分分析（第1主成分=最もばらつきが大きい方向）
pca = PCA(n_components=2)
pca.fit(X)
center = X.mean(axis=0)
pc1 = pca.components_[0]                 # 第1主成分の方向(単位ベクトル)
print("各主成分が説明する分散の割合:", pca.explained_variance_ratio_)
print("第1主成分の方向:", pc1)

# 第1主成分の方向の直線を散布図に重ねる
length = 3.0 * np.sqrt(pca.explained_variance_[0])
line = np.array([center - length * pc1, center + length * pc1])
plt.plot(line[:, 0], line[:, 1], color="black", linewidth=2, label="第1主成分")

plt.xlabel("petal length")
plt.ylabel("petal width")
plt.title("iris 2特徴量の散布図と第1主成分")
plt.legend()

# プログラムファイルと同じディレクトリのパスを取得して保存
script_dir = os.path.dirname(os.path.abspath(__file__))
save_path = os.path.join(script_dir, "ex10_iris_2d.png")
plt.savefig(save_path)

# 図を表示
plt.show()

実行結果の例

ヒント：共分散行列の対角成分(3.09, 0.58)は各特徴量の分散、非対角成分(1.29)は2つの特徴量の共分散である。相関係数が約0.96と1に近いことは、花びらの長さと幅がほぼ比例して増減することを意味する。主成分分析は、この「データが最もよく広がっている方向」を1本の軸（第1主成分）として取り出す手法であり、第1主成分だけで全体のばらつきの約99%（0.990）を説明できている。

考察ポイント：2次元データに特徴的なのは、ばらつきが「縦の分散」「横の分散」だけでなく「斜め方向の共分散」も持つ点である。相関が強い（散布図が斜めに伸びている）ほど、第1主成分1本でデータの大部分を表せる。これは、2つの数値で表したデータが、見かけ上は2次元でも、本質的にはほぼ1次元の広がりしか持たない場合があることを示している。色分けを見ると、setosaが他の2品種からはっきり離れて分布していることも読み取れる。

12. 演習11：３次元

テーマ：実在のデータ（scikit-learn内蔵のiris）から3つの特徴量を取り出し、クラスごとに色分けした3次元散布図を描く。そのうえで、3次元データに特徴的な分析として、3×3の共分散行列と、主成分分析（PCA）による各主成分の寄与率を求め、3次元のばらつきが実質的に何次元で表せるかを確認する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。irisの150個の花について、「がくの長さ（sepal length）」「花びらの長さ（petal length）」「花びらの幅（petal width）」の3つの特徴量を取り出す。(1)3つの品種ごとに色を変えて3次元散布図を描きex11_iris_3d.pngに保存・表示する、(2)3×3の共分散行列を計算する、(3)主成分分析で3つの主成分それぞれの寄与率（説明する分散の割合）を求める。

# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import os
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401  (3D投影の登録に必要)
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA

# Windows環境での日本語フォントの文字化け対策
plt.rcParams['font.family'] = 'Meiryo'

iris = load_iris()
# sepal length(列0), petal length(列2), petal width(列3) の3特徴量を使う
X = iris.data[:, [0, 2, 3]]
y = iris.target
names = iris.target_names

# (1) クラスごとに色分けした3次元散布図
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
for k in range(len(names)):
    ax.scatter(X[y == k, 0], X[y == k, 1], X[y == k, 2], s=20, label=names[k])
ax.set_xlabel("sepal length")
ax.set_ylabel("petal length")
ax.set_zlabel("petal width")
ax.set_title("iris 3特徴量の3次元散布図")
ax.legend()

# プログラムファイルと同じディレクトリのパスを取得して保存
script_dir = os.path.dirname(os.path.abspath(__file__))
save_path = os.path.join(script_dir, "ex11_iris_3d.png")
plt.savefig(save_path)
plt.show()

# (2) 3x3 共分散行列
cov = np.cov(X.T, ddof=0)
print("3x3 共分散行列:")
print(cov)

# (3) 主成分分析（各主成分の寄与率）
pca = PCA(n_components=3)
pca.fit(X)
print("各主成分が説明する分散の割合:", pca.explained_variance_ratio_)
print("第1+第2主成分の累積寄与率:", np.sum(pca.explained_variance_ratio_[:2]))

実行結果の例

ヒント：3次元では、共分散行列が3×3に増え、特徴量どうしの組み合わせ（sepal length×petal lengthなど）すべてについて共分散を持つ。主成分分析を行うと、第1主成分だけで全体のばらつきの約98%（0.982）、第1+第2主成分で約99.7%（0.997）を説明できる。これは、3次元で表したデータが、実質的にはほぼ2次元（むしろ1次元に近い）平面状の広がりしか持たないことを意味する。

考察ポイント：3次元データに特徴的なのは、見かけの次元（3）と、データが実際に広がっている本質的な次元が一致しないことが、寄与率という数値ではっきり分かる点である。3次元散布図を回して見ると、点がほぼ1枚の平面（さらにはほぼ1本の線）の近くに集まっていることが視覚的にも確認できる。主成分分析は、このような「高次元データの本当の広がりの次元」を測り、少ない軸でデータを要約する道具である。次元が2でも3でも、共分散行列と主成分分析の考え方がそのまま使える点は共通している。

13. 演習12：画像

テーマ：画像が「画素の値を並べた数値の配列（ベクトル）」であることを確認し、画素値の平均（明るさ）・分散（コントラスト）を、演習2・演習3と同じ内積の計算で求める。画像も表示する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。ランダムに置いた点から距離マップ（各画素から最も近い点までの距離）を作り、0〜255の値を持つ「画像」とする。画像を1列に並べ直した「画素値ベクトル」に対して、平均・分散を内積で計算し、numpyの結果と比較する。画像はex12_image.pngに保存・表示する。

# coding: utf-8
import random
import math
import numpy as np
import matplotlib.pyplot as plt
import os
from scipy import ndimage

# Windows環境での日本語フォントの文字化け対策
plt.rcParams['font.family'] = 'Meiryo'

def dot(x, y):
    s = 0.0
    for xi, yi in zip(x, y):
        s += xi * yi
    return s

# 点を生成し、距離マップから画像を作る（この演習だけで完結する）
random.seed(8)
N = 80
theta = math.radians(35)
lons, lats = [], []
for _ in range(N):
    s = random.uniform(-1, 1) * 0.25
    u = random.uniform(-1, 1) * 0.05
    dx = s * math.cos(theta) - u * math.sin(theta)
    dy = s * math.sin(theta) + u * math.cos(theta)
    lons.append(137.00 + dx)
    lats.append(35.00 + dy)

GRID = 256
a, b = min(lons), max(lons) - min(lons)
c, d = min(lats), max(lats) - min(lats)
im = np.ones((GRID, GRID))
for lo, la in zip(lons, lats):
    ix = min(int(GRID * (lo - a) / b), GRID - 1)
    iy = min(int(GRID * (la - c) / d), GRID - 1)
    im[ix, iy] = 0
dm = ndimage.distance_transform_edt(im)
img = (dm / dm.max() * 255).astype(np.uint8)

flat = img.astype(float).flatten().tolist()
n_px = len(flat)
mean_brightness = dot(flat, [1.0] * n_px) / n_px
centered = [p - mean_brightness for p in flat]
var_brightness = dot(centered, centered) / n_px
print("画像のサイズ(配列の形):", img.shape, " 画素数:", n_px)
print("内積で求めた平均(明るさ)    :", mean_brightness, " numpy.mean:", float(np.mean(img)))
print("内積で求めた分散(コントラスト):", var_brightness, " numpy.var:", float(np.var(img)))

plt.figure()
plt.imshow(img, cmap="gray")
plt.title("距離マップを画像として見る")
plt.colorbar()

# プログラムファイルと同じディレクトリのパスを取得して保存
script_dir = os.path.dirname(os.path.abspath(__file__))
save_path = os.path.join(script_dir, "ex12_image.png")
plt.savefig(save_path)

plt.show()

実行結果の例

ヒント：画像は、見た目は2次元の格子状だが、計算上は「すべての画素値を1列に並べたベクトル」として扱うことができる。65536個の画素を持つこの画像も、演習2・演習3で8個の要素を持つベクトルに対して行ったのと全く同じ式（内積を使った平均・分散）で処理できる。

考察ポイント：「明るさ」は画素値の平均、「コントラスト（濃淡の強さ）」は画素値の分散である、という対応が成立する。これらは演習2・演習3で確認した、平均・分散の計算であり、画像だからといって特別な計算が必要になるわけではない。要素数が数個から数万個に増えても、内積で表される計算の形は変わらない。

14. 演習13：フィルタ

テーマ：画像に「フィルタ」（ぼかしフィルタ、エッジ検出フィルタ）を適用し、画像がどのように変化するかを確認する。元画像・ぼかし後・エッジ検出後を並べて表示する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。距離マップから画像を作り、(1)5×5の範囲を平均する「ぼかしフィルタ」と、(2)輪郭を強調する「エッジ検出フィルタ」を、scipy.ndimage.convolveで適用し、結果を比較する。3枚の画像を並べてex13_filter.pngに保存・表示する。

# coding: utf-8
import random
import math
import numpy as np
import matplotlib.pyplot as plt
import os
from scipy import ndimage

# Windows環境での日本語フォントの文字化け対策
plt.rcParams['font.family'] = 'Meiryo'

# 点を生成し、距離マップから画像を作る（この演習だけで完結する）
random.seed(8)
N = 80
theta = math.radians(35)
lons, lats = [], []
for _ in range(N):
    s = random.uniform(-1, 1) * 0.25
    u = random.uniform(-1, 1) * 0.05
    dx = s * math.cos(theta) - u * math.sin(theta)
    dy = s * math.sin(theta) + u * math.cos(theta)
    lons.append(137.00 + dx)
    lats.append(35.00 + dy)

GRID = 256
a, b = min(lons), max(lons) - min(lons)
c, d = min(lats), max(lats) - min(lats)
im = np.ones((GRID, GRID))
for lo, la in zip(lons, lats):
    ix = min(int(GRID * (lo - a) / b), GRID - 1)
    iy = min(int(GRID * (la - c) / d), GRID - 1)
    im[ix, iy] = 0
dm = ndimage.distance_transform_edt(im)
img = (dm / dm.max() * 255).astype(np.uint8).astype(float)

# フィルタ（カーネル）の定義
blur_kernel = np.ones((5, 5)) / 25.0
edge_kernel = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]], dtype=float)

blurred = ndimage.convolve(img, blur_kernel, mode="nearest")
edges = ndimage.convolve(img, edge_kernel, mode="nearest")
print("元画像        の範囲:", img.min(), "〜", img.max())
print("ぼかし後の画像の範囲:", blurred.min(), "〜", blurred.max())
print("エッジ検出後の範囲  :", edges.min(), "〜", edges.max())

fig, axes = plt.subplots(1, 3)
axes[0].imshow(img, cmap="gray")
axes[0].set_title("元画像")
axes[1].imshow(blurred, cmap="gray")
axes[1].set_title("ぼかし")
axes[2].imshow(edges, cmap="gray")
axes[2].set_title("エッジ検出")

# プログラムファイルと同じディレクトリのパスを取得して保存
script_dir = os.path.dirname(os.path.abspath(__file__))
save_path = os.path.join(script_dir, "ex13_filter.png")
plt.savefig(save_path)

plt.show()

実行結果の例

ヒント：「フィルタ」とは、ある決まった大きさの数値の表（カーネル）を画像全体にずらしながら適用していく処理である。ぼかしフィルタのカーネルはすべて等しい正の値（合計1）なので、各画素の値が周囲の画素との平均に置き換わり、画像がなめらかになる。エッジ検出フィルタのカーネルは中心が正・周囲が負なので、周囲と値が違う場所（輪郭）で大きな値が出る。

考察ポイント：「フィルタを適用する」という処理が、各位置でカーネルと画像の一部を使って何らかの計算をしている、ということは想像できるが、その計算が具体的に何をしているのかは、演習14で「畳み込み」として内積に帰着させて確認する。

15. 演習14：畳み込み

テーマ：フィルタ（演習13）の正体が「カーネルを画像の各位置に当てはめて、そのたびに内積を取る」という畳み込み（convolution）であることを、自分で実装して確認する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。カーネルを180度回転（上下左右反転）させたうえで、画像の各位置の小さな範囲（パッチ）との内積を計算する関数を自作し、scipy.ndimage.convolveの結果と一致することを確認する。カーネルを反転させない場合（相関、correlation）との違いも、左右非対称なカーネルで確認する。

# coding: utf-8
import numpy as np
from scipy import ndimage


def dot(x, y):
    s = 0.0
    for xi, yi in zip(x, y):
        s += xi * yi
    return s


def convolve_true(image, kernel, mode="nearest"):
    # 数学的な定義どおり、核を180度回転（上下左右反転）してから
    # 各位置で局所パッチとの内積をとる（=畳み込み）
    kernel_flipped = kernel[::-1, ::-1]
    kh, kw = kernel.shape
    ph, pw = kh // 2, kw // 2
    padded = np.pad(image, ((ph, ph), (pw, pw)), mode="edge")
    out = np.zeros_like(image, dtype=float)
    flat_kernel = kernel_flipped.flatten().tolist()
    H, W = image.shape
    for i in range(H):
        for j in range(W):
            patch = padded[i:i + kh, j:j + kw].flatten().tolist()
            out[i, j] = dot(patch, flat_kernel)
    return out


# 簡単な30x30の画像と、左右非対称なカーネルで確認する
rng = np.random.default_rng(0)
small_img = rng.uniform(0, 255, size=(30, 30))
asym_kernel = np.array([[1, 2, 3], [0, 0, 0], [-1, -2, -3]], dtype=float)

manual = convolve_true(small_img, asym_kernel, mode="nearest")
by_scipy_convolve = ndimage.convolve(small_img, asym_kernel, mode="nearest")
by_scipy_correlate = ndimage.correlate(small_img, asym_kernel, mode="nearest")  # 反転しない場合

diff_with_flip = np.max(np.abs(manual - by_scipy_convolve))
diff_without_flip = np.max(np.abs(manual - by_scipy_correlate))

print("自作(核を反転)とscipy.ndimage.convolveの最大差 :", diff_with_flip)
print("自作(核を反転)とscipy.ndimage.correlateの最大差:", diff_without_flip)

実行結果の例

ヒント：「畳み込み（convolution）」は数学的にはカーネルを180度回転させてから内積をとる演算であり、回転させずに内積をとる演算は「相関（correlation）」と呼ばれ、両者は左右非対称なカーネルでは結果が異なる（演習13のぼかし・エッジ検出フィルタのカーネルは上下左右対称なので、回転させても結果が変わらず、フィルタとしては畳み込みと相関の違いが見えなかった）。

考察ポイント：フィルタ処理の正体は「カーネルとの内積を、画像の各位置で繰り返す」ことだった。これは演習1の考察ポイントで予告した「畳み込みも内積の一種である」という対応そのものである。画像の各画素を1つ計算するたびに、カーネルの大きさ分の小さな内積を1回計算している。

16. 演習15：音声

テーマ：音の波形が「時間軸上に並んだ数値の配列」であることを確認し、コード内で正弦波を合成して波形データを生成し、波形を表示する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。サンプリング周波数（1秒間に何回数値を記録するか）を200Hzとし、2秒間、5Hzの正弦波と18Hzの正弦波を足し合わせた波形を生成し、ex15_wave.pngに保存・表示する。

# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import os

# Windows環境での日本語フォントの文字化け対策
plt.rcParams['font.family'] = 'Meiryo'

fs = 200.0          # サンプリング周波数[Hz]（1秒間に200回、値を記録する）
duration = 2.0       # 録音する長さ[s]
t = np.arange(0, duration, 1.0 / fs)   # 時刻の配列（標本化）

f1, f2 = 5.0, 18.0   # 含めたい2つの周波数[Hz]
signal = 1.0 * np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
print("サンプリング周波数:", fs, "Hz")
print("サンプル数(配列の長さ):", len(t))
print("先頭5個の値:", signal[:5])

plt.figure()
plt.plot(t, signal)
plt.xlabel("時間 [s]")
plt.ylabel("振幅")
plt.title("生成した波形（5Hzと18Hzの正弦波の合成）")

# プログラムファイルと同じディレクトリのパスを取得して保存
script_dir = os.path.dirname(os.path.abspath(__file__))
save_path = os.path.join(script_dir, "ex15_wave.png")
plt.savefig(save_path)

plt.show()

実行結果の例

ヒント：実際の音声（CDやスマートフォン）は44100Hzや48000Hzといった、もっと高いサンプリング周波数で記録されているが、ここでは波形の形がグラフで見やすいように、人間の耳には聞こえない低い周波数（5Hz、18Hz）を使っている。重要なのは「音もまた、画像の画素値（演習12）と同じように、数値の配列である」という点である。

考察ポイント：グラフを見ると、ゆっくりした5Hzの波の上に、速い18Hzの波が重なっている様子が分かる。しかし、配列の数値を1つずつ見ているだけでは「5Hzと18Hzの成分が含まれている」ことは分からない。それを取り出す方法が、演習16の「スペクトル」である。

17. 演習16：スペクトル

テーマ：波形データと、いろいろな周波数の基準となる正弦波・余弦波との内積を計算することで、波形に含まれる周波数の成分（スペクトル）を取り出し、グラフに描く。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。5Hzと18Hzを合成した波形を生成し、1〜40Hzの範囲で候補の周波数を0.5Hzごとに用意する。それぞれの候補周波数について、基準となるcos波・sin波と元の波形との内積を計算し、その大きさをその周波数の「強さ」としてグラフにし、ex16_spectrum.pngに保存・表示する。

# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import os

# Windows環境での日本語フォントの文字化け対策
plt.rcParams['font.family'] = 'Meiryo'

def dot(x, y):
    s = 0.0
    for xi, yi in zip(x, y):
        s += xi * yi
    return s

fs = 200.0
duration = 2.0
t = np.arange(0, duration, 1.0 / fs)
f1, f2 = 5.0, 18.0
signal = 1.0 * np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)

freqs = np.arange(1, 41, 0.5)   # 候補の周波数[Hz]
spectrum = []
for f in freqs:
    ref_cos = np.cos(2 * np.pi * f * t)
    ref_sin = np.sin(2 * np.pi * f * t)
    c = dot(signal.tolist(), ref_cos.tolist())
    s = dot(signal.tolist(), ref_sin.tolist())
    amp = (2.0 / len(t)) * np.sqrt(c ** 2 + s ** 2)   # その周波数の強さ
    spectrum.append(amp)
spectrum = np.array(spectrum)

top2 = np.argsort(spectrum)[::-1][:2]
print("強さが大きい上位2つの周波数:", freqs[top2])
print("そのときの強さ            :", spectrum[top2])

plt.figure()
plt.plot(freqs, spectrum)
plt.xlabel("周波数 [Hz]")
plt.ylabel("強さ")
plt.title("内積で求めたスペクトル")

# プログラムファイルと同じディレクトリのパスを取得して保存
script_dir = os.path.dirname(os.path.abspath(__file__))
save_path = os.path.join(script_dir, "ex16_spectrum.png")
plt.savefig(save_path)

plt.show()

実行結果の例

ヒント：ある周波数のcos波・sin波と元の波形との内積が大きいということは、内積が「似ている度合いを測る演算」（演習1）であったことから、「元の波形にその周波数の成分がたくさん含まれている」ことを意味する。すべての候補周波数についてこの内積を計算してグラフにしたものが「スペクトル」である。

考察ポイント：スペクトルのグラフには、合成に使った5Hzと18Hzの位置に、ちょうど振幅1.0と0.5の大きさのピークが現れる（合成したときの係数と一致する）。これは演習1の考察ポイントで「波の周波数を取り出すのも内積である」と予告した内容そのものである。

18. 演習17：周波数

テーマ：波形のスペクトルから、波形に実際に含まれている周波数を特定し、生成時に指定した周波数と一致することを確認する。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。5Hzと18Hzを合成した波形を生成してスペクトルを求め、スペクトルの中で値が最大になる周波数（と、それを除いたうえでの2番目に大きい周波数）を求め、合成に使った周波数(5Hz、18Hz)と比較する。

# coding: utf-8
import numpy as np


def dot(x, y):
    s = 0.0
    for xi, yi in zip(x, y):
        s += xi * yi
    return s


fs = 200.0
duration = 2.0
t = np.arange(0, duration, 1.0 / fs)
f1, f2 = 5.0, 18.0
signal = 1.0 * np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)

freqs = np.arange(1, 41, 0.5)
spectrum = []
for f in freqs:
    c = dot(signal.tolist(), np.cos(2 * np.pi * f * t).tolist())
    s = dot(signal.tolist(), np.sin(2 * np.pi * f * t).tolist())
    spectrum.append((2.0 / len(t)) * np.sqrt(c ** 2 + s ** 2))
spectrum = np.array(spectrum)

best_idx = int(np.argmax(spectrum))
best_freq = freqs[best_idx]

second = spectrum.copy()
second[best_idx] = -1
second_idx = int(np.argmax(second))
second_freq = freqs[second_idx]

print("最も強い周波数  :", best_freq, "Hz  （実際に使った値:", f1, "Hz）")
print("2番目に強い周波数:", second_freq, "Hz  （実際に使った値:", f2, "Hz）")

実行結果の例

ヒント：「最も強い周波数を求める」処理は、np.argmax、つまり「スペクトルの値が最大になる候補を探す」という最適化（最大化）の操作である。これは演習18で扱う「関数の最適化」の具体例の1つである。

考察ポイント：コンピュータは波形の数値の並びしか持っていないが、内積を使ってスペクトルを計算し、その最大値を探すことで、元の波形を作ったときに使った周波数（5Hz、18Hz）を、誤差なく正確に復元できた。これは「内積」と「最適化」という2つの道具を組み合わせた結果である。

19. 演習18：関数と最適化

テーマ：本演習群に繰り返し登場した「最小化・最大化」という考え方を、「関数の最適化」として明確に取り出し、平均が二乗誤差の最小化として求められることを確認する。誤差関数のグラフも描く。

手順：

次のコードを実行（メモ帳を用いる場合は a.py のようなファイル名で保存して実行）する。8個の数値について、候補となるμをいろいろ試し、二乗誤差J(μ)=Σ(xᵢ-μ)²を計算してグラフにし、ex18_optimize.pngに保存・表示する。J(μ)を最小にするμが、平均(sum/n)と一致することを確認する。そのうえで、画像の距離マップ（最小化）や周波数探索（最大化）も同じ「関数の最適化」であったことを振り返る。

# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import os

# Windows環境での日本語フォントの文字化け対策
plt.rcParams['font.family'] = 'Meiryo'

x = [5.812, 1.948, 9.653, 9.24, 4.671, 6.635, 2.145, 2.217]
true_mean = sum(x) / len(x)

candidates = np.linspace(min(x) - 2, max(x) + 2, 400)   # μの候補
J = np.array([sum((xi - mu) ** 2 for xi in x) for mu in candidates])

best_idx = int(np.argmin(J))
best_mu = float(candidates[best_idx])
print("sum/nで求めた平均          :", true_mean)
print("J(μ)を最小にするμ(候補探索):", best_mu)
print("両者の差                  :", abs(true_mean - best_mu))

plt.figure()
plt.plot(candidates, J)
plt.axvline(true_mean, color="black", linestyle="--", label="平均(sum/n)")
plt.xlabel("μ")
plt.ylabel("二乗誤差 J(μ)")
plt.title("二乗誤差J(μ)を最小にするμを探す")
plt.legend()

# プログラムファイルと同じディレクトリのパスを取得して保存
script_dir = os.path.dirname(os.path.abspath(__file__))
save_path = os.path.join(script_dir, "ex18_optimize.png")
plt.savefig(save_path)

plt.show()

実行結果の例

ヒント：候補探索（候補をたくさん用意して一番良いものを選ぶ）によるμは、候補の細かさ（今回は400個）の分だけ、本当の最小値（平均）からわずかにずれる。候補をもっと細かくする、あるいは微分を使って正確に解く（dJ/dμ=0を解くとμ=平均になる）ことで、ずれをなくすことができる。

考察ポイント：本演習群を通じて、「ある値を求める」という操作の多くが、実は「ある関数を最小（または最大）にする値を探す」という最適化の問題だったことを確認した。平均は二乗誤差J(μ)=Σ(xᵢ-μ)²を最小にするμであり、画像の距離マップ（演習12）は各画素について点までの距離を最小にする処理であり、周波数の特定（演習17）は内積で求めたスペクトルの強さを最大にする周波数を探す処理であった。内積（演習1〜4）と最適化（本演習）という2つの道具を組み合わせることで、日時・位置・画像・音声という異なる種類のデータを、同じ枠組みで扱えることが分かる。