富士山が美しすぎたので、フィッティングをしてしまいました...

こんにちは、ほけきよです。あけましておめでとうございます。

新年ということで、おめでたいものを。おめでたいものといえば、富士山ですね！昨年、静岡に旅行したときに、山中湖周辺で撮った一枚がこちら

ううむ、完璧に近い形。外国人が"Oh, Fujiyama!!!"とエキサイテッドする理由もわかります。

かくいう私も、年始に写真フォルダを眺めていると、たまたまこの写真を見つけて、エキサイテッドしたわけです

うおお、フィッティングしよう

というわけで、今日は富士山をフィッティングしてみた話です。

手順
エッジ抽出
正規分布でフィッティング
なぜフィッティングできなかったのか。
外れ値に強いt分布
コード

手順

はじめに、簡単に手順をまとめておきます

写真から富士山の輪郭を抽出(エッジ抽出)

富士山の輪郭を数値化(分布)

富士山分布からデータを抽出(サンプリング)

それっぽい分布に当てはめる(最尤推定)

という感じ。データの変換の仕方とか細々としたところはコードを御覧ください。

エッジ抽出

画像の扱いはやっぱりOpenCVが便利なので、OpenCVを使います。流石に、エッジ抽出くらいだと関数が用意されているので、モノクロ化→輪郭抽出と進めていきます。

中間生成物と、取れた輪郭がこちら

めちゃめちゃキレイにエッジが取れたので、富士山の分布を作るのが簡単です。 画像の最上部の白い部分だけ取り出すという処理を書くだけで、こんな感じで輪郭が取れちゃいました。

正規分布でフィッティング

これだけきれいな形なので、まずは正規分布でフィッティングしてみましょう。

$N(x,\mu,\sigma)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\mathrm{exp}\left({-\frac{(x-\mu)^{2}}{2\sigma^{2}}}\right)$

正規分布の最尤推定は解析的に求められて、 $n$ 個のデータ点 ${x_1,x_2,...,x_n}$ があるとすると、平均と分散が

$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}x_i$

$\hat{\sigma}=\frac{1}{n}\sum_{i=1}^{n}(x_i-\mu)^{2}$

で表すことができます*1。普段我々が慣れ親しんでいる「平均」と「分散」に値が一致していることを確認してください。

さて、正規分布での平均と分散は山の中心と裾の広がり具合です。実際に富士山に重ねてみましょう。

おや？？全然フィッティングできていないように感じます。試しに、平均じゃなくて中央値を使ってみました(robust normalというやつ)が、頂点が正しくなったくらいで、裾の拡がり方などはあんまり変わらなかったです。

なぜフィッティングできなかったのか。

ここでもう一度富士山の写真を見てみましょう。特に両端。

裾の広がりが広い(ロングテールっぽい)

右側はなにかに邪魔されている(外れ値っぽい)

ことが見て取れます。 裾野が広く、外れ値に強い...？そんな分布ありましたね。

外れ値に強いt分布

そう、t分布です。

$St(x)=\frac{1}{\sqrt{\nu\pi}} \frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)} \left( 1+\frac{x^{2}}{\nu}\right)^{-\frac{\nu+1}{2}}$

自由度 $\nu$ というパラメータを持っていて、裾の広がり具合などが調整できます。ロバスト性といって、外れ値に強く、 $\nu \rightarrow \infty$ で正規分布と一致するという特徴を持つ分布です。

「そうか、t分布で最尤推定すれば、イケるんじゃない？」

と思ったまでは良いけど、、、案外t分布の最尤推定って見つからないもんですね。なんかライブラリでもあるかな？と軽い気持ちで試したかったのに。

調べてみると下記論文でt分布の最尤推定をやっているみたいです。正規分布みたいに解析的に求まらないので、EMアルゴリズムを使います。

ML ESTIMATION OF THE t DISTRIBUTION USING EM AND ITS EXTENSIONS, ECM AND ECME

多変数だし添字多いしでなかなか辛かったのですが、とりあえず論文をさらっと読んでふわっと実装してみました。

実際に最尤推定してフィッティングしてみると・・・？

おお!!

正規分布でのフィッティングより合っているように見えます！！さすがのロバスト性ですね！工場のデータとか眺めていても思うのですが、外れ値多いし、こういうロバストな推定方法を知っておくのは良いことなんじゃないのかな～と思います。かなり難易度上がるけど。

しかし富士山は美しいですね。実際はすこーしだけ非対称なようなので、そんな分布を持ってきてフィッティングしてみても良いのではないでしょうか？

コード

今回書いたコードです。実際に最尤推定が正しいのかどうか微妙なところなので、ツッコミ待ちです。*2

コードをみる

import cv2
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import gamma, digamma

