JP7438008B2 - 商品の販売数予測方法、装置及びプログラム、発注量決定方法、装置及びプログラム - Google Patents

Description

本発明は、店舗における商品の販売数を予測する技術及び発注量を決定する技術に関する。

従来、将来の商品販売数の予測や発注量の決定を行う手法が知られている（例えば特許文献１や非特許文献１参照）。特許文献１や非特許文献１に記載のものでは過去の販売実績や気象条件などに基づき状態空間モデルを用いて将来の販売数の予測を行っている。

周知のように状態空間モデルでは、観測可能な商品販売数が、長期的なデータ推移を示すトレンド成分や、季節的なデータ変動を示す季節調整成分、１週間など短期なデータ変動を示す短期調整成分などの複数の成分により推測されるものとしてモデル化したものである。状態空間モデルを用いた販売数の予測では、過去の商品販売数のほか前記各成分に関連する情報（例えば気象情報）を販売条件としてフィルタリング処理（１期先予測及びフィルタリング）を繰り返すことにより、現時点でのフィルタ分布を算出する。そして、当該フィルタ分布及び将来の販売条件を用いて１期先予測を繰り返すことにより、将来の任意の時刻（期）での予測フィルタを導出する。この予測フィルタは商品販売数の確率分布であり、これにより将来の任意の時刻（期）での商品販売数を予測する。

特許第６４３５８０３号公報

樋口知之，「予測にいかす統計モデリングの基本－ベイズ統計入門から応用まで」，講談社，２０１１年４月 佐藤・樋口，「ビッグデータ時代のマーケティング―ベイジアンモデリングの活用」，講談社，２０１３ 矢野・佐藤，「初期分布探索付き自己組織化状態空間モデルによる金融時系列解析の最前線：ｔ分布付き確率的ボラティリティ変動モデルへの応用」，ＦＳＡリサーチレビュー，２００６，Ｐ．１４３－１６６ Williams, C. K., & Rasmussen, C. E. (2006). Gaussian processes for machine learning (Vol. 2, No. 3, p. 4). Cambridge, MA: MIT Press CSATO, Lehel; OPPER, Manfred. Sparse on-line Gaussian processes. Neural computation, 2002, 14.3: 641-668. 中野・上野・中村，「Merging Particle Filterとその特性（特集 データマイニングと統計数理）」，統計数理，５６（２），２２５－２３４，２００８

しかし、従来の手法、特に特許文献１に記載されたようにフィルタリング処理においてカルマンフィルタを用いる手法では、システムモデルや観測モデルを表す式が複雑化すると何等かの近似を行う必要があり、推定精度に問題がある。すなわち、状態空間モデルを表す式が非線形化や非ガウス化すると、適用が困難になる。このため、従来は、商品販売数を理論的にモデル化できたとしても適用の困難さに起因して現実性に欠けるものとなり、結果としてより単純なモデル化が行われているだけであった。また単純なモデル化では、様々な商品の多様な販売挙動に対して適用することが困難であり、通年品や季節品といった挙動ごとに個別にモデルを構築する必要があった。

また、商品を販売する店舗では、将来の商品販売数の予測値に基づき当該商品の発注量を決定することが行われている。しかし、人手や上述のような単純なモデルによる予測値に基づいた発注量決定では、予測値の精度に起因し、欠品や過剰在庫が発生していた。このため、店舗別，商品別に適正な在庫量に基づいて発注量を決定することで、欠品や過剰在庫を抑制する必要があった。

本発明は上記事情に鑑みてなされたものであり、その目的とするところは、複雑なモデル化であっても計算コストが低く様々な商品に適用可能な販売数予測方法、装置及びプログラム、並びに発注量決定方法、装置及びプログラムを提供することにある。

上記目的を達成するために、本願発明は、商品販売数をコンピュータを用いて予測する販売数予測方法であって、前記コンピュータのデータ取得部が、入力された過去の販売実績数及び販売条件を取得する実績情報取得ステップと、前記コンピュータのフィルタ分布算出部が、状態変数を含むシステムモデル並びに前記状態変数及び外部変数を含む観測モデルからなり商品販売数を観測値とする非線形非ガウス型状態空間モデルを用い、過去の販売実績数及び販売条件に基づくフィルタリング処理により商品販売数のフィルタ分布を導出するフィルタ分布導出ステップとを備え、前記フィルタ分布導出ステップは、前記システムモデルに基づき前記状態変数を更新するステップと、オンラインガウス過程により前記外部変数及び前記オンラインガウス過程のモデルパラメータを更新するステップとを備えたことを特徴とする。

また、本願発明は、上記販売数予測方法において、前記コンピュータのデータ取得部が、入力された将来の販売条件を取得する将来情報取得ステップと、前記コンピュータの予測分布算出部が、前記非線形非ガウス型状態空間モデルを用い、前記フィルタ分布導出ステップにより導出されたフィルタ分布及び将来の販売条件に基づく予測処理により将来の商品販売数の予測分布を導出する予測分布導出ステップとを備え、前記予測分布導出ステップは、前記システムモデルに基づき前記状態変数を更新するステップと、オンラインガウス過程により前記外部変数を算出するステップとを備えたことを特徴とする。

また、本願発明は、商品の発注量をコンピュータを用いて決定する発注量決定方法であって、上記販売数予測方法における実績情報入力ステップ、フィルタ分布導出ステップ、及び予測分布導出ステップを備え、さらに、前記コンピュータの発注量算出部が、前記予測分布導出ステップにより算出された予測分布に基づき店舗における商品の必要在庫数を推定するとともに、必要在庫数に対する不足分を発注量として算出する発注量算出ステップを備えたことを特徴とする。

本発明によれば、非線形非ガウス型状態空間モデルとオンラインガウス過程の組み合わせによって販売数量予測モデルを構築することで様々な商品の多様な販売挙動に対して適用可能となり、商品ごとにモデルを構築することなく高精度の販売数予測が可能となる。また、本発明によれば、販売数予測を高精度で行うことができるので適正な在庫量の算出及びこれに基づく適正な発注量の決定が可能となり、これにより欠品や過剰在庫を抑制することができる。

販売数予測装置の機能ブロック図 販売数予測処理の概要図 システムモデルの式 気象データと特徴量の関係を説明する対応表 フィルタリング処理を説明するフローチャート 予測処理を説明するフローチャート 発注量決定処理の挙動を説明する図

本実施の形態に係る販売数予測装置は、非線形非ガウス型自己組織化状態空間モデルとオンラインガウス過程を組み合わせ、過去の販売実績数及び販売条件並びに将来の販売条件を入力として、店舗毎の将来の商品販売数の予測値を算出する。また、本実施の形態に係る販売数予測装置は、予測した商品販売数と店舗における在庫数に基づき当該店舗で必要となる将来の在庫数を算出するとともに、発注数量を算出する。

以下、本発明の一実施の形態に係る販売数予測装置について図１を参照して説明する。図１は販売数予測装置の機能ブロック図である。

販売数予測装置１は、図１に示すように、データ取得部１０と、過去実績データ記憶部２１と、将来条件データ記憶部２２と、在庫数記憶部２３と、フィルタ分布算出部３１と、予測分布算出部３２と、必要在庫数算出部３３と、発注量算出部３４と、分布データ記憶部４０と、必要在庫数及び発注量記憶部５０と、データ出力部６０とを備えている。

販売数予測装置１は、ＣＰＵ、メモリ、ＨＤＤ等を備えた周知のコンピュータからなる。販売数予測装置１の実装形態は不問である。例えば、販売数予測装置１は前記各部の一部又は全部を汎用のコンピュータにプログラムをインストールして実装してもよいし、ハードウェアとして実装してもよい。また、販売数予測装置１の各部は複数の装置に分散して実装してもよい。例えば、過去実績データ記憶部２１などの記憶手段をネットワーク上の他のコンピュータやストレージにより実現してもよい。

データ取得部１０は、キーボード等の所定の入力装置や所定の記憶媒体或いはネットワークを介して他の装置から各種データを取得して、取得したデータに対応する記憶部２１，２２，２３に当該データを記憶する。データ取得部１０によるデータの取得形態は不問である。

過去実績データ記憶部２１は、過去の販売実績数及び販売条件の時系列データを記憶する。過去の販売実績数及び販売条件は、フィルタ分布算出部３１での算出処理で用いられる。販売実績数は、店舗毎且つ商品毎の販売数である。販売条件は、商品販売に付帯するとともに商品販売数に影響を与えると考えられる情報である。販売条件としては、例えば、売価実績値，気象条件（気圧，降水量，気温，湿度，風速，日照時間，降雪量）などが挙げられる。

将来条件データ記憶部２２は、将来の販売条件を記憶する。将来の販売条件は、予測分布算出部３２で用いられる。将来の販売条件は、過去の平均的な価格、気象条件のみならず、価格設定を高めにした場合、安めにした場合、あるいは冷夏だった場合、猛暑だった場合といった想定し得る様々な状況に対応したデータを用意しておくことにより、多様な予測を行うことができる。

在庫数記憶部２３は、店舗における商品の在庫数を記憶する。在庫数は、発注量算出部３４で用いられる。在庫数は、店舗毎且つ商品毎の在庫数である。

フィルタ分布算出部３１は、後述するように、非線形非ガウス型自己組織化状態空間モデルとオンラインガウス過程を組み合わせて、過去の販売実績数及び販売条件に基づくフィルタリング処理により商品販売数のフィルタ分布を算出する。算出したフィルタ分布は、分布データ記憶部４０に記憶する。フィルタ分布算出部３１は、後述するように、フィルタ分布の算出過程において、オンラインガウス過程のモデルパラメータを逐次更新する。フィルタ分布算出部３１は、最終的なモデルパラメータを分布データ記憶部４０に記憶する。

予測分布算出部３２は、後述するように、非線形非ガウス型自己組織化状態空間モデルとオンラインガウス過程を組み合わせて、前記フィルタ分布算出部３１で算出されたフィルタ分布及び将来の販売条件に基づく予測処理により将来の商品販売数の予測分布を算出する。予測分布算出部３２は、販売数予測を所定の時間間隔（例えば、１日、１ヶ月、１年など）に対応した時間インデックスｉに対して複数行うことができ、これにより複数の時間幅に対する予測分布を算出することができる。算出した予測分布は、分布データ記憶部４０に記憶する。

分布データ記憶部４０は、フィルタ分布算出部３１により算出されたフィルタ分布及びモデルパラメータ並びに予測分布算出部３２で算出された予測分布を記憶する。フィルタ分布及び予測分布は、商品販売数の条件付き確率分布であり、本実施の形態では、粒子フィルタを用いるものとし、粒子の集合により各分布を表現・記憶する。

必要在庫数算出部３３は、商品販売数の予測分布に基づいて店舗毎且つ商品毎の必要在庫数を推定する。必要在庫数算出部３３は、予測分布算出部３２で複数の予測分布を算出した場合、各予測分布に対応する必要在庫量、換言すれば将来の条件情報に対応する必要在庫量を算出する。算出した必要在庫量は、必要在庫数及び発注量記憶部５０に記憶する。

発注量算出部３４は、店舗における商品の在庫数と前記必要在庫数算出部３３により算出された必要在庫数に基づき、当該必要在庫数に対する不足分を発注量として算出する。算出した発注量は、必要在庫数及び発注量記憶部５０に記憶する。

必要在庫数及び発注量記憶部５０は、店舗毎且つ商品毎の必要在庫数及び発注量を記憶する。

データ出力部６０は、分布データ記憶部４０や必要在庫数及び発注量記憶部５０に記憶されたデータを必要に応じて任意の所定の形式に変換して出力する。データ出力部６０の出力先は不問であり、例えばディスプレイなどの表示装置、記憶媒体、ネットワークを介した他の装置などが挙げられる。データ出力部６０は、予測分布やフィルタ分布を出力する際には、分布そのものをグラフ等により出力してもよいし、当該分布から商品販売数（スカラー値）の平均値やその他の統計値を算出して出力してもよい。

次に、本実施の形態に係る販売数予測装置における販売数予測処理の概要について図２を参照して説明する。図２は販売数予測処理の概要図である。

図２に示すように、本実施の形態では、状態変数（１０４）及び外部変数（１０５）並びにオンラインガウス過程パラメータ（１０６）からなるパラメータデータと、過去の販売数量，気象，価格等から成る過去実績データ（１０１）を用いて状態変数（１０４）及び外部変数（１０５）並びにオンラインガウス過程パラメータ（１０６）を最新の状態に更新するフィルタリングステップ（１０７）と、更新された状態変数（１０４）及び外部変数（１０５）並びにパラメータデータ（１０６）と気象および価格等からなる将来条件データ（１０８）を用いて将来の状態変数（１１１）及び外部変数（１１２）を更新し、販売数量予測値（１１４）を求める予測ステップ（１１３）と、販売数量予測値（１１４）に基づき、必要在庫量予測値（１１６）を求める必要在庫決定処理（１１５）と、必要在庫量予測値（１１６）と新規実績データ（１１８）を用い、発注量（１１９）を求める発注量決定処理（１１７）とを備えている。

次に、本実施の形態に係る販売数予測装置における販売数予測処理の詳細なアルゴリズムについて説明する。

［販売数量のモデル化］

ある店舗、ある商品のｉ日（ｉ＝１，２，３，…）における販売数量をｙ（ｉ）とすると、このような少数の販売数量の分布は以下のようなポアソン分布に従うと考えられる。

ここでλ（ｉ）はポアソン分布のパラメータを表し、日々変動するものとする。λ（ｉ）はｉ日における販売数量の分布の平均および分散に等しい。またｋ（ｋ＝１，２，３，…）は販売数量を表す。ここでλ（ｉ）は以下のような構造を持つと考える。

ここでｔ（ｉ）はトレンド成分、ｓ（ｉ）は季節成分、ｗ（ｉ）は曜日成分、β（ｉ）は日付成分、ｒ（ｉ）は価格成分をそれぞれ表す。

［各成分のモデル化］

＜トレンド成分＞
トレンド成分は販売数量のベースとなる成分であり、緩やかに変化する。以下のようなランダムウォークモデルで表す。

ｖ_ｔ（ｉ）は平均０、標準偏差σ_ｔ（ｉ）のガウス分布に従うノイズ項である。

＜季節成分＞
季節成分は季節とともに変動する成分である。また気象条件による変化もこの成分に含める。季節成分は以下のようなカーネル回帰でモデル化する。

後述する通り、ｋ（，）はカーネル関数、ｕ_ｉ，ｕ_ｊは日付，気象データから作成した特徴量ベクトルであり、ｍおよびａ_ｊはオンラインガウス過程により決定される。

ｖ_ｗ（ｉ）は平均０、標準偏差σ_ｗ（ｉ）のガウス分布に従うノイズ項である。

＜日付成分＞
日付成分は販売数量の各月の日付（１日～３１日）に対する依存分である。日付成分はｊ日（ｊ＝１，…，３１）の成分をβ_ｊ（ｉ）とし、以下のモデルで表す。

ｖ_βｊ（ｉ）は平均０、標準偏差σ_β（ｉ）のガウス分布に従うノイズ項、ｄ（ｉ）は日付インデックスｉの日付、δ（）はクロネッカーのデルタ関数である。なお、「ｖ_βｊ（ｉ）」の標記において、ｊはβの添え字である。

＜価格成分＞
価格成分は販売数量の価格への依存分である。価格成分をモデル化するために参照価格の概念を導入する（非特許文献２参照）。参照価格は、消費者が商品の高値、安値を判断する基準の価格であり、実売価格が参照価格よりも高ければ割高感を、実売価格が参照価格よりも安ければ割安感を消費者に与える。参照価格をｐ_ｒｅｆ（ｉ）、実売価格をｐ_ａｃｔ（ｉ）と表すと、ｐ_ｒｅｆ（ｉ）はｐ_ａｃｔ（ｉ）を１次指数平滑化した値として表される。平滑化係数をγ（ｉ）（０＜γ（ｉ）＜１）とすると、

(8)
価格成分ｒ（ｉ）はｐ_ａｃｔ（ｉ）とｐ_ｒｅｆ（ｉ）の差額に価格弾性係数ξ（ｉ）を掛けたものとして表される。

従って、γ（ｉ），ξ（ｉ）の時間変動を表現すればｒ（ｉ）の時間変動を表現できるが、０＜γ（ｉ）＜１、０＜ξ（ｉ）であることが保証されている必要があるため、以下のような形でγ（ｉ），ξ（ｉ）の時間変動を考える。

ｖ_α（ｉ），ｖ_ｑ（ｉ）はそれぞれ平均０、標準偏差σ_α（ｉ），σ_ｑ（ｉ）のガウス分布に従うノイズ項である。

［システムモデルおよび観測モデル］

直接観測できないトレンド成分，曜日成分，日付成分に関する成分（β_ｊ（ｉ），ｊ＝１，２，…，３１），価格成分に関連する成分（α（ｉ）およびｑ（ｉ））をまとめたベクトルｘ（ｉ）＝［ｔ（ｉ），ｗ（ｉ），ｗ（ｉ－１），…，ｗ（ｉ－５），β_１（ｉ），…， β_３１（ｉ）］^Ｔを状態ベクトルとし、これらの時間変動を表す式（３），式（５），式（６）をまとめたものがシステムモデルである。この状態ベクトルｘ（ｉ）に含まれる各状態変数は、特許請求の範囲の状態変数に相当する。

今、パラメータの時間変動による予測精度低下を回避するためにσ_ｔ（ｉ），σ_ｗ（ｉ），σ_β（ｉ），σ_α（ｉ），σ_ｑ（ｉ）の時間変動を以下のように考える。

ｕ_ｔ（ｉ），ｕ_ｗ（ｉ），ｕ_β（ｉ），ｕ_α（ｉ），ｕ_ｑ（ｉ）はそれぞれ［－Δｔ， Δｔ］，［－Δｗ， Δｗ］，［－Δβ， Δβ］，［－Δα， Δα］，［－Δｑ， Δｑ］の範囲の一様分布である。式（１４）～式（１８）は各成分の時間変動を決定するパラメータの時間変動であるが、前述の状態ベクトルｘ（ｉ）にこれらパラメータベクトルθ（ｉ）＝［α（ｉ），ｑ（ｉ），σ_ｔ（ｉ），σ_ｗ（ｉ），σ_β（ｉ），σ_α（ｉ），σ_ｑ（ｉ）］^Ｔを追加して拡大した状態ベクトルとして考えると、時変なパラメータの時間変動を示す拡張したシステムモデルを考えることができる。

また拡大した状態変数から出力である販売数量が生成される過程を記述したものが観測モデルであり、本実施の形態では販売数量実績値が式（１）および式（２）で示されるポアソン分布に従うことを示している。

以下、拡大した状態ベクトルをｚ（ｉ）＝［ｘ（ｉ）^Ｔ，θ（ｉ）^Ｔ］^Ｔとし、システムモデルおよび観測モデルをそれぞれ以下のように表現する。

ただし、

である。システムモデルの式を図３に示す。

なお、拡大した状態変数によってモデルのパラメータの時間変動を扱う手法は自己組織化状態空間モデルとして知られている（非特許文献３参照）。消費者の購買動向の変化などでモデルパラメータが変化した場合でも、常に良好な予測値が得られる状態を維持できることが期待できる。

本実施の形態では、季節成分については、状態空間モデルの外で別途推定（オンラインガウス過程）することを特徴としている。なお、季節成分に係る変数ｓ（ｉ）は、特許請求の範囲の外部変数に相当し、それ以外の成分に係る状態変数は、特許請求の範囲の状態変数に相当する。ここで、「外部変数」とは、状態空間モデルの観測モデルに含まれる変数であるがシステムモデルには含まれない変数を意味するものとする。

［オンラインガウス過程］

商品の年間の販売数量の推移に対して、従来は式（５）に示した季節調整モデルやＡＲモデル，ＡＲＩＭＡモデルといった時系列モデルが適用されてきたが、これらのモデルは気象条件などの外部要因を出力に反映させることが困難である。本実施の形態ではこの点を解決するために、販売数量における季節成分の推定にオンラインガウス過程を用いる。

本来のガウス過程は過去に得られた特徴量ｕ（ｉ）（ｉ＝１，…，Ｍ）とそれに対応する出力ｓ（ｉ）を用い、新たに得られた特徴量ｕ（Ｍ＋１）に対する出力の平均ｍ（ｕ（Ｍ＋１））と分散σ^２（ｕ（Ｍ＋１））を以下のように推定する（非特許文献４参照）。なお、「特徴量ｕ（ｉ）」はベクトルである点に留意されたい。

ここで、カーネル関数をｋ（，）とすると、Ｋはｉ次元をｋ（ｕ（ｉ），ｕ（Ｍ＋１））とするベクトル、Ｃ_Ｍは（ｉ，ｊ）要素をｋ（ｕ（ｉ），ｕ（ｊ））とする行列。ｓはｉ次元をｓ（ｉ）とするベクトルである。ガウス過程は非線形な対応関係を柔軟に表現する優れた手法であるが、計算量がｏ（Ｍ^３）に比例するためにデータ数Ｍが大きくなると計算の負荷が無視できなくなり、オンライン学習のような逐次データを取り込みながらモデルを更新するといった用途には適さない。また、消費者の嗜好や気象条件に対する反応が時間と共に変化する場合、学習開始時からの特徴量をすべて保持するガウス過程では新たな入力に対する予測精度が低下する恐れがある。本実施の形態ではこれらの問題をオンラインガウス過程を用いることで回避した。

オンラインガウス過程（非特許文献５参照）は特徴量ｕ（ｉ）（ｉ＝１，…，Ｍ）を逐次与え、その都度モデルを更新していくが、与えた特徴量を全て使用するのではなく、重要な特徴量のみを基底ベクトルセットＢＶに取り込み、ＢＶに取り込んだ特徴量のみを使用して式（２２），式（２３）を近似する。ＢＶに取り込む基底ベクトルの数には上限があり、取り込んだ基底ベクトルの数が上限を超えた場合は、最も出力に対する影響の小さい基底ベクトルがＢＶから取り除かれる。モデルの更新と基底ベクトルの取捨を逐次繰り返し、常に予測精度を最適に保持できるため、先に挙げた消費者の嗜好や気象条件に対する反応が変化するような場合においても、常に正確な予測を行えるようモデルの状態を保つことが期待できる。

以下では、ｉ番目までの特徴量を与えて学習したオンラインガウス過程ＧＰ_ｉに未知の特徴量ｕを与えて予測を行う場合を、

のように表し、ｕ（ｉ＋１），ｓ（ｉ＋１）を与えてＧＰ_ｉをＧＰ_ｉ＋１へ更新する場合を、

と表現する。

