Matplotlib その 1

ここでは,強力なプロットツール matplotlib を使って2次元プロットについて学びます.ダウンロードおよびマニュアルは,こちらへどうぞ → http://matplotlib.sourceforge.net/

matplotlib.pyplot に親しむ

matplotlib では pylab というモジュールを使った簡便な(Matlabライクな)プロット方法が用意されていますが,オブジェクト指向的な記法が推奨されています.こちらを使うには,matplotlib.pyplot をインポートしてやります.

以下では,実用上必要になりそうな手法をいくつか具体例を通して学んでいきます.

関数のプロット

まずは普通の関数のプロットから.矩形波とその三角近似をプロットします.細かい方法はスクリプト内のコメントを参照してください.

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

# plotting functions

import numpy as np
import matplotlib.pyplot as plt

# 矩形波
f = lambda x: np.where(np.abs(np.floor(x/np.pi)%2) < 0.001, 1.0, -1.0)
# 三角近似
g = lambda x: (4.0/np.pi) * ( np.sin(x) + np.sin(3*x)/3.0 + np.sin(5*x)/5.0 )

x = np.linspace(0.0, 6.0 * np.pi, 1000)

# プロットが描かれる土台のウインドウを作る
fig = plt.figure()

# fig の中にプロット領域を作る
ax = fig.add_subplot(111)

# プロット領域の範囲を決める
ax.set_xlim(0.0, 6.05*np.pi)
ax.set_ylim(-1.8, 1.8)

# ax 上にプロットする
fplot = ax.plot(x, f(x), 'k-')
gplot = ax.plot(x, g(x), 'r-')

# 凡例を入れる
ax.legend((fplot,gplot),('Rectangular Wave', 'Fourier Approximation'))

# 軸の名前を入れる
ax.set_xlabel('time')
ax.set_ylabel('amplitude')

# 横軸のティックを決める
ax.set_xticks(np.arange(7)*np.pi)
ax.set_xticklabels(['$0$','$pi$', '$2pi$','$3pi$','$4pi$','$5pi$','$6pi$'],
                    size='15')

# タイトルを挿入
ax.set_title('simpleplot.py')

# 保存と表示
fig.savefig('simpleplot.png', dpi=96)
fig.show()

実行すると次のような図が得られます.

等高線プロット

次に,Cobb-Douglas型生産関数の等生産量曲線をプロットしましょう.ここでは,関数 $latex f$ がCobb-Douglas型とは,

 f(x, y) = A x^a y^b

が成り立つ$latex A, a, b > 0$ が存在することとします.2通りのパラメータについて等高線プロットを行います:

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

# plotting isoquant curves

import numpy as np
import matplotlib.pyplot as plt

def cobb_douglas(A, a, b):
    "Cobb-Douglas Class of Functions"
    func = lambda x, y: A * (x ** a) * (y ** b)
    return func

#
# plot
# 

A1, a1, b1 = 1.0, 0.4, 0.6
A2, a2, b2 = 1.5, 0.6, 0.3

title1 = 'A = %.1f, a = %.1f, b = %.1f' % (A1, a1, b1)
title2 = 'A = %.1f, a = %.1f, b = %.1f' % (A2, a2, b2)

func1 = cobb_douglas(A1, a1, b1)
func2 = cobb_douglas(A2, a2, b2)

x = np.linspace(0.0, 10.0, 250)
y = np.linspace(0.0, 10.0, 250)

X, Y = np.meshgrid(x,y)
Z1 = func1(X, Y)
Z2 = func2(X, Y)

fig = plt.figure()

ax1 = fig.add_subplot(121, aspect='equal')
IQ = ax1.contour(X, Y, Z1)
ax1.set_title(title1)

ax2 = fig.add_subplot(122, aspect='equal')
IQ_bw = ax2.contour(X, Y, Z2, colors = 'k')
ax2.set_title(title2)

fig.savefig('isoquant.png', dpi = 96)
fig.show()

出力される図は次のような感じ

ax1 = fig.add_subplot(121, aspect='equal')
IQ = ax1.contour(X, Y, Z1)

この部分は,1行2列のプロットの表をつくって,1番目(左側のプロット領域)に,Z1の等高線をプロットするということを意味しています.縦横の比はデフォルトでは横長ですが,ここでは1対1にしています(aspect=’equal’)

ax2 = fig.add_subplot(122, aspect='equal')
IQ_bw = ax2.contour(X, Y, Z2, colors = 'k')

ここでは,1行2列のプロットの表の,2番目(右側のプロット領域)に,Z2の等高線をプロットしています.今回は,色を黒にしました(colors=’k')

領域を表示する

ある条件を満たすパラメータの範囲を図示したい場合はしょっちゅう起こります.ここでは,問題を簡単にするため,
パラメータ x, y が

 { 2(x- 1/2)^2 < y < 1/10 + x/8 }

を満たす領域を図示したいとしましょう.結果は次の図になります.

この作図は fill_between() メソッドを使えばよいです.

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

# showing region

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0.0, 1.0, 100)

f_upper = lambda x: 0.1 + 0.8 * x
f_lower = lambda x: 2*(x - 0.5)**2

fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal')
upper, lower = ax.plot(x,f_upper(x),'k--', x, f_lower(x),'k-.')

ax.legend((upper, lower), ('upper bound','lower bound'), loc = 2)
ax.fill_between(x, f_upper(x), f_lower(x), where = f_upper(x) > f_lower(x), color='#cccccc')

fig.savefig('region.png', dpi=96)
fig.show()

fill_between() メソッドと,where オプションの使い方に注意してください.

One response so far

  • Kenji Sato says:

    スクリプトの最終行 fig.show() は,IPython で RUN するとうまくいきますが,端末で実行するとうまく行かないかもしれません.そういう場合は,この部分は,plt.show() に書き換えてください

Leave a Reply