フィッティングパラメータが2個の場合、カイ自乗は
      　Δχ²=Σ_ijB_ijθ_iθ_j
で表されます。 B=[[1,-1],[-1,4]] について計算した結果が図7.5です。このcode7.3.py では、Berendsen氏オリジナルの contour を用いて等高線上の点を求めています。 点は autoplotp でプロットすればよいとありますがそれがみつかりません。そこでポピュラーな plot を使って右図を得ることができました。氏のコードcode7.3.pyに対する主な変更点は data ではなく xlist, ylist に分けて返すようにしたことです。作画のプログラムを以下にしまします。

1. def contour:
2. ･･････（元々の code 7.3。ただし 修正点がいくつかあります）
3. def fxy(x,y):
4.    v = np.array([[x,y]])
5.    w = np.dot(v, matrx @ v.T)
6.    return w
7. plt.xlim(-5., 5.)
8. plt.ylim(-3., 3.)
9. xycenter= np.array([0.0, 0.0])
10. for i in range(1,6):
11.    z = i
12.    xlist,ylist = contour(fxy,z,xycenter)
13.    plt.plot(xlist,ylist,color='blue')
14. font1 = {'family':'serif','color':'blue','size':11}
15. font2 = {'family':'serif','color':'darkred','size':12}
16. plt.title(r"$\Delta \chi^2$ contours (Fig.7.5)", fontdict = font1)
17. plt.xlabel(r"$\Delta \theta_1$", fontdict = font2)
18. plt.ylabel(r"$\Delta \theta_2$", fontdict = font2)
19. plt.grid()
20. plt.show()

＜コメント＞

fxy(x,y)=v B v^T　において v は [x,y] という横ベクトル、v^Tは縦ベクトルです。 B行列を指定し、展開すると x²-2xy+4y² ですが、左では練習のためにマトリックス計算で求めています。

元々の code 7.3 への修正点は以下のとおりです。

• import numpy as np を追加
• import matplotlib.pyplot as plt を追加
• return data を削除
• return xlist,ylist を追加

こちらを使うと描画までやってくれます。ベクトルvに対してv  B v^T を計算しています。

1. import numpy as np
2. import matplotlib.pyplot as plt
3. matrx = np.array([[1.0, -1.0],[-1.0, 4.0]])/3.0
4. nxpt = 201; nypt = 121
5. x = np.linspace(-5., 5.,nxpt)
6. y = np.linspace(-3., 3.,nypt)
7. z = np.zeros((nypt,nxpt))
8. xx, yy = np.meshgrid(x,y)
9. def findz(x,y):
10.    for i in range(nypt):
11.      yi = y[i,0]
12.      for j in range(nxpt):
13.        xj = x[0,j]
14.        v = np.array([xj,yi])
15.        w = matrx @ v.T
16.        z[i,j] = np.dot(v, w)
17.    return z
18. zz = findz(xx, yy)
19. # zz = (xx**2 - 2.*xx*yy + 4.*yy**2)/3.0
20. xycenter= np.array([0.0, 0.0])
21. figure = plt.contourf(xx,yy,zz,levels=[1.,2.,3.,4.,5.],\
22.    colors=['#FF6666','#FFB266','#FFFF66','#B2FF66', \
23.    '#66FF66','#66FFFF'])
24. font1 = {'family':'serif','color':'blue','size':11}
25. font2 = {'family':'serif','color':'darkred','size':12}
26. plt.title(r"$\Delta \chi^2$ contours (Fig.7.5)", fontdict = font1)
27. plt.xlabel(r"$\Delta \theta_1$", fontdict = font2)
28. plt.ylabel(r"$\Delta \theta_2$", fontdict = font2)
29. plt.grid()
30. plt.show()

＜コメント＞

19行目の＃を18行目に移せば展開計算になります。