00_PIP.md

パターン情報処理

2017年度第13回&第14回の宿題（の一部）

http://www-tlab.math.ryukoku.ac.jp/wiki/index.php?PIP%2F2017%2Fhw1234

• gendat.py
• leastsquares.py

gendat.py

# -*- coding: utf-8 -*-

from __future__ import division
from __future__ import print_function

import numpy as np
import os

def trueFunction(x):

    return 0.2 * x + 0.5


def gendat(n, seed = 0, nsig = 0.0):

    np.random.seed(seed)
    x = np.linspace(0.0, 10.0, num = n)
    y = trueFunction(x) + nsig * np.random.randn(n)

    return np.vstack((x, y)).T


if __name__ == '__main__':

    N = 21

    login = os.getlogin()
    try:
        idnum = int(login[1:])
    except:
        idnum = 0
        
    # N 個のデータを生成
    X = gendat(N, seed = idnum, nsig = 0.2)

    # テキストファイルに保存
    np.savetxt('hoge.txt', X, fmt = '%.3f')
    

leastsquares.py

# -*- coding: utf-8 -*-

from __future__ import division
from __future__ import print_function

import numpy as np


### 最小二乗法（正規方程式を解いてパラメータを求める）
#
# x と y は 同じ長さのベクトル
#
def solve(x, y):

    N = len(x)  # N: データ数
    X = np.vstack((np.ones(N), x))  # X: 2 x N 行列
    
    A = np.dot(X, X.T)   # A = XX^T
    b = np.dot(X, y)     # b = XY^T

    # 連立方程式 Aw = b の解 w を求める
    w = np.linalg.solve(A, b)

    return w


### 最小二乗法による直線あてはめを実行し，結果を表示
#
# x と y は 同じ長さのベクトル
#
def estimate(x, y):

    N = x.shape[0]
    print('# N =', N)
    print('# x =', x)
    print('# y =', y)

    w = solve(x, y)
    print('# 推定されたパラメータ: ', w)

    b, a = w  # w[0] is b and w[1] is a
    print('# 推定された直線の式: y =', a, '*x +', b)
    

if __name__ == '__main__':

    ### 講義資料 Q1 ###
    x = np.array([0, 4, 8, 20])
    y = np.array([0, 1, -1, -2])
    estimate(x, y)
    print()
    
    ### ファイルからデータを読み込み ###
    fn = 'hoge.txt'
    try:
        data = np.loadtxt(fn)
    except FileNotFoundError:
        print(fn + ' is not found')
        quit()
        
    x = data[:, 0]  # 1列目が x 
    y = data[:, 1]  # 2列目が y
    estimate(x, y)
    print()