# 定数定義
CANNY_WINDOW_NAME = "canny"
ORG_FILE_NAME = "fujiyama.JPG"
CANNY_FILE_NAME = "fujiyama_canny.png"
# 元の画像を読み込む
org_img = cv2.imread(ORG_FILE_NAME, cv2.IMREAD_UNCHANGED)
# グレースケールに変換
gray_img = cv2.imread(ORG_FILE_NAME, cv2.IMREAD_GRAYSCALE)
# エッジ抽出
canny_img = cv2.Canny(gray_img, 50, 110)
# ファイルに保存
cv2.imwrite(CANNY_FILE_NAME, canny_img)
# 画像から数値を抽出
ns = []
max_height = canny_img.shape[0]
for i, c in enumerate(canny_img.T):
    flg = True
    ns.append(max_height - np.argmax(c)) #座標反転用(転置すると上下反転していた)
ns = np.array(ns)
offset = ns.min()
fujiyama = ns - offset #山の稜線を出すfujiyamaデータ。

# データ配列用に並べ替える
sample = []
for i,n in enumerate(fujiyama):
    sample_ = [i]*n # 山の高さ=データの数とする
    sample.extend(sample_)
sample = np.array(sample)


class Gaussian(object):
    """
    fit : ハイパーパラメータ(mean, sigma)の調整だけでOK
    predict : 正規分布
    """
    def fit(self, x):
        self.mu = np.array(x).mean()
        self.sigma = np.array(x).std()
    
    def predict(self, x):
        return 1/np.sqrt(2*np.pi*self.sigma**2)*(np.exp(-(x-self.mu)**2/(2*self.sigma**2)))
    
class RobustGaussian(object):
    """
    平均の代わりにmedianを使う
    """
    def fit(self, x):
        self.mu = np.median(np.array(x))
        self.sigma = np.sqrt(np.dot(np.array(x)- self.mu, np.array(x)-self.mu)/len(x))

    def predict(self,x):
        return 1/np.sqrt(2*np.pi*self.sigma**2)*(np.exp(-(x-self.mu)**2/(2*self.sigma**2)))

class StudentsT(object):

    def __init__(self, mu=0, sigma=1, nu=1, iteration=10, solve_nu=False):
        self.mu = mu
        self.sigma = sigma
        self.nu = nu
        self.iteration = iteration
        self.solve_nu = solve_nu
        
    def fit(self, x, N):
        self.n = len(x)
        for _ in range(N):
            # E step
            self.w = (self.nu+1)/(self.nu+(x-self.mu)**2/self.sigma)
            # M step
            self.mu = np.sum(self.w*x)/np.sum(self.w)
            self.sigma = 1/self.n*(np.sum(self.w*(x-self.mu)*(x-self.mu)))
            # nuを求める(二分法)
            if self.solve_nu:
                _a = 1
                _b = 3
                for i in range(self.iteration):
                    _x = (_a+_b)/2.0 # 中点を出す  
                    if self._expect_nu(_a)*self._expect_nu(_x) < 0:
                        _b = _x
                    if self._expect_nu(_b)*self._expect_nu(_x) < 0:
                        _a = _x
                self.nu = _x
    
    def _expect_nu(self, nu):
        """
        最尤推定でnuを求めるときの方程式
        """
        x = -digamma(nu/2)+np.log(nu/2)
        a1 = 1/self.n*np.sum(np.log(self.w)-self.w)+1
        a2 = digamma(self.nu/2+1/2)-np.log(self.nu/2+1/2)
        return x+a1+a2
    
    def predict(self, x):
        """
        studentsのt分布
        """
        return (1/np.sqrt(np.pi*self.nu*self.sigma))*gamma((self.nu+1)/2)/gamma(self.nu/2)*np.power(1+1/(self.nu*self.sigma)*(x-self.mu)**2,-(self.nu+1)/2)

f =1000 #縮小率
data = sample/f
print("plot x sample")
# prepare model
students_t = StudentsT(nu=1)
gaussian = Gaussian()
gaussian_robust = RobustGaussian()
# maximum likelihood estimate
gaussian.fit(data)
print("Fit Gaussian")
gaussian_robust.fit(data)
print("Fit Robust Gaussian")
students_t.fit(data,100)
print("Fit Students T Distribution")

print(gaussian.mu, gaussian_robust.mu)

# plot results
x = np.linspace(0, 1920/f, 1920)
# 補正 : x方向にはf倍(圧縮の都合)。y座標は1/f(正規化の都合)*元データのスケールに合わせる
fujiyama_student = students_t.predict(x)/f*fujiyama.sum()+offset
fujiyama_gaussian = gaussian.predict(x)/f*fujiyama.sum()+offset
fujiyama_robust = gaussian_robust.predict(x)/f*fujiyama.sum()+offset

# plot
plt.figure(figsize=(9.6,6.4))
plt.plot(ns,color="blue",label="fujiyama") #富士山関数を正規化
plt.plot(x*f,fujiyama_student, label="students t",color="orange")
plt.plot(x*f,fujiyama_gaussian, label="normal",color="red")
plt.plot(x*f,fujiyama_robust, label="robust normal",color="green")

plt.xlim(0,1920)
plt.ylim(0,1280)
plt.legend()
plt.savefig("fitting.png",transparent=True)