ヒストスプライン平滑化をする際の節点の決め方ついて教えてください。

Question

与えられたヒストグラムをスプライン平滑化し、ヒストグラムの背後にある確率密度関数の推定ができるという論文を読みました。
具体的な内容を説明すると、
サンプル数n、区間[x_0, x_n]で等間隔の節点x_j (j=0, 1, ..., N)を決め、度数v_j∈[x_j-1, x_j] (j=1, 2, ..., N)のヒストグラムの面積s_j=v_j / n (j=1, 2, ..., N)を得る。このとき、ヒストグラムについての累積分布関数G_j (j=0, 1, ..., N)は、
・G_0 = 0
・G_j = Σ[j=1→i] s_j (j=1, 2, ..., N)
となり、このときS(x_j)=G_j (j=0,1, ..., N)を満たす3次スプライン関数S(x)が存在する。
このスプライン関数S(x)を微分するとヒストグラムの密度推定量φ(x)が得られ、この関数を未知の確率密度関数の近似関数としてみなすことができる。

とありました。
これを見る限り、ヒストグラムの節点の取り方が変わるとS(x)、φ(x)の値が変わると思います。
データの背後にある確率密度関数をより精度よく推定するための最適な節点の個数の決め方として何かアイデアはないでしょうか。
スタージェスの公式などは思いつきましたが、この公式を使って節点の個数を決めることが直接最適なφ(x)の推定につながるとは思えませんでした。

stomachman · Accepted Answer

No.8へのコメントについて。
　その論文は専ら「データから、スプライン曲線で表現された近似的な確率密度関数を構成する方法」だけを論じた論文であって、その確率密度関数をどう使うかについてはろくすっぽ考察されていない、ということなんでしょうかね？ もしそうなら、実践では大抵bootstrap法を使うんで、確率密度関数なんぞ要らない。ですから、その論文では少なくともbootstrap法を含む他の手段を批判して、スプラインを使う方がどこがどう優れているか、を論じていなくてはならない。それがキッチリ書いてないようなら、単に論文の著者が不勉強なだけであり、読むに値しません。ポイッ。
　あるいは、そうではなくて、その論文ではある具体的な課題を解決する手段として「確率密度関数をスプライン曲線で代用する」ということをやっている、というのであれば、じゃあどう使っているのか、をNo.8でお尋ねしたのですから、その答は得られなかったということで。

stomachman · Answer

No.7へのコメントについて。

ここまで拝見してきて、どうも、ご質問の課題に取り組むには、少なくとも1〜2箇月かけて何冊か教科書（近似理論、スプライン関数、最小二乗法など）を読み、基礎的概念を習得なさる必要があるだろう、という印象を持っています。
　ですが、何か急ぐ理由でもあるのでしょうね。だとすると、やはり一番重要なポイントをお尋ねすべきでしょう。すなわち「確率密度関数の近似曲線を描いたとして、さて、それを何にどう使うお積りですか？」ということです。というのも、用途によっては、そもそも確率密度関数の近似曲線を描く必要があるのかどうか、あるいは、そうすることが適切かどうか、が大いに疑わしいからです。（この質問にお答えがないようなら、今回で回答は終わりにします。）

> 「フィッティング」というのは