本実施の形態では、季節成分の予測にオンラインガウス過程を適用する。特徴量としては日付データおよび気象データを用いる。

日付データは以下の換算値を用いる。ｄ’（ｉ）はその年の１月１日からの経過日数を表す。

気象データは、例えば、平均気圧、降雨量、平均気温、平均湿度、平均風速、日照時間、降雪量を含む。気象データは各日の気象データ項目から、それぞれの過去数年の平均値を引き、その値をそれぞれの標準偏差で割って正規化し、更に６日前～当日の値の平均をとったものを使用する。

特徴ベクトルの各次元への対応を図４に示す。これらをまとめたベクトルｕ（ｉ）＝［ｕ_１（ｉ），ｕ_２（ｉ），ｕ_３（ｉ），ｕ_４（ｉ），ｕ_５（ｉ），ｕ_６（ｉ），ｕ_７（ｉ），ｕ_８（ｉ）］^Ｔが季節成分の特徴量である。

なお気象成分ｓ（ｉ）は実際には観測されないので、モデルの学習の際に必要となるｓ（ｉ）の値は後述の通りフィルタリングステップの中で得られる季節成分の事後推定値を用いる。

［フィルタリングステップ］

フィルタリング処理は、図１に示すフィルタ分布算出部３１により実施される。

ここでは、ｉ＝１，…，Ｍ（Ｍ＝３６５，２×３６５，…）における販売数量，単価，気象データが既に得られており、ｉ＝１からｉ＝Ｍまでフィルタリングを行うことで状態ベクトルｚ（ｉ），オンラインガウス過程モデルＧＰ_ｉを逐次求めるという前提で説明を行う。

＜状態変数，パラメータ，オンラインガウス過程モデルの初期化＞
ｉ＝１からフィルタリングを開始するにあたり、状態変数，パラメータ，オンラインガウス過程モデルの初期値ｘ（０），θ（０），ＧＰ_０が必要となる。これらはそれぞれ以下のように求める。

状態ベクトルｘ（ｉ）は販売数量の挙動をトレンド成分ｔ（ｉ），曜日成分ｗ（ｉ），日付成分β_ｊ（ｉ）（ｊ＝１，…，３１）に分解したものなので、各成分の初期値については対数をとった販売数量実績値のデータに平滑化処理や差分等の処理を行って求める。例えばトレンド成分ｔ（０）であれば販売数量実績値の対数をとった後、ｉ＝１，…，Ｍの期間の平均をとって求めることができる。また、販売数量データの対数をとった後に長さ７の移動平均をとれば曜日成分を除去することができるので、移動平均前のデータから移動平均後のデータを引けば、残差データには曜日成分のみが残ることになる。残差データの値に対して曜日毎の平均をとることで曜日成分の初期値である｛ｗ（０），…，ｗ（－６）｝を得ることができる。同様に、日付成分については長さ３１の移動平均データ、および長さ７の移動平均データを元のデータから引くことで得ることができる。

θ（０）についてはフィルタリング期間における尤度を最大にすることで最適なθ（０）を決定する手法（非特許文献３参照）などが考案されており、これらの手法を用いて求めることができる。

なおパーティクルフィルタ内においては、拡大した状態ベクトルｚ（ｉ）＝｛ｘ（ｉ）^Ｔ，θ（ｉ）^Ｔ｝^ＴはＮ個の粒子で表現されている。以下、フィルタ分布（事後分布）に従う粒子を

予測分布に従う粒子を

のように表す。先に挙げた初期化手法では

が得られることとなる。

ＧＰ_０の学習を行うにあたっては、まず対数をとった販売数量データからトレンド成分，曜日成分，日付成分を引いた残差データｓ_０（ｉ）（ｉ＝１，…，Ｍ）を求める。このデータと前述の季節成分特徴量との組

をオンラインガウス過程に逐次与え、ＧＰ_０の学習を行うことができる。

＜フィルタリング＞
本実施の形態では、パーティクルフィルタとオンラインガウス過程を組み合わせたフィルタリング処理を行う。フィルタリングステップのフローチャートを図５に示す。

現時点をｉとし、既に（ｉ－１）における拡大した状態変数のフィルタ分布に従う粒子

とオンラインガウス過程のモデルＧＰ_ｉ－１が得られているとする（図５のステップＳ１０２，Ｓ１０３）。フィルタリングは以下の手順で行う。

（ステップ１）
ＧＰ_ｉ－１により季節成分平均値ｓ^～（ｉ）を求める（図５のステップＳ１０６）。なお、「ｓ^～」における記号チルダは記号ｓの上方に記載されていることを表している。

（ステップ２）
それぞれＮ（０，σ_ｔ（ｉ－１）），Ｎ（０，σ_ｗ（ｉ－１）），Ｎ（０，σ_β（ｉ－１）），Ｎ（０，σ_α（ｉ－１）），Ｎ（０，σ_ｑ（ｉ－１）），Ｕ（－Δｔ，Δｔ），Ｕ（－Δｗ，Δｗ），Ｕ（－Δβ，Δβ），Ｕ（－Δα，Δα），Ｕ（－Δｑ，Δｑ）に従う粒子をＮ個ずつ発生させて、これを

とする。これと、

をシステムモデルに与え、 １期先予測分布に従う粒子

を得る（図５のステップＳ１０７）。

（ステップ３）
１期先予測の粒子一つ一つについて尤度

を求める（図５のステップＳ１０８）。

（ステップ４）
得られた尤度を以下のように正規化する（図５のステップＳ１０９）。

（ステップ５）
１期先予測の粒子

を

の重みでリサンプリングし、得られた粒子を

と表す（図５のステップＳ１１０）。これは状態ベクトルｚ（ｉ）のフィルタ分布（ｙ（ｉ）を観測した上でのｚ（ｉ）の事後分布）を表している。

なお、リサンプリングについてはＭｅｒｇｉｎｇ Ｐａｒｔｉｃｌｅ Ｆｉｌｔｅｒ（非特許文献６参照）等の手法を適用し、高精度化を図ることができる。

（ステップ６）

より求めたトレンド成分

，曜日成分

，日付成分

，価格成分

より、季節成分の期待値ｓ^＾（ｉ）を求め、ＧＰ_ｉの学習を行う（図５のステップＳ１１１）。なお、「ｓ^＾」における記号ハットは記号ｓの上方に記載されていることを表している。

Ｅ｛｝は期待値をとることを表す。

（ステップ７）
ｉ＝Ｍまでステップ１～６を繰り返す。

（ステップ８）
ｉ＞Ｍとなったら、

及びＧＰ_Ｍでパラメータデータを更新し、パラメータデータを記載したファイル，データベース等に書出す（図５のステップＳ１１３）。

なお、ステップ２～５は通常のパーティクルフィルタによるフィルタリングであり、フィルタリングのループにステップ１および６の処理を含めたことが独自の手法である。

［予測ステップ］

予測処理は、図１に示す予測分布算出部３２により実施される。

ここでは、ｉ＝Ｍまでのフィルタリングが既に完了して

及びＧＰ_Ｍが得られており、ｉ＝Ｍ＋１からｉ＝Ｌ（Ｍ＜Ｌ）までの予測を行うことで将来の販売数量の予測値ｙ^＾（ｉ）（ｉ＝Ｍ＋１，…，Ｌ）を逐次求めるという前提で説明を行う。なお、「ｙ^＾」における記号ハットは記号ｙの上方に記載されていることを表している。

＜将来の価格データおよび気象データ＞
予測においては、ｉ＝Ｍ＋１，…，Ｌにおける単価，気象データが必要となるが、未来の単価や気象データについては事前に入手することはできない。従って販売計画値や平均気象データ等の値をあらかじめ決定しておく必要がある。複数の施策案による価格設定パターンや、温暖／寒冷といった気象パターンを複数与え、条件による販売数量の変化を得ることもできるので、目的に応じてこれらの値を事前に決定しておく必要がある。

＜予測＞
本実施の形態では、パーティクルフィルタとオンラインガウス過程を組み合わせた予測処理を説明する。予測ステップのフローチャートを図６に示す。

現時点をｉ（ｉ＝Ｍ＋１，Ｍ＋２，…）とし、既に（ｉ－１）における拡大した状態変数の予測分布に従う粒子

と、ｉ＝Ｍにおけるオンラインガウス過程のモデルＧＰ_Ｍが得られているとする（図６のステップＳ２０２，Ｓ２０３）。予測は以下の手順で行う。

（ステップ１）
ＧＰ_Ｍにより季節成分平均値ｓ^～（ｉ）を求める（図６のステップＳ２０６）。なお、「ｓ^～」における記号チルダは記号ｓの上方に記載されていることを表している。

（ステップ２）
それぞれＮ（０，σ_ｔ（ｉ－１）），Ｎ（０，σ_ｗ（ｉ－１）），Ｎ（０，σ_β（ｉ－１）），Ｎ（０，σ_α（ｉ－１）），Ｎ（０，σ_ｑ（ｉ－１）），Ｕ（－Δｔ，Δｔ），Ｕ（－Δｗ，Δｗ），Ｕ（－Δβ，Δβ），Ｕ（－Δα，Δα），Ｕ（－Δｑ，Δｑ）に従う粒子をＮ個ずつ発生させて、これを

とする。これと、

をシステムモデルに与え、予測分布に従う粒子

を得る（図６のステップＳ２０７）。

（ステップ３）
季節成分平均値ｓ^～（ｉ）と、

より求めたトレンド成分

，曜日成分

，日付成分

，価格成分

より販売数量の予測値を求める。ただし、式（１）に示すように販売数量はポアソン分布に従うとしているので、ここではポアソン分布のパラメータの分布を近似する粒子

を求める（図６のステップＳ２０８）。

（ステップ４）
ｉ＝Ｌまでステップ１～３を繰り返す。

（ステップ５）
ｉ＞Ｌとなったら、

を、予測値を記載したファイル，データベース等に書出す（図６のステップＳ２１０）。

なお、ステップ２および３は通常のパーティクルフィルタによる予測であり、予測のループにステップ１の処理を含めたことが独自の手法である。

［必要在庫決定処理］

必要在庫決定処理は、図１に示す必要在庫数算出部３３により実施される。

上述の予測ステップにより

が得られた後，これらのデータに基づき日々の必要在庫量を予測する。

の期待値をλ^￣ _ｉ｜Ｍと表すと、ｉ日における販売数量予測値ｙ^～（ｉ）はλ＝λ^￣ _ｉ｜Ｍのポアソン分布

に従う。なお、「λ^￣」における記号オーバーラインは記号λの上方に記載されていることを表している。また、「ｙ^～」における記号チルダは記号ｙの上方に記載されていることを表している。

例えば累積確率Ｐ＝９９％を与えるｙ^～（ｉ）の値を累積確率関数等から求めると、ｉ日において９９％の確率で欠品が生じない販売数量ｙ_Ｐ（ｉ）を求めることができる。Ｐの値は目的に応じて変更可能である。店頭において最小限確保しておきたい数量がある場合、この値をτとして必要在庫量ｙ_ｕ（ｉ）を以下のように決定することができる。

ｙ_ｕ（ｉ）は、ｉにおいて最小限確保しておきたい数量τを確率Ｐで保証する在庫数量である。

