組合せバンディットを用いたコグニティブ無線における
グループ形成方策
Coalition Formation Policy in Cognitive Radio Networks
by Combinatorial Bandits
飯塚 翔
1∗川原 純
1笠原 正治
1Sho Iizuka
1Jun Kawahara
1Shoji Kasahara
11
奈良先端科学技術大学院大学情報科学研究科
1
Graduate School of Information Science, Nara Institute of Science and Technology
Abstract: In cognitive radio networks (CRN), secondary users (SUs) find idle times of a wire-less channel and utilize them secondarily in order to eliminate the shortage of usable frequency bands. Cooperative spectrum sensing (CSS), with which multiple secondary users form a coalition and sense a wireless channel cooperatively, is a scheme to improve sensing performance by SUs. In order to reduce sensing errors by CSS, it is necessary to form a coalition in consideration of the sensing performance of each SU. In this study, we discuss coalition formations in the absence of prior information about the performance of each SU, that is, each SU knows neither his/her sensing performance nor other SUs’ one. In order to minimize sensing errors under this constraint, we propose a statistical coalition formation policy based on the multi-armed bandit problem. In the proposed method, coalition selections are regarded as arm selections, and Thompson sampling performs coalition selections by considering exploitation-exploration trade-offs. Numerical exam-ples show that the proposed method particularly reduces sensing errors in a scenario where there are many SUs and the performance difference among SUs is large.
1
はじめに
無線資源の有効利用を目的としてコグニティブ無線シ ステム(CRN) [1, 2]が提案されている.CRNではチャ ネルを割り当てられた一次ユーザ(PU)によるチャネル の利用状況を二次ユーザ(SU) がセンシングし,PUに よってチャネルが利用されていない期間にSUがチャネ ルを利用することによって,チャネルの利用効率を改善 することを目的としている. SUがチャネルをセンシングする際には誤検知(missdetection)と誤警報(false alarm) という二種類の誤り
が生じる.誤検知は,実際にはチャネルが利用中(busy) であるにも関わらず,SUはチャネルが利用中でない (idle)と判断する誤りである.誤警報は,実際にはチャ ネルが利用中でない(idle)でないにも関わらず,SUは チャネルが利用中である(busy)と判断する誤りである. ∗奈良先端科学技術大学院大学情報科学研究科 〒630-0192 奈良県生駒市高山町 8916-5 E-mail: [email protected] 誤検知が生じると,PUの通信との衝突が生じ,PUの 通信品質が悪化する.誤警報が生じると,SUが通信可 能であったはずの機会を失い,SUのスループットが低 下する. SUによるチャネルのセンシング精度を改善するため に,複数のSUがグループを形成し,個々のSUのセン シング結果を共有することでセンシング精度を改善する ことを目的とした協調センシング (CSS) [3, 4]が提案 されている.協調センシングのためのグループ形成手法 として,ゲーム理論に基づくアプローチが提案されてい る[5, 6].このアプローチでは,グループ内の個々のSU のセンシング精度から計算されるグループのセンシング 精度を用いて効用が定義され,個々のSUは独立に自分 自身の効用を最大化するようにグループを形成する.し かしながら,このアプローチは個々のSUは自分自身の みならず他のSUのセンシング精度も把握しているとい うことを仮定している. 本稿ではSUのセンシング精度についての事前情報を 必要としない統計的グループ形成法を提案する.この手 人工知能学会研究会資料 SIG-FPAI-B509-15
法ではグループ形成の問題を多腕バンディット問題とし て定式化し,誤検知・誤警報を削減することを目的とし てグループ形成の試行錯誤が行われる.試行錯誤の各ス テップでは,SU自身が得ることができる情報のみから 個々のSUのセンシング精度を推定した後,アルゴリズ ムによって提案されたグループで協調センシングを行 い,通信を試みる.その後,協調センシングの結果と通 信の結果を用いて個々のSUのセンシング精度について の信念を更新する.提案手法はコグニティブ無線システ ムにおいてSUの設置場所の電波環境が不明な場合や, SUが劣悪な環境に置かれることでセンシング精度が劣 化するような場合に有用であると考えられる. 本稿の構成は以下のとおりである.第2節ではモデル の説明と問題の定式化を行う.第3節ではベースライン であるε-first戦略に基づくグループ形成方策とその問 題点,および提案手法である組合せバンディットに基づ く方策について述べる.第4節では数値実験により提案 手法の評価を行う.第5節ではまとめと今後の課題を述 べる.
2
モデルと定式化
本稿では,SUがグループを1つ形成し,そのグルー プでチャネルの協調センシングを行い,チャネルを用 いた通信を試みるという手続きを時間ステップの集合 T = {1, 2, . . . , T }それぞれで行うことを考える.N を システム全体のSU数,N = {1, 2, . . . , N}をSUの集合 として,形成されるグループはN中のM (1≤ M ≤ N) 台のSUからなる.N の中で1はグループ形成を主導 するSUであり,必ずグループに含まれると仮定する. すなわち,t回目の時間ステップに形成されるグループ Ctは Ct={1} ∪ C, where C ⊆ N \ {1} and |C| = M − 1 と表される. t回目の時間ステップにおけるグループCtの形成後, SUは協調センシングを行い,チャネルを用いた通信を 試みる.協調センシングの結果,チャネルがidleであ ると判断されたならば,グループはチャネルを利用した 通信を試みる.一方で,チャネルがbusyであると判断 されたならば,グループはその時間ステップ中は通信を 試みない.グループが通信を試みたとき,PUがチャネ ルを利用していなかったならば通信は成功する一方で, PUがチャネルを利用していたならばPUの通信と衝突 が生じることにより通信は失敗する.また,PUがチャ ネルを利用していなかったにも関わらずグループがチャ ネルをbusyと判断した場合,グループは本来は通信可 能であったはずの通信機会を失う.PUによるチャネル 利用はパラメータρ (0≤ ρ ≤ 1)のベルヌーイ過程に従 う.すなわち,各時間ステップにおいてチャネルは確率 ρでbusyであり,確率1− ρでidleである.ここでρ をPUのチャネル利用率と呼ぶ.SUはρについて事前 情報を持たないと仮定する. i (i∈ N )番目のSUの検知率をλ(i)D (0≤ λ(i)D ≤ 1), 誤警報率をλ(i)F (0≤ λ(i)F ≤ 1) と表す.ここでλ(i)D はPU がチャネルを使用していた場合にi番目の SUが チャネルをbusyと判断する確率であり,λ(i)F はPUが チャネルを使用していなかった場合にi番目のSUが チャネルをbusyと判断する確率である.SUはすべて のi∈ N についてλ(i)D とλ(i)F について事前情報を持た ないと仮定する. 協調センシングによる意思決定はグループ内の個々の SUのセンシング結果を用いてk-out-of-Nルール [7]に よって行われる.これは,グループ内のk台以上のSU がチャネルをbusyと判断した場合に,グループがチャ ネルをbusyと判断するというルールである.誤検知と 誤警報がSU間で独立に生じると仮定すると,グループ の検知率ΛD(Ct)と誤警報率ΛF(Ct)はそれぞれ ΛD(Ct) = ∑ C⊆Ct,|C|≥k ∏ i∈C λ(i)D ∏ i∈Ct\C ( 1− λ(i)D ) ΛF(Ct) = ∑ C⊆Ct,|C|≥k ∏ i∈C λ(i)F ∏ i∈Ct\C ( 1− λ(i)F ) と表される. SUは全時間ステップT の中でPUの通信と衝突す る回数を減らしつつ,通信に成功する回数を増やすこと を目的として振る舞う.よって,時間ステップの集合T におけるグループ形成の試行錯誤の良さは,効用関数 U (T ) = αNS(T ) − (1 − α)NMD(T ) によって評価される.ここでNS(T )は通信に成功した 回数,NMD(T )はPUの通信と衝突が発生した回数で ある.αはNS(T )とNMD(T )のバランスをとるパラ メータであり,SUにとって既知の情報である.この効 用U (T )の期待値は,1ステップあたりの期待値 u(Ct) = α(1− ρ)(1 − ΛF(Ct))− (1 − α)ρ(1 − ΛD(Ct)) を用いることで E[U(T )] =∑ t∈T u(Ct) と表すことができる.すなわち,ρ, λ(i)D, λ(i)F が既知の 場合は1ステップあたりの期待値u(Ct)を最大にするグ ループを選び続けることでU (T )の期待値を最大にする ことができる.
3
提案手法
個々のSUのセンシング精度についての事前情報が 得られず,1 ステップあたりの効用の期待値u(Ct)を 最大にするグループが不明な場合の素朴な方策として, ε-first戦略[8]に基づく方策が考えられる.この方策で は全体の時間ステップT を探索ステップ{1, 2, . . . , T′} と活用ステップ{T′+ 1, T′+ 2, . . . , T}に分割する.探 索ステップではM 台のSUを一様にランダムに選択し, そのグループで協調センシングを行い,通信を試みる ことで,SUのセンシング精度を推定するために必要な データを得る.探索ステップが終了した時点で,それま でに得られたデータからρ, λ(i)D, λ(i)F を推定する.活用 ステップでは推定されたρ, λ(i)D , λ(i)F を用いてu(Ct)を 最大にすると推定されたグループを選択する. この方策は,SU自身が得ることができる情報のみか らSUのセンシング精度を推定し,そこから高いセンシ ング精度をもつと推定されたグループを選択することが できる.しかしながら,探索ステップの長さを決めるT′ の設定において探索と活用のトレードオフが存在する. すなわち,T′を大きく設定した場合,推定に利用可能な データの増大によってρ, λ(i)D, λ(i)F の推定は正確になる が,活用ステップが短くなることによってu(Ct)を最大 にすると推定されたグループを選択することができる回 数が減少する.一方でT′を小さく設定した場合,u(Ct) を最大にすると推定されたグループを多く選択すること ができるが,学習データの不足によりρ, λ(i)D, λ(i)F の推 定は不正確になる. このような探索と活用のトレードオフがある状況にお いて最適な方策を考える問題として多腕バンディット 問題[9]がある.多腕バンディット問題とは,報酬の確 率分布が異なる多数のスロットマシンのアームがあり, アームを動かすことができる回数に制限があるときに, 得られる報酬の合計を最大化するようなアームの選択 を行う方策を考える問題である.この問題においても, 個々のアームから得られる報酬の確率分布を正確に推定 しようとすると,報酬の期待値が最大と推定されたアー ムを選択できる回数が減少してしまうというトレードオ フがある. 提案手法では,多腕バンディット問題をグループ形成 の問題に応用するために,形成するグループの選択を アームの選択と考え,多腕バンディット問題に対する方 策を適用する.しかしながら,形成可能なグループの組 合せのそれぞれを単純に1つのアームとみなす方法には 2つの問題点がある.1つめは,全SU数N やグルー プ内に含まれるSU数M の増加にともなって形成可能 なグループの組合せの数が急激に増加するため,アーム の本数が増大してしまうことである.2つめは,あるグ ループを形成し,そのグループに含まれるSUのセンシ ング精度についてのデータが得られたとき,そのグルー プ以外の他のグループについての情報が部分的に得ら れているにも関わらず,それを活用していないことであ る.たとえば,SU{1, 2, 3}からなるグループで協調セ ンシングを行ってSUのセンシング精度に関する情報が 得られた場合に,いくつかのSUを共有する{1, 2, 4}や {1, 3, 5}といったグループについての情報も部分的に得 ることができる.1つのグループを1つのアームとみな す単純な方法では,あるアームを選んだ際に他のアーム について得られる情報を活用しないため,本稿の問題設 定において効果的に機能しないと考えられる. 提案手法では,組合せバンディットに対する Thomp-son sampling [10]に基づいてグループの選択を行う. Thompson sampling [11, 12]では,i番目のアームを引 いたときの過去の報酬の履歴xiからアームの平均報酬 µi の推定値の事後分布P (µi | xi)を推定する.この事 後分布から,それぞれのアームiについてそのアームが 最適である確率P (µi > maxj̸=iµj | xi, xj)を計算し, その確率に一致するようにアームを選択する.ここで P (µi> maxj̸=iµj | xi, xj)を解析的に求めることは必 ずしも簡単ではないが,それぞれのアームiについて事 後分布P (µi | xi)からµiのサンプルを生成し,µi の サンプルが最も大きいアームiを選択することで等価な 方策が実現できる.Thompson samplingは組合せバン ディットに対しても優れた性能を示すことが実験的に示 されており [13],単純な場合では性能の最適性が理論的 に示されている[14]. 提案手法の擬似コードをAlgorithm 1に示す.この 手続きでは,1ステップの中で最初にEstimateVB(·) でρ, λ(i)D , λ(i)F の事後分布の推定を行い,事後分布から の乱数サンプルを生成したあと,生成された乱数サンプ ルを用いてFindOptimal(·)でu(Ct)を最大にするグ ループCtを探索する.その後,Collect(·)で協調セ ンシングのためにグループCtに含まれる個々のSUの センシング結果を収集し,Communicate(·)で協調セ ンシングの結果に応じてチャネルを用いて通信を試み る.グループCtに含まれる個々のSUのセンシング結 果s(i)t は,t回目の時間ステップにi番目のSUがチャ ネルをbusyと判断した場合には1,チャネルをidleと 判断した場合には0となる.Communicate(·)の結果 ytは,グループがチャネルをidleと判断し,通信を試み た結果,通信に成功した場合には1,通信を試みたがPU の通信と衝突が発生した場合には2,グループがチャネAlgorithm 1 組合せバンディットに基づくグループ 形成
procedure CoalitionFormation(N , M , T , α)
Create lists CS [1, T ], S[1, T ], and Y [1, T ] to storeCt, st, and yt
for t = 1, 2, . . . , T do
(a[·], b[·]) ← EstimateVB(CS, S, Y ) Sample ρ∼ β(ρ; a [ρ] , b [ρ])
for i = 1, 2, . . . , N do
Sample λ(i)D ∼ β(λ(i)D ; a[λ(i)D], b[λ(i)D]) Sample λ(i)F ∼ β(λ(i)F ; a[λ(i)F ], b[λ(i)F ])
end for Ct← FindOptimal(N, M, α, ρ, λD, λF) ▷ 乱数サンプルρ, λ(i)D, λ(i)F を用いてu(Ct)を最大化す るグループCtを探索する for all i∈ Ctdo s(i)t ← Collect(i) ▷ i番目のSUのセン シング結果を取得する end for yt← Communicate(st) ▷ センシング結果 stから通信の意思決定を行い通信を試みる CS [t]← Ct S[t]← st Y [t]← yt end for end procedure ルをbusyと判断して通信を行わなかった場合には3と なる.
EstimateVB(·)でのρ, λ(i)D, λ(i)F の事後分布の推定
では,PUのチャネル利用状況の一部を隠れ変数と扱っ て推定する必要がある.yt = 1の場合にはチャネルが idle,yt = 2の場合にはチャネルがbusyであることが 確定する一方で,yt= 3の場合にはチャネルがbusyで あることをSUが正しくセンシングしたのか,チャネル がidleであるにも関わらず誤警報によってSUが誤って busyと判断したのか不明であるためである.このため, 事後分布の推定にはマルコフ連鎖モンテカルロ法を用い た推定[15]や隠れ変数を考慮した変分ベイズ推定[16] が必要となる.本稿では後者の変分ベイズ推定を用いた 方法を用いる.この方法ではρ, λ(i)D, λ(i)F の事後分布を ベルヌーイ分布の共役事前分布 [17]であるベータ分布
β(·, ·)として,それぞれβ(a [ρ] , b [ρ]), β(a[λ(i)D], b[λ(i)D]),
β(a[λ(i)F ], b[λ(i)F ])とする.ここでa[·], b[·]は角括弧中の 変数が従うベータ分布のパラメータを表す.推定では最 初にa[·], b[·]を1で初期化したあと,a[·], b[·]の変化が 十分小さくなるまで変分Eステップと変分Mステップ を繰り返す.変分Eステップではyt= 3となるtそれ ぞれについて rt= 1/(1 + exp(−ηt)) を計算する.ここで ηt= (ψ(a[ρ])− ψ(b[ρ])) +∑ i∈Ct s(i)t ( ψ ( a[λ(i)D] )
− ψ(a[λ(i)D] + b[λ(i)D] )) +∑ i∈Ct (1− s(i)t ) ( ψ ( b[λ(i)D ] )
− ψ(a[λ(i)D ] + b[λ(i)D ] )) −∑ i∈Ct s(i)t ( ψ ( a[λ(i)F ] )
− ψ(a[λ(i)F ] + b[λ(i)F ] )) −∑ i∈Ct (1− s(i)t ) ( ψ ( b[λ(i)F ] )
− ψ(a[λ(i)F ] + b[λ(i)F ] )) である.ψはディガンマ関数 [18]を表す.変分Mス テップではρ, λ(i)D, λ(i)F の事後分布のパラメータをそれ ぞれ a[ρ] =|T2| + ∑ t∈T3 rt+ 1 b[ρ] =|T1| + ∑ t∈T3 (1− rt) + 1 a[λ(i)D] = ∑ t∈T2(i) s(i)t + ∑ t∈T3(i) rts (i) t + 1 b[λ(i)D] = ∑ t∈T2(i) (1− s(i)t ) + ∑ t∈T3(i) rt(1− s (i) t ) + 1 a[λ(i)F ] = ∑ t∈T1(i) s(i)t + ∑ t∈T3(i) (1− rt)s (i) t + 1 b[λ(i)F ] = ∑ t∈T1(i) (1− s(i)t ) + ∑ t∈T3(i) (1− rt)(1− s (i) t ) + 1 と計算する.ここでTj (j = 1, 2, 3)は{t | t ∈ T ∧ yt= j}を表し,T(i) j (i∈ N , j = 1, 2, 3)は{t | t ∈ Tj∧ i ∈ Ct}を表す. FindOptimal(·)でのu(Ct)を最大にするグループ Ctの探索では,N が小さい場合には全探索によって探 索が行えるが,N が大きい場合には全探索での探索は 時間計算量の観点で難しい.そこで,本稿ではN が大 きい場合にAlgorithm 2に示すヒューリスティックを 利用する.これは最初に Copt = N の状態から始め,
u(Copt\ {i})を最大にするSU iをCoptから取り除くこ
Algorithm 2 u(Copt)を最大にするグループCoptを探
索するヒューリスティック
function FindOptimal(N , M , α, ρ, λD, λF)
Copt← N
while |Copt| > M do
CS ← {Copt\ {i} | i ∈ Copt\ {1}}
Copt← arg maxC∈CSu(C)
end while returnCopt end function 表1 評価実験におけるパラメータ設定 パラメータ名 設定値 N 50 M 5 k 3 ρ 0.5 λ(i)D (0.7, 0.95, 0.95, 0.95, 0.95, 0.7× 45) λ(i)F (0.3, 0.05, 0.05, 0.05, 0.05, 0.3× 45) α 0.2
4
評価実験
本節では,提案手法の評価を行う.ベースラインは第 3節で述べたε-first戦略に基づく方策である.SUは時 間ステップ集合T のそれぞれの各ステップにおいて方 策に従ってグループを形成し,通信の成功回数NSと PUの通信との衝突回数NMDからなる効用関数U (T ) を最大化することを目的とする. 評価実験におけるパラメータ設定を表1に示す.この 設定はSUの台数が多く,SU間のセンシング精度の差が 大きい状況を想定している.SU{2, 3, 4, 5}のセンシン グ精度は検知率・誤警報率ともにSU{1, 6, 7, . . . , 50}の ものよりも優れているため,SUのセンシング精度が既知 であればSU{1, 2, 3, 4, 5}からなるグループがE[U(T )] を最大にする.このグループの検知率は0.995,誤警報 率は= 0.00454である.一方でSU{1, 6, 7, 8, 9}からな るグループはE[U(T )]を最小にするが,このグループ の検知率は0.837,誤警報率は0.163である.ベースラ インであるε-first戦略に基づく方策ではT′の値を200, 300, 400に設定した3種類で実験を行う.実験は乱数 生成器のシード値を変更して70回行う. 図1に効用関数の変化を示す.この図から,すべての 時間ステップにおいて提案手法はベースラインを上回っ ていることがわかる.ε-first戦略に基づく方策では時間 0 200 400 600 800 1000 # of steps 0 20 40 60 80 cummulative reward -first (200) -first (300) -first (400) bandit 図1 効用の時間変化 0 200 400 600 800 1000 # of steps 0 5 10 15 20 25 30 35# of occurences of miss detection
-first (200) -first (300) -first (400) bandit 図2 誤検知の累積発生回数 ステップT′ において傾向が変化しているが,これは時 間ステップT′ において探索ステップから活用ステップ に移行し,探索ステップで得られたデータを用いて最 適と推定されたグループを形成するようになるからで ある. 図 2, 3にそれぞれ誤検知と誤警報の累積発生回数を 示す.この図から,本実験設定においては,1000回の時 間ステップが経過した時点で,提案手法はベースライン と比較して誤検知と誤警報の発生回数を約3分の1に削 減していることがわかる.
5
おわりに
本稿では,個々のSUのセンシング精度についての事 前情報が得られない状態での協調センシングのためのグ ループ形成法について議論し,組合せバンディットに基 づくグループ形成方策を提案した.提案手法では組合せ バンディットに対するThompson samplingに基づいて グループの選択を行うため,変分ベイズ推定によりSU のセンシング精度の事後分布の推定を行い,事後分布か0 200 400 600 800 1000 # of steps 0 5 10 15 20 25 30 35
# of occurences of false alarm
-first (200) -first (300) -first (400) bandit 図3 誤警報の累積発生回数 ら生成した乱数サンプルを用いて,効用関数の1ステッ プあたりの期待値を最大化するグループをヒューリス ティックによって探索した.今後の課題としては,本稿 の問題設定においては1グループのみの形成を目的とし ていたが,これを複数グループの形成に拡張することが 考えられる.
謝辞
本研究の一部は,科研費基盤 (B) 15H04008および SCAT研究助成による支援を受けて実施している.参考文献
[1] J. Mitola and G.Q. Maguire, “Cognitive radio: making software radios more personal,” IEEE Pers. Commun., vol.6, no.4, pp.13–18, 1999. [2] S. Haykin, “Cognitive radio: brain-empowered
wireless communications,” IEEE J. Sel. Areas Commun., vol.23, no.2, pp.201–220, 2005. [3] A. Ghasemi and E.S. Sousa, “Collaborative
spec-trum sensing for opportunistic access in fad-ing environments,” IEEE DySPAN, pp.131–136, 2005.
[4] I.F. Akyildiz, B.F. Lo, and R. Balakrishnan, “Cooperative spectrum sensing in cognitive ra-dio networks: A survey,” Phys. Commun., vol.4, no.1, pp.40–62, 2011.
[5] W. Saad, Z. Han, T. Basar, M. Debbah, and A. Hjorungnes, “Coalition formation games for col-laborative spectrum sensing,” IEEE Trans. Veh. Technol., vol.60, no.1, pp.276–297, 2011.
[6] T. Nishida, M. Sasabe, and S. Kasahara, “Max-imizing communication opportunity for
collabo-rative spectrum sensing in cognitive radio net-works,” ITNAC, pp.1–6, 2017.
[7] C. Sun, W. Zhang, and K. Ben Letaief, “Cooper-ative spectrum sensing for cognitive radios under bandwidth constraints,” IEEE WCNC, pp.1–5, 2007.
[8] J. Vermorel and M. Mohri, “Multi-armed ban-dit algorithms and empirical evaluation,” ECML, pp.437–448, 2005.
[9] H. Robbins, “Some aspects of the sequential de-sign of experiments,” Bull. Amer. Math. Soc., vol.58, no.5, pp.527–536, 1952.
[10] W. Chen, Y. Wang, and Y. Yuan, “Combinato-rial multi-armed bandit: General framework and applications,” ICML, pp.151–159, 2013.
[11] O. Chapelle and L. Li, “An empirical evaluation of Thompson sampling,” NIPS, pp.2249–2257, 2011.
[12] S. Agrawal and N. Goyal, “Analysis of Thompson sampling for the multi-armed bandit problem,” COLT, vol.23, pp.39.1–39.26, 2012.
[13] A. Gopalan, S. Mannor, and Y. Mansour, “Thompson sampling for complex online prob-lems,” ICML, pp.100–108, 2014.
[14] J. Komiyama, J. Honda, and H. Nakagawa, “Op-timal regret analysis of Thompson sampling in stochastic multi-armed bandit problem with mul-tiple plays,” ICML, pp.1152–1161, 2015.
[15] 飯塚 翔,川原 純,笠原正治,“k-out-of-Nルー ルによる協調センシングのためのマルコフ連鎖モ
ンテカルロ法を用いたパラメータ推定法,” 信学技
報,vol.117,no.204,pp.67–72,2017.
[16] M.J. Beal and Z. Ghahramani, “The variational bayesian EM algorithm for incomplete data: with application to scoring graphical model struc-tures,” Bayesian Stat., vol.7, pp.453–464, 2003. [17] S.J.D. Prince, Computer vision: models,
learn-ing, and inference, Cambridge University Press, 2012.
[18] M. Abramowitz and I.A. Stegun, Handbook of mathematical functions with formulas, graphs, and mathematical tables, Dover, 1972.