[Python]Matplotlibで2次元ヒストグラムを作成する方法

２変数のヒストグラムを描画する場合、「３次元のグラフ」にするか、「２次元で頻度を色または濃度に割り当てる」かどちらかになる。
３次元グラフは頻度方向の変化がわかりやすい反面、隠れてしまう部分が存在するために全体の分布がわかりづらくなることもある。
一方２次元で頻度を色や濃度に割り当てた場合は、頻度方向の微妙な差はわかりづらいが、全体の分布がどうなっているか把握しやすい。

## データ
２つの２次元正規分布をデータとして使う。分布は下記の通り。

```python
import numpy as np
x, y = np.vstack((np.random.multivariate_normal([0, 0], [[10.0, 0],[0,20]], 5000) 
                 ,np.random.multivariate_normal([0,15], [[10.0, 0],[0, 5]], 5000))).T
```
![hist2d_01.png](https://qiita-image-store.s3.amazonaws.com/0/100523/33dbb7d5-aa73-9e22-a3cb-0b823645dbcb.png)

## 2D histogram 
２次元ヒストグラムははmatplotlibのhist2dを使用する。
ヒストグラムの頻度は戻り値として得られる。
戻り値は、`counts, xedges, yedges, Image`となっている。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hist2d(x,y, bins=40, cmap=cm.jet)
ax.set_title('1st graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.colorbar(H[3],ax=ax)
plt.show()
```
![hist2d_02.png](https://qiita-image-store.s3.amazonaws.com/0/100523/fb796181-d3be-a556-a331-5db6a08fb2c0.png)

## ビン数の指定
ビン数のしては、パラメータbinsで行う。スカラーで指定すると、縦横とも同じビン数になる。
別々に指定したい場合は、例のようにする。
１次元のヒストグラム同様にedgeを指定することもできる。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hist2d(x,y, bins=[40,10], cmap=cm.jet)
ax.set_title('2nd graph')
ax.set_xlabel('x')
ax.set_ylabel('y')

fig.colorbar(H[3],ax=ax)
plt.show()
```
![hist2d_03.png](https://qiita-image-store.s3.amazonaws.com/0/100523/e683b115-3460-64c2-f137-e862f5b1af21.png)

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hist2d(x,y, bins=[np.linspace(-30,30,61),np.linspace(-30,30,61)], cmap=cm.jet)
ax.set_title('3rd graph')
ax.set_xlabel('x')
ax.set_ylabel('y')

fig.colorbar(H[3],ax=ax)
plt.show()
```
![hist2d_04.png](https://qiita-image-store.s3.amazonaws.com/0/100523/fe9629a1-8349-0c12-990f-c9aace680bbc.png)

## 正規化
ヒストグラムを正規化したいときは、パラメータnormedをTrueにする。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hist2d(x,y, bins=[np.linspace(-30,30,61),np.linspace(-30,30,61)], normed=True, cmap=cm.jet)
ax.set_title('4th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')

fig.colorbar(H[3],ax=ax)
plt.show()
```
![hist2d_05.png](https://qiita-image-store.s3.amazonaws.com/0/100523/8062ade5-9cf9-fed8-c35f-8d2bde3d0f6e.png)

## カラーマップを変更する
カラーマップを変更する場合は、パラメータcmapに例のように指定する。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hist2d(x,y, bins=[np.linspace(-30,30,61),np.linspace(-30,30,61)], normed=True, cmap=cm.gray)
ax.set_title('5th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.colorbar(H[3],ax=ax)
plt.show()
```
![hist2d_06.png](https://qiita-image-store.s3.amazonaws.com/0/100523/7e3fdb49-99f6-8756-cec3-2a9eb161b323.png)

## カラーマップの範囲を指定する
複数のヒストグラムを比較する場合にカラーマップの範囲を指定したいことがある。
この場合は、例のようにset_climを使用する。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hist2d(x,y, bins=[np.linspace(-30,30,61),np.linspace(-30,30,61)], normed=True, cmap=cm.jet)
ax.set_title('6th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
H[3].set_clim(0,0.05)
fig.colorbar(H[3],ax=ax)
plt.show()
```
![hist2d_07.png](https://qiita-image-store.s3.amazonaws.com/0/100523/7d6cf3ad-9405-9693-86b0-5cfdcfefd15d.png)

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hist2d(x,y, bins=[np.linspace(-30,30,61),np.linspace(-30,30,61)], normed=True, cmap=cm.jet)
ax.set_title('7th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
H[3].set_clim(0,0.01)
fig.colorbar(H[3],ax=ax)
plt.show()
```
![hist2d_08.png](https://qiita-image-store.s3.amazonaws.com/0/100523/bd6d3a16-3cbb-3e8f-c8cb-2ca1f351bd8e.png)

## Logスケールにする
ヒストグラムをLogスケールにする場合は、matplotlib.colors.LogNorm()を使用する。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hist2d(x,y, bins=[np.linspace(-30,30,61),np.linspace(-30,30,61)], norm=matplotlib.colors.LogNorm(), cmap=cm.jet)
ax.set_title('8th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')

fig.colorbar(H[3],ax=ax)
plt.show()
```
![hist2d_09.png](https://qiita-image-store.s3.amazonaws.com/0/100523/8dd4f188-bfb4-561a-660e-9657b27e041f.png)

## ヒストグラムの上に等高線（Contour)を重ねる

等高線はcontourを使って書く。
このとき下記の点に注意した
1. histogramと水平、垂直位置が一致するように、extentで最大値、最小値を指定する
2. countsは垂直方向にx軸のデータが、水平方向にy軸のデータが入っているので転置(transpose)する
3. LogNormを使って対数表示にする。これはなくてもいいけど無いと線がほとんど見えない。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

counts, xedges, yedges, Image= ax.hist2d(x,y, bins=[np.linspace(-30,30,61),np.linspace(-30,30,61)], norm=matplotlib.colors.LogNorm(), cmap=cm.jet)
ax.contour(counts.transpose(),extent=[xedges.min(),xedges.max(),yedges.min(),yedges.max()])
ax.set_title('8th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.colorbar(Image,ax=ax)
plt.show()
```
![hist2d_10.png](https://qiita-image-store.s3.amazonaws.com/0/100523/65de57a2-ccc2-eec1-5b89-c9e804df1581.png)


## ヒストグラムのビンの形を六角形にする
ビンの形を六角形にしたいとき。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = ax.hexbin(x,y, gridsize=20, extent=[-30, 30, -30, 30], cmap=cm.jet)
ax.set_title('8th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.colorbar(H,ax=ax)
plt.show()
```
![hist2d_11.png](https://qiita-image-store.s3.amazonaws.com/0/100523/8275df98-94f7-2e34-e9f2-84b59b50ca03.png)

## ヒストグラムを滑らかにする(kernel density estimation)

ヒストグラムはビンの間隔が広い時は特にBlockyになる。
これを滑らかにつなげる方法にカーネル密度推定というものがある。
これはランダム変数の標本分布から確率密度関数を推定するためのものである。

kernel density estimationは、scipyやscikit-learnなどに入っている。
scipyのkernel density estimationはband幅を標準偏差で正規化して指定する。
これによりデータの分布が変わってもそれなりに滑らかに推定される。
したがって、パラメータbw_method=1.0とした場合は、value.std(ddof=1)がband幅として使われる。(valueはデータ)
ここでddofはdelta degrees of freedom で標準偏差を計算するときにN-ddofで割られる。

`kernel = gaussian_kde(value)`でkernelはgaussian_kdeオブジェクトなので、実際に値を得る場合は、
`kernel(positions)`のようにx,y座標の位置を渡す。

実際にはmgridでメッシュを作って、ravelでx,yを別々に一次元にしてvstackでくっつけてからkernelに渡している。

`kernel(positions)`で戻ってきた値は1次元になっているのでreshapeで２次元にしている。

最終的にはcontourfでグラフとして表示する。

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
import matplotlib.cm as cm

xx,yy = np.mgrid[-30:30:1,-30:30:1]
positions = np.vstack([xx.ravel(),yy.ravel()])
 
value = np.vstack([x,y])

kernel = gaussian_kde(value)

f = np.reshape(kernel(positions).T, xx.shape)

fig = plt.figure()
ax = fig.add_subplot(111)

ax.contourf(xx,yy,f, cmap=cm.jet)
ax.set_title('11th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()
```
![hist2d_12.png](https://qiita-image-store.s3.amazonaws.com/0/100523/7bf8af8d-f7cf-32f0-51e3-42b779038e47.png)

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
import matplotlib.cm as cm

xx,yy = np.mgrid[-30:30:1,-30:30:1]
positions = np.vstack([xx.ravel(),yy.ravel()])
 
value = np.vstack([x,y])

kernel = gaussian_kde(value, bw_method=0.5)

f = np.reshape(kernel(positions).T, xx.shape)

fig = plt.figure()
ax = fig.add_subplot(111)

ax.contourf(xx,yy,f, cmap=cm.jet)
ax.set_title('12th graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()
```
![hist2d_13.png](https://qiita-image-store.s3.amazonaws.com/0/100523/331fa06c-a993-65a3-ad35-f4a6646c5d8d.png)

## numpy のhistogram2dを使った場合

numpyのhistogram2dを使ってヒストグラムを作ることができる。
この場合は、頻度とx,yのedgeが得られるので、自分でグラフにする必要がある。
ここでは表示にimshowを使用した。
この時ヒストグラムのデータはmatplotlibのhist2dと同じように垂直方向にｘ軸方向のデータ、水平方向にｙ軸方向のデータが入っているので、
転置をして更に始点をlower leftにしている。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

fig = plt.figure()
ax = fig.add_subplot(111)

H = np.histogram2d(x,y, bins=[np.linspace(-30,30,61),np.linspace(-30,30,61)])
im = ax.imshow(H[0].T, interpolation='nearest', origin='lower', extent=[-30,30,-30,30], cmap=cm.jet)
ax.set_title('13st graph')
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.colorbar(im, ax=ax)
plt.show()
```
![hist2d_14.png](https://qiita-image-store.s3.amazonaws.com/0/100523/83cd76a6-0bc8-f558-75b5-6166fdc306f1.png)