上記の方法は日別に必要在庫量を求める方法であるが、通常はｍ日毎に店舗へ商品が送られることが多い。このような場合、ｍ日間での必要在庫量を考える必要がある。ここでＰｏｉｓｓ（ｘ_１＝ｋ｜λ_１）に従う確率変数ｘ_１とＰｏｉｓｓ（ｘ_２＝ｋ｜λ_２）に従う確率変数ｘ_２との和ｘ_１＋ｘ_２はＰｏｉｓｓ（ｘ_１＋ｘ_２＝ｋ｜λ_１＋λ_２）に従うという性質から、ｉ～ｉ＋ｍ－１のｍ日間の平均販売数量λをΛ^￣ _ｉ｜Ｍと表せば、

で求められるので、必要在庫量はλ＝Λ^￣ _ｉ｜Ｍに従うポアソン分布の累積確率関数から求めたｍ日間の販売数量Ｙ_Ｐ（ｉ）として求めることができる。この場合の必要在庫量をＹ_ｕ（ｉ）、最小限確保したい数量をτとすると、

Ｙ_ｕ（ｉ）は、ｉ～ｉ＋ｍ－１のｍ日間において最小限確保しておきたい数量τを確率Ｐで保証する在庫数量である。

以上、日別で必要在庫量を管理する場合はｙ_ｕ（ｉ）を、ｍ日間の必要在庫量を管理する場合はＹ_ｕ（ｉ）をそれぞれの場合の必要在庫量予測値とする。

［発注量決定処理］

発注量決定処理は、図１に示す発注量算出部３４により実施される。

発注量決定処理は必要在庫量予測値と店舗における販売数量実績値を比較し、次回必要在庫量予測値と店舗在庫の差分を発注量とする。次回必要在庫量予測値＜（今回必要在庫量－今回販売実績値（ｍ日間集計））の場合は次回発注を行わない。図７に示す。なお、図中のＹ（ｉ）は、

とする。また、ここでは発注分は翌日に店舗へ届けられるものとする。

本実施の形態に係る販売数予測装置によれば、非線形非ガウス型状態空間モデルとオンラインガウス過程の組み合わせによって販売数量予測モデルを構築しているので、様々な商品の多様な販売挙動に対して適用可能となり、商品ごとにモデルを構築する必要がない。また、本実施の形態に係る販売数予測装置によれば、店舗別，商品別に適正な在庫量に基づく出荷量を決定するので、顧客に対する欠品等が生ずることがなく、また店舗にとって過剰在庫にならない商品の在庫量，商品の出荷量あるいは供給量の決定を低コストで提供することができる。

また、本実施の形態に係る販売数予測装置によれば、非線形非ガウス型自己組織化状態空間モデルとオンラインガウス過程の組み合わせによって販売数量予測モデルを構築しているので、モデルパラメータを自ら調整し、常に高い予測精度を維持できるとともに、メンテナンス負担を小さく抑えることができる。

また、販売数量予測モデルにより推定される販売数量予測分布は、最小の時間間隔（例えば日別）に対応した時間インデックスｉについて推定されるが、時間インデックスｉに対する販売数量予測分布から、時間インデックスｉよりも大きな任意の期間（例えば週や月など）に対する販売数量予測分布を統計的な手法によって求めることで、任意の発注間隔に対応することができる。

また、上記実施の形態では、販売数量予測モデルにより推定される販売数量予測分布は店舗別・商品別に推定しているが、更に、複数の店舗または複数の商品またはその双方によるグループ化を行い、店舗別・商品別の販売数量予測分布から前述のグループに対する販売数量予測分布を統計的な手法によって求めることで、任意の店舗群・商品群に対応することができる。

また、本実施の形態に係る販売数予測装置によれば、日別，週別，月別等の任意の期間の必要在庫量に対する不足分を手動または自動により発注することで、過剰在庫および欠品を同時に抑制することができる。

以上、本発明の一実施の形態について詳述したが、本発明は上記実施の形態に限定されるものではなく、本発明の主旨を逸脱しない範囲において、種々の改良や変更をしてもよい。

例えば、上記実施の形態では、粒子フィルタを用いて状態変数を推定していたが、カルマンフィルタやその拡張手法など他の非線形なフィルタ手法を適用することもできる。ただし、モデルが適用条件を満たしていることが条件となる。もっとも、予測精度の向上の観点から非線形非ガウス型状態空間モデルを用いることに鑑みると、計算量等の観点から粒子フィルタを用いることが好適である。

上記実施の形態では、状態空間モデルに含まれる季節成分をオンラインガウス過程による機械学習処理の対象とした例について説明したが、他の成分を対象としてもよい。ただし、オンラインガウス過程の対象にできる条件は、入力ｘ_ｎ（ｎ＝１，２，…）に対する出力ｙ_ｎの同時分布が常にガウス分布となるデータである。なお、季節商品などでは気象条件によってオンシーズン時の販売数のピークの高さが変わっており、従ってモデルも与えた気象条件に応じた変化を表現できる必要があるものの、この対応関係を人手でモデル化するのが困難であった。本実施の形態では、この問題を解決するためにガウス過程のような機械学習手法を導入している。したがって、同様に人手によるモデル化が困難な成分についてオンラインガウス過程の対象とすると好適である。

また、上記実施の形態では、トレンド成分，季節成分，曜日成分，価格成分をモデル化して状態空間モデルを構成していたが、任意の成分をモデル化し且つ各モデルを任意に組み合わせて状態空間モデルを構成することができる。例えば、山積みやチラシ、特売日などのデータがあるのであれば、それらをモデルに取り込むことが可能である。また、上記実施の形態における各成分のモデル化は一例であり、他のモデルを採用してもよい。例えば、上記実施の形態ではトレンド成分としてランダムウォークを採用したが、２階差分モデルを用いてもよい。

また、上記実施の形態では、自己組織化状態空間モデルを用いることにより、消費者の購買動向の変化などでモデルパラメータが変化した場合でも、常に良好な予測値が得られる状態を維持できるように構成しているが、自己組織化ではない状態空間モデルであってもよい。

また、上記実施の形態では、販売数量を様々な挙動の成分に分解するという方向でモデル化したが、消費者の消費行動自体をモデル化できるのであれば、そちらの方向からの状態空間モデル構築を行うこともできる。

また、上記実施の形態では、過去の販売数量実績値の分布がポアソン分布に比較的よく適合したため観測モデルにポアソン分布を用いたが、異なる販売環境（例えばネット販売など）で他の分布の方がよく適合するという結果が得られれば、分布の型を変えたモデルを作成してもよい。

１…販売数予測装置
１０…データ取得部
２１…過去実績データ記憶部
２２…将来条件データ記憶部
２３…在庫数記憶部
３１…フィルタ分布算出部
３２…予測分布算出部
３３…必要在庫数算出部
３４…発注量算出部
４０…分布データ記憶部
５０…必要在庫数及び発注量記憶部
６０…データ出力部

Claims (12)

1. 商品販売数をコンピュータを用いて予測する販売数予測方法であって、
   前記コンピュータのデータ取得部が、入力された過去の販売実績数及び販売条件を取得する実績情報取得ステップと、
   前記コンピュータのフィルタ分布算出部が、状態変数を含むシステムモデル並びに前記状態変数及び外部変数を含む観測モデルからなり商品販売数を観測値とする非線形非ガウス型状態空間モデルを用い、過去の販売実績数及び販売条件に基づくフィルタリング処理により商品販売数のフィルタ分布を導出するフィルタ分布導出ステップとを備え、
   前記フィルタ分布導出ステップは、前記システムモデルに基づき前記状態変数を更新するステップと、オンラインガウス過程により前記外部変数及び前記オンラインガウス過程のモデルパラメータを更新するステップとを備えた
   ことを特徴とする販売数予測方法。
2. 前記コンピュータのデータ取得部が、入力された将来の販売条件を取得する将来情報取得ステップと、
   前記コンピュータの予測分布算出部が、前記非線形非ガウス型状態空間モデルを用い、前記フィルタ分布導出ステップにより導出されたフィルタ分布及び将来の販売条件に基づく予測処理により将来の商品販売数の予測分布を導出する予測分布導出ステップとを備え、
   前記予測分布導出ステップは、前記システムモデルに基づき前記状態変数を更新するステップと、オンラインガウス過程により前記外部変数を算出するステップとを備えた
   ことを特徴とする請求項１記載の販売数予測方法。
3. 前記フィルタ分布導出ステップでは、粒子フィルタリングにより前記システムモデルに基づき前記状態変数を更新する
   ことを特徴とする請求項１又は２何れか１項記載の販売数予測方法。
4. 前記フィルタ分布導出ステップにおける粒子フィルタリングでは、前記状態変数及び前記外部変数に基づき尤度を算出し、算出した尤度で粒子を重み付けするステップを含む
   ことを特徴とする請求項３記載の販売数予測方法。
5. 前記予測分布導出ステップでは、粒子フィルタリングにより前記システムモデルに基づき前記状態変数を更新する
   ことを特徴とする請求項２記載の販売数予測方法。
6. 前記非線形非ガウス型状態空間モデルは非線形非ガウス型自己組織化状態空間モデルである
   ことを特徴とする請求項１乃至５何れか１項記載の販売数予測方法。
7. 前記販売条件は気象情報を含み、
   前記非線形非ガウス型状態空間モデルは前記外部変数として季節成分を含み、
   前記オンラインガウス過程では前記気象情報に基づき算出した特徴量を用いて前記季節成分及び前記モデルパラメータを更新する
   ことを特徴とする請求項１乃至６何れか１項記載の販売数予測方法。
8. 商品の発注量をコンピュータを用いて決定する発注量決定方法であって、
   請求項２記載の実績情報取得ステップ、フィルタ分布導出ステップ、及び予測分布導出ステップを備え、さらに、
   前記コンピュータの発注量算出部が、前記予測分布導出ステップにより算出された予測分布に基づき店舗における商品の必要在庫数を推定するとともに、必要在庫数に対する不足分を発注量として算出する発注量算出ステップを備えた
   ことを特徴とする発注量決定方法。
9. 商品販売数を予測する販売数予測装置であって、
   所定の記憶手段に記憶されている過去の販売実績数及び販売条件を入力として、状態変数を含むシステムモデル並びに前記状態変数及び外部変数を含む観測モデルからなり商品販売数を観測値とする非線形非ガウス型状態空間モデルを用い、過去の販売実績数及び販売条件に基づくフィルタリング処理により商品販売数のフィルタ分布を導出するフィルタ分布導出手段を備え、
   前記フィルタ分布導出手段は、前記システムモデルに基づき前記状態変数を更新するとともに、オンラインガウス過程により前記外部変数及び前記オンラインガウス過程のモデルパラメータを更新する
   ことを特徴とする販売数予測装置。
10. 所定の記憶手段に記憶されている将来の販売条件を入力として、前記非線形非ガウス型状態空間モデルを用い、前記フィルタ分布導出手段により導出されたフィルタ分布及び将来の販売条件に基づく予測処理により将来の商品販売数の予測分布を導出する予測分布導出手段を備え、
    前記予測分布導出手段は、前記システムモデルに基づき前記状態変数を更新するとともに、オンラインガウス過程により前記外部変数を算出する
    ことを特徴とする請求項９記載の販売数予測装置。
11. 請求項９又は１０記載のフィルタ分布導出手段及び予測分布導出手段を備え、さらに、
    前記予測分布導出手段により算出された予測分布に基づき店舗における商品の必要在庫数を推定するとともに、必要在庫数に対する不足分を発注量として算出する発注量算出手段を備えた
    ことを特徴とする発注量決定装置。
12. コンピュータを請求項９又は１０記載の販売数予測装置或いは請求項１１記載の発注量決定装置として機能させることを特徴とするプログラム。
    

Images