サンプルたちの全てをピッタリ通るような曲線を描けば、それは補間です。この場合、曲線を表すパラメータのうち、自由に決められるパラメータの個数(これを「自由度」という）は（特殊なタマタマの例を除いて）サンプルの個数と同じかそれ以上になります。たとえば4次スプライン関数は、節点の位置と、各節点間の多項式の0次,1次, 2次, 3次, 4次の係数、これらすべてがパラメータですけれども、節点において3次以下の微係数が連続であるべし、という制約が付いていますので、自由に決められるパラメータの個数は制約の数だけ少ない。これが自由度ということです。
　さて、自由度をサンプルの個数未満のある値に制限した曲線だと、サンプルたちの全てをピッタリ通ることはできない。しかし、サンプルたちの近くを通るようにならできる。これが近似です。なるべく良い近似になるようにパラメータを決めることをfittingと言います。実際には、近似の悪さを測るための定量的な尺度を決めておいて、この尺度を最小にするパラメータの値を決定する。
　Fittingに使われる曲線が理論的な根拠を持っているのなら、それはモデルと呼ばれます。「本来なら理論から導かれた通りの曲線（ただしパラメータの値は未知）になるはずなのに、サンプルが測定ノイズを含んでいるせいで、曲線がサンプルをピッタリ通らない」という考え方をしているわけです。ですから逆に言えば、fittingによってモデルのパラメータを決めることは「サンプルが持つ測定ノイズを取り除く」という意味を持ちます。たとえば、ある一つの物体Xの重量をN回測ったとすると、データ(x[k], y[k])（x[k]=k, k=1,2,...,N）が得られます。さて「Xの重量は変化しないはずだ」という理論に基づけば、 y[k] = c + ε[k]というモデルが決まる。cは（値が未知の）パラメータであり、だからこの曲線の自由度は1です。 ε[k]は残差で、このモデルの場合には「測定誤差の推定値」に当たります。このモデルを測定値{x[k], y[k]} にfittingしてcを決定したい。近似の悪さを測るための定量的な尺度として残差二乗和を使うことにすると、それは最小二乗法にほかならず、すなわち Σ(ε[k])^2 を最小化するようなcを探すわけです。（そうやって得られるcには「平均値」という名前が付いています。）
　 スプライン曲線は「弾性のある針金（その名前をsplineと言います）の上の何箇所かの位置を固定したとき、針金が自然に形作る形態（弾性エネルギーが最小になる形）」という理論的根拠を持っている。ですから、針金の形を計算するのに使うのなら、スプライン曲線はモデルです（もちろん、節点とは、位置を固定した場所のことです）。
　しかし（累積分布をスプラインで近似するという話のように）曲線がスプラインであるべき理論的根拠なんか全くない、という場合には、単に「似てればいいじゃん」「スプラインならいろんな曲線が作れるから便利じゃん」というだけのことです。（節点の個数や位置を定まらないのも、理論がないせいです。）なので、スプライン関数はグラフィクスや制御工学、機械工学などの技術者にとっては知っておくと良い知識ですが、それ以外の分野では重要ではありません。Fittingをやるのは、本来は、その現象に関する固有の理論とそれに導かれる固有のモデルがあるからこそであり、ならば当然、スプラインなんぞに出番はないわけです。

> そもそもスプライン近似と最小二乗法は同じことなのか

違います。まるで「外国旅行に行くことと飛行機に乗ることは同じか」というぐらいにトンチンカンなご質問です。
● fittingにおいて、サンプルたちをどんな関数で近似するか。もしスプライン関数で近似するのなら、それはスプライン近似です。
● fittingにおいて、近似の悪さの度合いをどう測るか。もし、残差2乗和を尺度にするのなら、それは最小二乗法です。

> （残差）というのは上の座標のy座標と何の残差を指す

ヒストグラムを作らずに、最初っから累積分布を使う方の話だとして回答します。ご質問の課題の場合、小さい方からk番目のサンプルのy座標はk-1/2であって、誤差なんかありません。誤差を持つのはx座標の方ですね。だから、この場合のfittingは、サンプル(x[k], y[k])(k=1,2,...,N)について
　　x[k] = (y[k]の区分的多項式）+ ε[k]  （ε[k]は残差）
において、残差が何らかの意味で（例えば最小二乗法ならΣ(ε[k])^2が）最小になるように区分的多項式を決定する、というのが適切です。もちろん、この区分的多項式には(0,0)を通るという条件がつく。それだけでなく、0以上1以下であって単調増加だという条件も付きますから、最小二乗法だけでうまく行くとは限りません。確率密度関数を計算するには、まずはこの区分的多項式の逆関数（単調なので、必ず一価関数になります）を作って、その導関数を計算することになります。もちろん、その結果は（区分的ではあるが）多項式にはなりません。
　それじゃダメなんだ確率密度関数を区分的多項式で表したいんだ、というのがもし絶対の要求事項ですと：
(A) 上記の式の (y[k]の区分的多項式）の部分を (x[k]の区分的多項式の逆関数）に置き換える。（すると、非線形最小二乗法の問題になって、反復解法で解く必要があり、計算量がかなり大きくなります。）
(B)いやそんな本気にならなくたって、どうせ「単に似てればいいじゃん」ってだけだし、ということなら
  y[k] = (x[k]の区分的多項式）+ ε[k]  （ε[k]は残差）
を使う。

stomachman · Answer

No.6へのコメントひとつ目については、ナンセンスにしか見えんです。

> 区間k=5,x1~x5とする

ということはヒストグラムを描いて（情報を捨てて）からfittingをやる話の方でしょうかね。

> 節点2個(x1<m1<x2, x3<m2<x4)

わずか5個のサンプル点に対して節点が2つもあるなんてのは、パラメータの個数の方がサンプル点の個数より多いわけですから明らかに過剰ですが、ま、それはさておき

> 最適な節点数を求める際に計算する(残差)^2/分散は、

本当にワケわからん。節点数を2個だとすでに決めちゃった上での話ですよね？「最適な節点数」もヘッタクレもないでしょうに。

> [{s1(x1)-F(x1)}^2 / (x1の分散)] +

　なんだこれ。まず分母の (x1の分散)を推定する方法はどうするんですか。もしかして、「区間内に入るサンプル数の分散は、区間内に入るサンプル数の期待値に等しい」ということから、「実際に区間内に入ったサンプル数」で代用する、ってことでしょうか？だったら単に「その区間の頻度」と言えば良いのにな。
　分子は一層の謎で、

> si(x)はスプライン関数, Fはデータ点を重み付き最小二乗法によって解いた関数

全然意味わからんです。「重み付き最小二乗法によって解いた関数」なるものこそがsi(x)じゃないんですかい？？？？「データ点」って一体何ですか？

ーーーーーーーーーーーーーーーーーーーーーーーー

No.6へのコメント二つ目について。

こっちは、ヒストグラムを使わない話の方ですかね。（一体どっちをやりたいんだか、どうにも要領を得ないなあ。）

> {k,g(k)}

というのがそもそも意味わからんです。「g(k)はkまでの累積度」とかおっしゃる「kまで」って何？

　例えば、データを小さい順に並べたら
　　0.4, 0.8, 1.5, 2.1, 2.2, ...
になったとします。この場合に{k,g(k)}なるものは具体的にどうなるんでしょうかね。

　ちなみに、No.6のグラフの一番上のやつをこのデータに基づいて描くと、タテ軸方向の段差の高さはどれも1です。そして段差は、ヨコ軸が0.4, 0.8, 1.5, 2.1, 2.2, ...のところに存在します。小さい順でm番目までのデータの累積度数がmであることは言うまでもありません。グラフの一番下のやつにプロットしてある左からj番目の点の座標を(x[j], y[j])とすると、もちろん x[j]は「小さい順でj番目のデータの値そのもの」であり、y[j] = j-1/2 ですから、j=1,2,3,4,5,…について
　　(0.4, 0.5), (0.8, 1.5), (1.5, 2.5), (2.1, 3.5), (2.2, 4.5), …
となります。

　それで、kとかg(k)とかは一体何を表しているのかをご説明くださいな。

stomachman · Answer

No.5へのコメントについて。

> 累積度数)のスプライン関数S(x)の節点を最適化

ご質問は、ヒストグラムを作ってから（従って、各区間の中での分布がどう偏っているかの情報は失ってしまった後で）曲線をfittingする、という話だったと思う。もしそうなら、累積度数分布に曲線をfittingしたってしょうがなくって、最初っからヒストグラムにfittingして確率密度関数を推定すればいいんです。

一方、そういうコトじゃなくて、元のデータから確率密度関数を推定するのが目的であって、ヒストグラムなんかどうでもよろしい、というご質問なのだとしたら、ハナシが違ってきて、累積密度関数を考えることに意味があります。
　この場合、（人工的に決めた）等間隔の区間、なんてものはない。データひとつひとつが階段を形作る折れ線グラフです。階段の横幅は不揃いになる。最終的にn次スプラインで確率密度を表現したいのなら、この（登りしかない）階段をn+1次スプラインで近似しておいて、その導関数を確率密度関数の推定値とすることになります。このとき、

> 分布関数g(i)は区間ごとに不連続な定数関数

いいえ。階段の横幅が広いところは「データがないために平たい」んですから、そんなもんに曲線をfitさせる理由はありません。「定数関数」になる部分は完全に無視して良いということです。なので、こんな風にやる↓

stomachman · Answer

> ① プライン近似に使用するヒストグラムの各区間の大きさが異なる

いいえ。ヒストグラムの各区間の幅と、節点の位置とは無関係です。

> ②まとめると全体の流れとしては、
> 1.AICのような尺度を使って節点の個数を求める
> 2.節点の最適な配置を求める
> 3.フィッティング

違います。

節点の個数のいろいろなバリエーションについて、以下を繰り返す。
　　節点の配置をテキトーに決める。
　　フィッティングし、残差の程度を評価する。
　　以下を繰り返す。
　　　　評価を改良するように、節点の配置を修正する。
　　　　フィッティングし、(AICのような尺度を使って)評価する。
　　　　評価が収束したら繰り返しを抜ける。
最も良い結果になった節点の個数と配置、多項式のパラメータを報告する。
終わり。

kamiyasiro · Answer

等間隔の節点って、「連続修正点」のことではないでしょうか。

ヒストグラムは、詳細な情報を欠落させてしまっています。

しかし、ビンとビンの中点（境界線）だけは、正しい累積確率を持ちます。

ですから、スプライン近似曲線を引くのであれば、その正確な点を通すように引くべきです。

図で説明します。

ヒストグラムの累積分布曲線を描くと、ヒストグラムの階級値のところにデータが溜まっているので、垂直に累積確率が上昇しますが、次の階級値までは一定の累積確率（水平）を維持します。つまり階段状のグラフになります。

その階段のステップの中点が「連続修正点」（呼び名は色々あります）になります。

図では赤い点です。

あと、どれくらいの刻みにするかは、ノイズ（観測ばらつき）との兼ね合いでしょう。

まあ、細かく取っても、累積確率なので、少なくとも単調増加にはなります。

stomachman · Answer

No.1もう一個ミスプリ修正。
「N^3に反比例」じゃなくて、「N^(1/3)に反比例」。

stomachman · Answer

No.1ミスプリ修正。
　{k, n(k)} とすべきところが　{k,h(k)} になっちゃってました。

stomachman · Answer

> スタージェスの公式
　お書きの通り、この「公式」はヒストグラムとはまるで関係ありません。ヒストグラムのヨコ軸の刻み幅を決めるのに使う、とか言う人もいるようですが、全くの誤謬です。（そんなことしたら、N数が大きいほど幅が過大になり、折角のデータを無駄にします。）ヒストグラムの幅をどうするかは、一見ツマラナイことのようでも実務上はかなり重要で、ご質問の場合にも、幅をどうするかによって、得られる結果が全然違ってくるだろうことに留意すべきでしょう。少なくとも、その幅は(N^3)に反比例するんでなくちゃいけません。と、それはさておき。

そもそも、何のためにスプライン近似をやるんでしょうか。滑らかな確率密度関数が欲しい理由として考えられることと言えば、与えた確率密度関数に従う乱数（ベイズ統計学で必須）を生成したい、ということかしらん、と想像します。確率密度関数が1変数関数 φ(x)のとき、その分布関数 Φ(x) = ∫{-∞～x} φ(t) dt の逆関数Φ*を用意して、一様乱数XをΦ*(X)で変換するだけで、φ(x)に従う乱数が生成できますね。この場合には滑らかなφなんてものは不要です。でも残念ながら、φが2変数以上の場合にはこの手は使えない。ところがMCMC法等を使えばφに従う乱数が効率よく生成できる。ただし、その際にφの偏微分を計算する必要があり、だからφは滑らかでないとまずい…という文脈でなら、φが滑らかであることに意味があるでしょう。
　しかしこれは、φをサンプルデータのヒストグラムから作る理由にはなりません。なぜなら単に「サンプルデータセットから等確率でランダムに1個のサンプルを選んでくる」ということをやるだけで、「サンプルデータのヒストグラムに従う乱数」が生成できる（その上、ヒストグラムを作る必要すらないので、「横軸の刻みの幅」の心配も無用だ）からです。φなんか作るまでもありません。（ノンパラメトリックな統計では、このようにしてサンプルデータセットそのものを母集団の分布の代用として利用するのはよくやることで、総じて「ブートストラップ法」と呼ばれます。）
　そうしますと、「ヒストグラムを作った上でスプライン関数で平滑化する」ことのメリットと言えば、ヒストグラムを記述するのに必要なbit数が節約できることぐらいしか思いつかない。コンピュータの能力がごく限られていた時代の発想のように思われます。

ところで、スプライン近似することには理論的根拠がないので、統計的推論や論証に使うと言い訳が残る（例えば「論証に不都合な凸凹を隠すための恣意的な加工じゃないのか？」と批判される）。しかし、予測や制御に使う分にはご随意に、という位置付けになるでしょうかね。スプライン関数の節点を漫然と等間隔に並べるとか、サンプル点のどれかに一致させるとか、そういうヘボなやり方をするのも、「予測や制御に使う分にはご随意に」です。と、用途の話はさておき。

ともあれホンキでプライン近似するんであれば、ご質問でご指摘の通り、近似関数とデータとのズレを定量評価して、多項式の次数や節点の個数やその配置を最適化すべきでしょう。
　お読みの文献はそこまで論じてはいないようで、そんなのだったら「ちょっとやってみた」だけのこと。誰でも考えるし、誰でもすぐやれる。そういうものを「確率密度関数の推定」だなんて称するのはオコガマシイですねー。（なお、分布関数の微分が確率密度関数なのは、分布関数の定義から自明です。）だから、その文献のその部分に、さしたる意味があるようには思われません。その上、近似することにこれと言って明確なメリットがあるわけでもなさそうで、一体どんな目的でそんな近似をやろうというのでしょうか。となると、その文献全体の意義もどうなんでしょうかね。というのも、これまたさておき、ようやく本題。

ヒストグラム全体の曲線下面積はサンプル数Nに一致しています。で、ある区画kの面積がn(k)だとします。N個中のn(k)個が区画kに入った、ってことですが、仮に同じ実験をたくさん繰り返したとすれば、同じ区画に入るサンプル数は毎回バラツクでしょう。そこで、「確率p(k)で区画kに入る」という独立な事象をN回繰り返したらn(k)はどんな分布に従うか、と考える。するとよく知られている通り、n(k)は二項分布 B(N,p(k); n(k)) = (NCn(k)) (p(k)^n(k))(1-p(k))^(N-n(k)) に従う。その平均はNp(k), 分散はNp(k)(1-p(k))であり、Nが大きいときには正規分布でよく近似できる。そこでp(k)をn(k)/Nと推定すれば、区画kに入る確率p(k)[いわゆる「真値」]は平均n(k)/N、分散n(k)(N-n(k))/Nの正規分布に従う、と近似できます。かくして「分散n(k)(N-n(k))/Nの正規分布に従う測定誤差を持つデータ {k, h(k)}が与えられた」と捉え直せます。そして、kごとに誤差の幅が異なることを考慮して、スプライン（であろうが、他のナニカであろうが）によってこのデータ{k, h(k)}を近似することを考える…というのがご質問の趣旨かと思います。

じゃあどうするか。方針だけ考えてみます。スプライン関数の節点と次数が与えられたとき、「分散がkに依って異なる正規分布に従う誤差を持つデータ {k, h(k)}」に最小二乗法で多項式をフィッティングするのは線形の問題なので、簡単な線形代数の計算で解けます。次に、節点の個数が決まっているとして、その最適な配置は非線形の問題であり、探索で反復的に改善していく必要があります。隣接する節点間の区間での(残差^2)/分散の程度が、どの区間についても均一になるように、ちょっとずつ調節していくわけです。
　さらに、節点の数と多項式の次数、という離散的なパラメータをどう選ぶのが最適か。これらを変えるとスプライン全体の自由度が変化しますから、近似の良さを単純には比較できない。そこで情報量の観点から、例えばAICのような尺度を使って評価する必要があるでしょう。

で、えーと、そこまで手を掛けるだけのメリットがあったっけ、と思い返してみると、さ、どうなんでしょうかね。

ヒストスプライン平滑化をする際の節点の決め方ついて教えてください。

No.8へのコメントについて。

No.7へのコメントについて。

No.6へのコメントひとつ目については、ナンセンスにしか見えんです。

No.5へのコメントについて。

> ① プライン近似に使用するヒストグラムの各区間の大きさが異なる

等間隔の節点って、「連続修正点」のことではないでしょうか。

No.1もう一個ミスプリ修正。

No.1ミスプリ修正。

> スタージェスの公式

関連するカテゴリからQ&Aを探す

デイリーランキングこのカテゴリの人気デイリーQ&Aランキング

マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング