sobasobasobaの備忘録

numpy で2変数の回帰分析

programming python numpy matplotlib

polyfit.png:

np.polyfit(x, y, n)
- n次式で2変数の回帰分析

np.polyval(p, t):
- pで表される多項式に t を代入し、値を計算して返す
- p[0]*t**(N-1) + p[1]*t**(N-2) + ... + p[N-2]*t + p[N-1]

#!/usr/bin/env python
#coding: utf-8

import matplotlib.pyplot as plt
import numpy as np


def read():
    labels = (("car_price", "gas_price", "weight", "co2", "power", "duration"))
    
    ret = []
    ret.append((8895,33,2560,97,113,12))
    ret.append((7402,33,2345,114,90,8))
    ret.append((6319,37,1845,81,63,12))
    ret.append((6635,32,2260,91,92,13))
    ret.append((6599,32,2440,113,103,13))
    ret.append((8672,26,2285,97,82,12))
    ret.append((7399,33,2275,97,90,13))
    ret.append((7254,28,2350,98,74,8))
    ret.append((9599,25,2295,109,90,13))
    ret.append((8748,29,2390,97,102,13))
    ret.append((6488,35,2075,89,78,13))
    ret.append((9995,26,2330,109,100,9))
    ret.append((11545,20,3320,305,170,8))
    ret.append((9745,27,2885,153,100,8))
    ret.append((12164,19,3310,302,225,8))
    ret.append((11470,30,2695,133,110,9))
    ret.append((9410,33,2170,97,108,13))
    ret.append((13945,27,2710,125,140,13))
    ret.append((13249,24,2775,146,140,10))
    ret.append((10565,23,2640,151,110,8 ))
    ret.append((10320,26,2655,133,95,9  ))
    ret.append((10945,25,3065,181,141,12))
    ret.append((9483,24,2750,141,98,9   ))
    ret.append((12145,26,2920,132,125,13))
    ret.append((12459,24,2780,133,110,12))
    ret.append((10989,25,2745,122,102,13))
    ret.append((17879,21,3110,181,142,12))
    ret.append((11650,21,2920,146,138,13))
    ret.append((9995,23,2645,151,110,9  ))
    ret.append((11499,23,2935,135,130,13))
    ret.append((11588,27,2920,122,115,13))
    ret.append((18450,23,2985,141,114,10))
    ret.append((24760,20,3265,163,160,13))
    ret.append((13150,21,2880,151,110,10))
    ret.append((12495,22,2975,153,150,9 ))
    ret.append((16342,22,3450,202,147,10))
    ret.append((15350,22,3145,180,150,9 ))
    ret.append((13195,22,3190,182,140,10))
    ret.append((14980,23,3610,232,140,8 ))
    ret.append((23300,21,3480,180,158,13))
    ret.append((17899,22,3200,180,160,13))
    ret.append((13150,21,2765,151,110,9 ))
    ret.append((21498,23,3480,180,190,10))
    ret.append((16145,23,3325,231,165,10))
    ret.append((14525,18,3855,305,170,8 ))
    ret.append((17257,20,3850,302,150,10))
    ret.append((15395,18,3735,202,150,10))
    ret.append((12267,18,3665,182,145,10))
    ret.append((14944,19,3735,181,150,13))
    
    return labels, np.array(ret)


def work(x, y, xLabel, yLabel):
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)

    poly = np.polyfit(x, y, 1) # 1次式で回帰分析
    yy = np.polyval(poly, x)   # poly[0]: 傾き,  poly[1]: 切片
    ax.plot(x,yy)

    poly = np.polyfit(x, y, 2) # 2次式で回帰分析
    xx = np.arange(x.min(),x.max(),1000)
    yy = np.polyval(poly,xx)   # poly[0] * t^2 + poly[1] * t^1 + poly[0]
    ax.plot(xx,yy)
    
    ax.plot(x,y,linestyle="None",marker="o")
    ax.set_xlabel(xLabel)
    ax.set_ylabel(yLabel)

    fig.savefig("polyfit.png",bbox_inches="tight")


if __name__ == "__main__":
    labels, X = read()
    work(X[:,0], X[:,4], labels[0], labels[4])

参考:
- http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html
- http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyval.html

numpy で3変数の回帰分析

programming python numpy matplotlib

polyfit.png:

np.polyfit(x, zip(y,z), n)
- n次式で3変数の回帰分析

3Dデータのプロット
1. ax = fig.add_subplot(1, 1, 1, projection="3d", azim=<角度>)
2. ax.plot(x,y,z)
  - azim に値を代入して視点を変える

#!/usr/bin/env python
#coding:utf-8

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

def read():
    labels = (("car_price", "gas_price", "weight", "co2", "power", "duration"))
    
    ret = []
    ret.append((8895,33,2560,97,113,12))
    ret.append((7402,33,2345,114,90,8))
    ret.append((6319,37,1845,81,63,12))
    ret.append((6635,32,2260,91,92,13))
    ret.append((6599,32,2440,113,103,13))
    ret.append((8672,26,2285,97,82,12))
    ret.append((7399,33,2275,97,90,13))
    ret.append((7254,28,2350,98,74,8))
    ret.append((9599,25,2295,109,90,13))
    ret.append((8748,29,2390,97,102,13))
    ret.append((6488,35,2075,89,78,13))
    ret.append((9995,26,2330,109,100,9))
    ret.append((11545,20,3320,305,170,8))
    ret.append((9745,27,2885,153,100,8))
    ret.append((12164,19,3310,302,225,8))
    ret.append((11470,30,2695,133,110,9))
    ret.append((9410,33,2170,97,108,13))
    ret.append((13945,27,2710,125,140,13))
    ret.append((13249,24,2775,146,140,10))
    ret.append((10565,23,2640,151,110,8 ))
    ret.append((10320,26,2655,133,95,9  ))
    ret.append((10945,25,3065,181,141,12))
    ret.append((9483,24,2750,141,98,9   ))
    ret.append((12145,26,2920,132,125,13))
    ret.append((12459,24,2780,133,110,12))
    ret.append((10989,25,2745,122,102,13))
    ret.append((17879,21,3110,181,142,12))
    ret.append((11650,21,2920,146,138,13))
    ret.append((9995,23,2645,151,110,9  ))
    ret.append((11499,23,2935,135,130,13))
    ret.append((11588,27,2920,122,115,13))
    ret.append((18450,23,2985,141,114,10))
    ret.append((24760,20,3265,163,160,13))
    ret.append((13150,21,2880,151,110,10))
    ret.append((12495,22,2975,153,150,9 ))
    ret.append((16342,22,3450,202,147,10))
    ret.append((15350,22,3145,180,150,9 ))
    ret.append((13195,22,3190,182,140,10))
    ret.append((14980,23,3610,232,140,8 ))
    ret.append((23300,21,3480,180,158,13))
    ret.append((17899,22,3200,180,160,13))
    ret.append((13150,21,2765,151,110,9 ))
    ret.append((21498,23,3480,180,190,10))
    ret.append((16145,23,3325,231,165,10))
    ret.append((14525,18,3855,305,170,8 ))
    ret.append((17257,20,3850,302,150,10))
    ret.append((15395,18,3735,202,150,10))
    ret.append((12267,18,3665,182,145,10))
    ret.append((14944,19,3735,181,150,13))
    
    return labels, np.array(ret)


def work(x, y, z, xLabel, yLabel, zLabel):
    fig = plt.figure()
    ax = fig.add_subplot(2,2,1, projection="3d")
    bx = fig.add_subplot(2,2,2, projection="3d", azim = 30)
    cx = fig.add_subplot(2,2,3, projection="3d", azim = 120)
    dx = fig.add_subplot(2,2,4, projection="3d", azim = 210)
    
    poly = np.polyfit(x, zip(y,z), 1) # 1次式で回帰分析
    yy = np.polyval(poly[:,0], x)   # poly[0][0]: y傾き,  poly[1][0]: y切片
    zz = np.polyval(poly[:,1], x)   # poly[0][1]: z傾き,  poly[1][1]: z切片
    ax.plot(x,yy,zz)
    bx.plot(x,yy,zz)
    cx.plot(x,yy,zz)
    dx.plot(x,yy,zz)

    poly = np.polyfit(x, zip(y,z), 2) # 2次式で回帰分析
    xx = np.arange(x.min(),x.max(),1000)
    yy = np.polyval(poly[:,0],xx)
    zz = np.polyval(poly[:,1],xx)
    ax.plot(xx,yy,zz)
    bx.plot(xx,yy,zz)
    cx.plot(xx,yy,zz)
    dx.plot(xx,yy,zz)

    ax.plot(x, y, z, linestyle="None",marker="o")
    ax.set_xlabel(xLabel)
    ax.set_ylabel(yLabel)
    ax.set_zlabel(zLabel)
    ax.set_title("default (azim = -60)")

    bx.plot(x, y, z, linestyle="None",marker="o")
    bx.set_xlabel(xLabel)
    bx.set_ylabel(yLabel)
    bx.set_zlabel(zLabel)
    bx.set_title("azim = 30")
    
    cx.plot(x, y, z, linestyle="None",marker="o")
    cx.set_xlabel(xLabel)
    cx.set_ylabel(yLabel)
    cx.set_zlabel(zLabel)
    cx.set_title("azim = 120")
    
    dx.plot(x, y, z, linestyle="None",marker="o")
    dx.set_xlabel(xLabel)
    dx.set_ylabel(yLabel)
    dx.set_zlabel(zLabel)
    dx.set_title("azim = 210")

    fig.tight_layout()
    fig.savefig("polyfit.png",bbox_inches="tight")
    

if __name__ == "__main__":
    labels, X = read()

    work(X[:,0], X[:,2], X[:,4], labels[0], labels[2], labels[4])

参考:
- http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html
- http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyval.html

matplotlibでサーファイスプロットする

programming python matplotlib numpy

schwefel.png:

surf = ax.plot_surface(x, y, z)
- 三次元プロット
fig.colorbar(surf)
- カラーバーの表示

plot_surface関数のキーワード(一部)
- cmap: カラーマップの指定 (hot, gray, coolwarm など)
- cstride: x方向の色のプロットのステップ数(デフォルト値: 10)
- rstride: y方向の色のプロットのステップ数(デフォルト値: 10)

colorbar関数のキーワード
- shrink: カラーバーを縦の長さを何倍にするか (デフォルト値: 1.0)
- aspect: カラーバーの横の長さ(値が大きいほど細くなる) (デフォルト値: 20)

#!/usr/bin/env python
#coding:utf-8


from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np


def schwefel(x, y):
    return 418.9829 * 2 - (x * np.sin(np.sqrt(abs(x))) + \
                           y * np.sin(np.sqrt(abs(y))))


def work():
    N = 100
    MIN_X = -500
    MAX_X =  500

    x = np.linspace(MIN_X, MAX_X, N)
    y = np.linspace(MIN_X, MAX_X, N)
    x, y = np.meshgrid(x, y)
    z = schwefel(x, y)
    
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection="3d")

    surf = ax.plot_surface(x, y, z, cstride=1, rstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
    
    fig.colorbar(surf, shrink=0.5, aspect=10)

    ax.set_title("Schwefel")

    fig.savefig("schwefel.png", bbox_inches="tight")
    

if __name__ == "__main__":
    work()

参考:
- http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#mpl_toolkits.mplot3d.Axes3D.plot_surface
- http://www-optima.amp.i.kyoto-u.ac.jp/member/student/hedar/Hedar_files/TestGO_files/Page2530.htm

matplotlibで等高線図を書く

programming python matplotlib numpy

schwefel.png:

ax.contour(x, y, z, zdir=, offset=<平面のどの位置に描画するか>)

contour関数のオプション(一部)
- zdir: どの平面に投射するか(x, y, zのいずれか)
- offset: 平面のどの位置に描画するか
- cmap: カラーマップ

#!/usr/bin/env python
#coding:utf-8

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np


def schwefel(x, y):
    return 418.9829 * 2 - (x * np.sin(np.sqrt(abs(x))) + \
                           y * np.sin(np.sqrt(abs(y))))


def work():
    N = 100
    MIN_X = -500
    MAX_X =  500

    x = np.linspace(MIN_X, MAX_X, N)
    y = np.linspace(MIN_X, MAX_X, N)
    x, y = np.meshgrid(x, y)
    z = schwefel(x, y)
    
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection="3d")

    ax.contour(x, y, z, zdir='z', offset=0, cmap=cm.coolwarm)
    ax.contour(x, y, z, zdir='x', offset=MIN_X, cmap=cm.coolwarm)
    ax.contour(x, y, z, zdir='y', offset=MAX_X, cmap=cm.coolwarm)

    ax.set_title("Schwefel")

    fig.savefig("schwefel.png", bbox_inches="tight")
    

if __name__ == "__main__":
    work()

参考:
- http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#contour-plots
- http://www-optima.amp.i.kyoto-u.ac.jp/member/student/hedar/Hedar_files/TestGO_files/Page2530.htm

matplotlibで等高線を書くその2

programming python matplotlib numpy

schwefel.png:

CS = ax.contour(x, y, z, <等高線の領域数>)
- 等高線図をプロット
ax.clabel(CS)
- 各等高線に値を記入

contour関数のオプション (一部)
- cmap: カラーマップ

clabel関数のオプション (一部)
- fontsize: フォントの大きさ
- inline: True なら、値の周りの等高線を消す (デフォルトはTrue)
- colors: 値の色を設定 (何も設定しない場合、等高線と同じ色となる)

#!/usr/bin/env python
#coding:utf-8


from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np


def schwefel(x, y):
    return 418.9829 * 2 - (x * np.sin(np.sqrt(abs(x))) + \
                           y * np.sin(np.sqrt(abs(y))))


def work():
    N = 100
    MIN_X = -500
    MAX_X =  500

    x = np.linspace(MIN_X, MAX_X, N)
    y = np.linspace(MIN_X, MAX_X, N)
    x, y = np.meshgrid(x, y)
    z = schwefel(x, y)
    
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)

    CS = ax.contour(x, y, z, 5, cmap=cm.coolwarm)
    ax.clabel(CS, fontsize=9, inline=False, colors="black")

    ax.set_title("Schwefel")

    fig.savefig("schwefel.png", bbox_inches="tight")
    

if __name__ == "__main__":
    work()

参考:
- http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#contour-plots
- http://www-optima.amp.i.kyoto-u.ac.jp/member/student/hedar/Hedar_files/TestGO_files/Page2530.htm

Graphviz インストール手順

Graphvis

Graphvis を使えば、データ構造としてのグラフの可視化を簡単に行える。まずはインストールするための手順をメモ。(※自分の環境: Windows7, Cygwin (8.15-1))

http://www.graphviz.org/Download_windows.php より、graphviz-2.34.msi をダウンロード
graphvis-2.34.msi をダブルクリック。インストール先は C:\cygwin\usr\local\graphviz を指定。
PATH を通す。~/.bash_profile に以下の文を追加。

PATH="/usr/local/graphviz/bin:${PATH}"

.bash_profile を読み込みなおす。

$ source ~/.bash_profile

動作確認。以下のコマンドを実行して、バージョンが表示されたらOK。

$ dot -v

参考: http://homepage3.nifty.com/kaku-chan/graphviz/chapter_001.html

Graphvis 使用例

Graphvis

以下の内容で、sample.dot を作成

digraph {
    rankdir=LR;
    
    0 [ shape = doublecircle, label = "0 \n generate=5" ];
    1 [ shape = doublecircle, label = "1 \n generate=2" ];
    3 [ shape = rect, label = "3 \n consume=2" ];
    4 [ shape = rect, label = "4 \n consume=1" ];
    5 [ shape = rect, label = "5 \n consume=4" ];
    
    0 -> 0 [ label = "1" ];
    0 -> 1 [ label = "2" ];
    0 -> 2 [ label = "5" ];
    1 -> 0 [ label = "1" ];
    1 -> 2 [ label = "8" ];
    2 -> 3 [ label = "1" ];
    2 -> 4 [ label = "7" ];
    3 -> 5 [ label = "2" ];
    3 -> 6 [ label = "5" ];
    4 -> 2 [ label = "7" ];
    4 -> 3 [ label = "5" ];
    4 -> 5 [ label = "1" ];
    6 -> 0 [ label = "5" ];
}

以下のコマンドを入力

$ dot -Tpng sample.dot -o sample.png

画像 sample.png が出力される

参考:
- 書いたグラフ(サンプルの2番目): http://poj.org/problem?id=1459
- 属性リスト: http://www.graphviz.org/content/attrs
- ノードの形状リスト: http://www.graphviz.org/content/node-shapes

Graphvis 使用例その2

Graphvis

以下の内容で、sample.dot を作成

graph {
    1 -- 6;
    1 -- 12;
    1 -- 15;
    1 -- 16;
    2 -- 3;
    2 -- 8;
    2 -- 10;
    3 -- 4;
    3 -- 6;
    3 -- 9;
    5 -- 6;
    7 -- 15;
    11 -- 16;
    13 -- 15;
    14 -- 15;
}

以下のコマンドを入力

$ dot -Tpng sample.dot -o sample.png

画像 sample.png が出力される

参考:
- 書いたグラフ(サンプルの2番目): http://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=1184
- 属性リスト: http://www.graphviz.org/content/attrs
- ノードの形状リスト: http://www.graphviz.org/content/node-shapes