"測位用信号を受信することと、前記測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことと、を含む位置算出方法。過去の前記位置算出演算の結果に基づいて前記変動の基準を更新することと、前記基準を用いて前記変動を算出することと、を更に含む、請求項1に記載の位置算出方法。前記位置算出演算は、少なくとも前記演算結果の変動を確率変数とした正規分布モデルに基づいて定められた演算であり、過去の前記位置算出演算の結果に基づいて前記正規分布モデルを更新することを更に含む、請求項1又は2に記載の位置算出方法。前記受信信号に基づいて擬似距離を測定することを更に含み、前記位置算出演算は、前記演算結果の変動を確率変数とした正規分布モデルと、前記擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた演算である、請求項3に記載の位置算出方法。前記擬似距離に含まれる誤差を確率変数とした確率分布モデルは、正規分布モデル又は正規混合分布モデルである、請求項4に記載の位置算出方法。前記位置算出演算は、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する演算であり、前記正規分布モデルは、当該演算結果の各要素の変動を確率変数としたモデルである、請求項3~5の何れか一項に記載の位置算出方法。前記位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、前記受信信号を用いた所定の収束演算を行って算出することを更に含み、前記位置算出演算は、過去に前記収束演算によって算出された前記所定のパラメーターの値を用いる演算であることを含む、請求項1~6の何れか一項に記載の位置算出方法。測位用信号を受信する受信部と、前記受信部による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部と、を備えた位置算出装置。従来の位置算出演算では、位置算出演算を行う時刻以下、演算時刻と称す。毎に、最小二乗法を利用するなどして位置算出演算を1回1回個別に行うのが通常である。従って、例えば、測定した擬似距離に瞬間的に大きな誤差が混入した場合、そのときの位置算出演算の精度が急激に低下するという欠点があったが、以降の位置算出演算には影響しないという利点もあった。但し、瞬間的とは言え、真の位置から大きく離れた位置が演算結果として得られる、いわゆる位置飛びが発生することは、場合によって大きな問題となり得た。本発明は上述した課題に鑑みて為されたものであり、その目的とするところは、位置算出の正確性を向上させるための新しい位置算出手法を提案することにある。以上の課題を解決するための第1の形態は、測位用信号を受信することと、前記測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことと、を含む位置算出方法である。また、他の形態として、測位用信号を受信する受信部と、前記受信部による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部と、を備えた位置算出装置を構成してもよい。この第1の形態等によれば、測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う。これにより、各演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行うことが可能となる。この場合、測定量に瞬間的に大きな誤差が混入したとしても、位置算出精度が低下することを防止できる。また、第2の形態として、第1の形態の位置算出方法において、過去の前記位置算出演算の結果に基づいて前記変動の基準を更新することと、前記基準を用いて前記変動を算出することと、を更に含む、位置算出方法を構成することとしてもよい。この第2の形態によれば、過去の位置算出演算の結果に基づいて演算結果の変動の基準を更新することで、過去の位置算出演算の結果を今回の位置算出演算の結果に反映させることが可能となる。更新した基準を用いて演算結果の変動を算出することで、位置算出の正確性が向上する。また、第3の形態として、第1又は第2の形態の位置算出方法において、前記位置算出演算は、少なくとも前記演算結果の変動を確率変数とした正規分布モデルに基づいて定められた演算であり、過去の前記位置算出演算の結果に基づいて前記正規分布モデルを更新することを更に含む、位置算出方法を構成することとしてもよい。この第3の形態によれば、少なくとも演算結果の変動を確率変数とした正規分布モデルに基づいて定められた位置算出演算を行う。つまり、演算結果の変動の分布が正規分布に従うと仮定する。その上で、正規分布モデルを固定的に設定するのではなく、過去の位置算出演算の結果に基づいて正規分布モデルを更新することで、時間変化に対応した適切な正規分布モデルに基づく位置算出が可能となる。また、第4の形態として、第3の形態の位置算出方法において、前記受信信号に基づいて擬似距離を測定することを更に含み、前記位置算出演算は、前記演算結果の変動を確率変数とした正規分布モデルと、前記擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた演算である、位置算出方法を構成することとしてもよい。この第4の形態によれば、測位用信号の受信信号に基づいて擬似距離を測定する。そして、演算結果の変動を確率変数とした正規分布モデルと、擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた位置算出演算を行う。演算結果の変動を確率変数とした正規分布モデルに加えて、擬似距離に含まれる誤差を確率変数とした確率分布モデルを用いることで、位置算出演算を適切に行うことができる。より具体的には、第5の形態として、第4の形態の位置算出方法における前記擬似距離に含まれる誤差を確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとする、位置算出方法を構成することができる。通常の測位環境においては、擬似距離に含まれる誤差の分布は正規分布に従うと仮定することができる。しかし、例えばマルチパス環境のように、擬似距離に含まれる誤差の分布を正規分布と仮定することが必ずしも適切でない場合もある。そこで、第5の形態のように、擬似距離に含まれる誤差を確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとして位置算出演算を行うことで、測位環境に見合った適切な位置算出を実現することが可能となる。また、第6の形態として、第3~第5の何れかの形態の位置算出方法において、前記位置算出演算は、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する演算であり、前記正規分布モデルは、当該演算結果の各要素の変動を確率変数としたモデルである、位置算出方法を構成することとしてもよい。この第6の形態によれば、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する位置算出演算を行う。その際、演算結果の各要素の変動を確率変数とした正規分布モデルを用いることで、各要素に時間的な相関を持たせて位置算出を行うことが可能となる。また、第7の形態として、第1~第6の何れかの形態の位置算出方法において、前記位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、前記受信信号を用いた所定の収束演算を行って算出することを更に含み、前記位置算出演算は、過去に前記収束演算によって算出された前記所定のパラメーターの値を用いる演算であることを含む、位置算出方法を構成することとしてもよい。この第7の形態によれば、位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、受信信号を用いた所定の収束演算を行って算出する。そして、過去に収束演算によって算出された所定のパラメーターの値を用いる位置算出演算を行う。これにより、特定の算出パラメーターについては、時間的な相関を持たせずに、所定の収束演算を行って値を算出することができる。そして、その算出結果を、所与の確率分布モデルに基づいて定められた位置算出演算に反映させることで、位置算出を効果的に行うことが可能となる。以下、図面を参照して、本発明の好適な実施形態の一例について説明する。本実施形態は、衛星測位システムの一種であるGPSGlobalPositioningSystemを適用した実施形態である。なお、本発明を適用可能な実施形態が以下説明する実施形態に限定されるわけでないことは勿論である。1.原理本実施形態における位置算出方法位置算出の原理について説明する。GPSを利用した衛星測位システムにおいて、測位用衛星の一種であるGPS衛星は、エフェメリスやアルマナックといった衛星軌道データを含む航法データを、測位用の衛星信号の一種であるGPS衛星信号に乗せて発信している。GPS衛星信号は、拡散符号の一種であるC/ACoarseandAcquisitionコードによって、スペクトラム拡散方式として知られるCDMACodeDivisionMultipleAccess方式によって変調された1.57542[GHz]の通信信号である。C/Aコードは、コード長1023チップを1PNフレームとする繰返し周期1msの擬似ランダム雑音符号であり、各GPS衛星に固有のコードである。GPSでは、複数のGPS衛星の位置や各GPS衛星からGPS受信機位置算出装置までの擬似距離といった、GPS衛星信号の受信信号に係る測定量を用いて、GPS受信機の位置位置座標や時計誤差クロックバイアスを演算する位置算出演算を行う。擬似距離は、衛星サーチを行ってGPS衛星信号を捕捉し、その結果として得られるメジャメント情報を用いて測定される。衛星サーチには、GPS衛星信号を受信した信号に対する位相方向の相関演算以下、位相サーチと称す。と、周波数方向の相関演算以下、周波数サーチと称す。とが含まれる。位相サーチを行うことで、GPS受信機がGPS衛星信号を受信した際のC/Aコードの位相以下、コード位相と称す。をメジャメント情報として取得できる。また、周波数サーチを行うことで、GPS受信機が受信したGPS衛星信号の受信周波数やドップラー周波数をメジャメント情報として取得できる。観念的には、GPS衛星とGPS受信機との間には複数のC/Aコードが並んでいると考えることができる。GPS衛星とGPS受信機との間の距離は、ちょうどC/Aコードの整数倍の長さになるとは限らず、端数部分が生じ得る。この擬似距離の端数部分に相当するのがコード位相である。擬似距離の整数部分は、GPS受信機及びGPS衛星の概略位置から求められる。本実施形態では、GPS受信機位置算出装置の位置ベクトルP及びクロックバイアス時計誤差bを未知数とし、これらを成分とするGPS受信機の状態ベクトルXを次式1のように定義する。式1において、x,y,zは、測位を行う座標系におけるGPS受信機位置算出装置の位置ベクトルPの各軸の位置成分である。後々の便宜のため、位置成分x,y,z及びクロックバイアスbを、状態Xの成分という意味でXi=x,y,z,bT=X1,X2,X3,X4Tと表記する。但し、下付きの添え字のi=1,2,3,4は、状態Xの成分の番号を示す。上付きの添え字のTは行列の転置を示す。また、本実施形態の数式では、ベクトルを意味する矢印の表記を省略して簡潔に記載する。GPS受信機の真の位置ベクトルをP=P0δP、真のクロックバイアスをb=b0δbとおく。但し、δP及びδbは、それぞれ位置ベクトル及びクロックバイアスの初期値に適用される未知の変化量修正量であり、それぞれ位置変化量ベクトル及びクロックバイアス変化量と呼称する。また、位置変化量ベクトルδP及びクロックバイアス変化量δbを纏めて、次式2のように状態変化量ベクトルδXを定義する。従来から用いられているGPSの測位計算では、例えば次式3の観測方程式を利用して、状態変化量ベクトルδXを求める。本明細書では、GPS衛星のうちGPS受信機が捕捉に成功した衛星以下、捕捉衛星と称す。の番号をkと表記し、原則的に各諸量について下付きの添え字で表記する。また、本実施形態では、GPS受信機が合計K個のGPS衛星の捕捉に成功したものとして説明する。つまり、k=1,2,,Kとして説明する。式3の左辺のδρは擬似距離修正量ベクトルであり、次式4のように表される。つまり、擬似距離修正量ベクトルδρは、各捕捉衛星の擬似距離修正量δρkを成分とするベクトルである。また、式3の右辺の行列Gは、捕捉衛星の衛星配置を決定付ける幾何行列であり、次式5ように表される。但し、lkは、GPS受信機からk番目の捕捉衛星に向かう視線方向のベクトル以下、視線ベクトルと称す。である。また、式3の右辺のEは擬似距離誤差ベクトルであり、次式6のように表される。但し、εkは、k番目の捕捉衛星からGPS受信機までの擬似距離に含まれる誤差以下、擬似距離誤差と称す。である。擬似距離誤差ベクトルEは、各捕捉衛星の擬似距離誤差εkを成分とするベクトルである。式3における未知数は、位置変化量ベクトルδP=δx,δy,δz及びクロックバイアス変化量δbの合計4個である。従って、式3は、K≧4の場合に解くことができる。K>4の場合は、式3は過決定方程式となる。この場合は、K個の捕捉衛星の擬似距離誤差εの二乗和を最小とする解を求める最小二乗法が一般的に用いられる。具体的には、式3の最小二乗解は、例えば次式7のように求められる。式3の観測方程式を利用した位置算出演算では、演算時刻毎に、GPS衛星信号を受信して擬似距離を測定し、式7の最小二乗解を解析的に求める。この場合、状態変化量ベクトルδXを、演算時刻毎に独立として計算している。しかし、GPS受信機位置算出装置は、通常は連続的に移動する。例えば、GPS受信機を携行したユーザーが通常歩行や通常走行しているのであれば、演算時刻毎のGPS受信機の移動距離にそれほど大きな差は現れないはずである。つまり、ユーザーが突飛な動作をしない限り、GPS受信機の位置変化量ベクトルδPの変動はほぼ一定となることが期待される。GPS受信機を自動車等の移動体に設置した場合も同様である。このことから、本願発明者は、位置算出演算の演算結果の変動を確率変数とした所与の確率分布モデルを定義することを考えた。例えば、ユーザーが通常歩行しており、1秒間隔で位置算出演算を行うとすると、1秒毎のGPS受信機の移動距離は数メートル程度となり、個人差を考慮しても、その誤差幅は数十センチ程度になると考えられる。従って、歩行を想定した場合、位置算出演算の演算結果のうち、位置変化量ベクトルδPについては、期待値が数メートル程度、標準偏差が数十センチ程度と仮定できるであろう。また、GPS受信機を自動車に設置する場合を考えると、位置変化量ベクトルδPは、期待値が数十メートル程度、標準偏差が数メートル程度と仮定できるであろう。クロックバイアスbについても同様に考えると、短い時間間隔で考えた場合、クロックバイアスbがそれほど大きく変動するとは考えにくい。そのため、クロックバイアス変化量δbについても確率分布モデルを仮定できるであろう。そこで、本願発明者は、位置ベクトルP及びクロックバイアスbを成分とする状態ベクトルXについて、その変動である状態変化量ベクトルδXを確率変数とする正規分布型の確率分布モデル以下、状態変動正規分布モデルと称す。を仮定した。図1は、状態変動正規分布モデル状態変動GMGaussianModelの説明図である。図1には、状態変化量ベクトルδXの1つの成分δXiに着目したモデルを図示している。横軸は状態変化量の成分δXiであり、縦軸はその確率密度pδXiである。状態変化量δXiの期待値μを状態変化量期待値と呼称し、状態変化量δXiの標準偏差σを状態変化量標準偏差と呼称する。状態変化量δXiを確率変数とし、その分布が正規分布Nμ,σ2に従うと仮定する。しかし、本実施形態では、状態変化量ベクトルδXには、位置変化量ベクトルPの各成分δx,δy,δz及びクロックバイアス変化量δbが含まれる。そこで、一次元の正規分布モデルを拡張して、次式8で与えられる多次元状態変動正規分布モデルpδXδXを用いることとする。式8において、上付きの括弧書きのmは確率変数として取り扱う算出パラメーターの数を示す。算出パラメーターとは、位置算出演算によって算出可能なパラメーターであり、本実施形態ではδx,δy,δz,δbの4つのパラメーターがこれに相当する。4つの算出パラメーター全てを確率変数とするのであればm=4となる。詳細は後述するが、4つの算出パラメーターのうちの一部の算出パラメーターを確率変数とすることも可能である。この場合は、確率変数とする算出パラメーターの数がmとなる。式8の多次元状態変動正規分布モデルpδXδXは、非独立のm個の確率変数の同時発生確率を与える確率密度関数である。μδXは、状態変化量ベクトルδXの各成分の期待値でなる状態変化量期待値ベクトルである。また、SδXは、状態変化量ベクトルδXの成分間の共分散でなる状態変化量共分散行列である。これらについては、数式を用いて詳細に後述する。本実施形態では、上記の多次元状態変動正規分布モデルpδXδXを適用して、最尤推定法を利用して状態変化量ベクトルδXの最尤推定解を求める。これは、測位用信号を受信し、測位用信号の受信信号を用いた演算であって、位置算出演算の演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことに相当する。状態変化量ベクトルδXの最尤推定解が求まれば、当該最尤推定解を状態ベクトルの初期値X0に加算することで、真の状態ベクトルXが推定できる。さて、本実施形態の位置算出演算では、GPS衛星信号の受信信号に基づいて擬似距離を測定する。そして、上記の状態変動正規分布モデルと、擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた位置算出演算を行う。つまり、擬似距離誤差εも確率変数として取り扱い、所与の確率分布モデルを仮定して演算を行う。本実施形態では、擬似距離誤差εを確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとする。図2は、擬似距離誤差正規分布モデル擬似距離誤差GMの説明図である。図2において、横軸は擬似距離誤差εであり、縦軸はその確率密度pεεである。擬似距離誤差εの期待値μεを擬似距離誤差期待値と呼称し、擬似距離誤差εの標準偏差σεを擬似距離誤差標準偏差と呼称する。擬似距離誤差εが、正規分布Nμε,σε2に従うと仮定して観測方程式の解を求める手法が最小二乗法である。広く用いられるのはμε=0,σε=1とする擬似距離誤差正規分布モデルであり、この場合の最小二乗解が式7と一致する。上記の擬似距離誤差正規分布モデルpεεは、GPS受信機の天空が開けた環境例えばオープンスカイ環境において有用である。オープンスカイ環境では、GPS衛星から品質の良いGPS衛星信号を受信することができることから、擬似距離誤差εはゼロ近傍の値となることが期待されるためである。しかし、GPS受信機の天空が開けていない環境では、擬似距離誤差εの挙動が変化する。その典型例はマルチパス環境である。図3は、マルチパス環境の説明図である。直接波信号を点線で示し、間接波信号を一点鎖線で示す。間接波信号の存在により、GPS受信機が測定する擬似距離には誤差が生ずる。擬似距離は、GPS受信機が位相サーチを行って取得したコード位相を用いて算出されるが、マルチパス環境では、このコード位相に誤差が重畳する。コード位相の誤差は、真のコード位相に対して、正の誤差となる場合及び負の誤差となる場合がある。図4は、コード位相検出の原理の説明図である。図4では、横軸をコード位相、縦軸を相関値として、C/Aコードの自己相関値の概略例を示している。なお、以下の説明では、相関値というときは、相関値の大きさ絶対値を意味するものとする。C/Aコードの自己相関値は、理想的にはピーク値を頂点とする左右対称の略三角形の形状で表される。この場合に、相関値のピーク値以下、相関ピーク値と称す。に対応する位相が、受信したGPS衛星信号のC/Aコードの位相である。相関ピーク値を検出するために、例えば、あるコード位相に対して、一定量だけ進んだコード位相以下、進み位相と称す。と、一定量だけ遅れたコード位相以下、遅れ位相と称す。とを考える。そして、進み位相の相関値と遅れ位相の相関値とが等しくなる相関値を相関ピーク値として検出する。図4では、相関値は左右対称の略三角形の形状であるため、進み位相の相関値と遅れ位相の相関値との中心のコード位相の相関値が相関ピーク値となり、対応するコード位相がピーク位相として検出される。この検出されたピーク位相のことを検出ピーク位相と称する。図4の相関値の形状は理想形状であるが、マルチパス環境では相関値の形状が変化する。図5及び図6は、マルチパス環境における相関値の形状の一例を示す図である。図5は、間接波信号が直接波信号と同位相で到達した場合の相関値のグラフの一例であり、図6は、間接波信号が直接波信号と逆位相で到達した場合の相関値のグラフの一例である。これらの図には、直接波信号と、間接波信号と、この直接波信号と間接波信号とを合成したマルチパス信号とのそれぞれに対応する相関値のグラフを示している。横軸はコード位相、縦軸は相関値である。間接波信号に対する相関値は、直接波信号に対する相関値と同様に略三角形の形状をなしているが、間接波信号の相関ピーク値の大きさは、直接波信号の相関ピーク値よりも小さい。これは、GPS衛星から送出されたGPS衛星信号が、建物や地面に反射したり障害物を透過することによって、送出された時点における信号強度が、受信時には弱められていることによるものである。また、間接波信号のピーク位相は、直接波信号のピーク位相よりも遅れている。これは、GPS衛星から送出されたGPS衛星信号が、建物や地面に反射したり障害物を回折することによって、GPS衛星からGPS受信機までの伝搬距離が長くなったことによるものである。マルチパス信号に対する相関値は、直接波信号の相関値と間接波信号の相関値との和となるため、三角形状が歪んでピーク値を中心とした左右対称とはならない。間接波信号が直接波信号と同位相でGPS受信機に到達した場合は、直接波信号と間接波信号とは互いに強め合う。このため、合成波信号の相関値は、直接波信号に対する相関値の大きさと間接波信号に対する相関値の大きさとの合算値となる。この場合、相関値の形状は、例えば図5に示すような形状となる。図5において、検出ピーク位相は、真のピーク位相よりも遅れた位相となる。それに対し、間接波信号が直接波信号よりも、例えば半周期以上1周期未満の範囲等の逆位相となった場合は、直接波信号と間接波信号とは互いに弱め合う。このため、合成波信号の相関値は、直接波信号に対する相関値の大きさから間接波信号に対する相関値の大きさを減じた値となる。この場合、相関値の形状は、例えば図6に示すような形状となる。図6において、検出ピーク位相は、真のピーク位相よりも進んだ位相となる。ここで、検出ピーク位相と真のピーク位相との位相差をコード位相誤差と定義する。そして、便宜的に、検出ピーク位相が真のピーク位相よりも遅れている場合のコード位相誤差の符号を正、検出ピーク位相が真のピーク位相よりも進んでいる場合のコード位相誤差の符号を負と定義する。この場合、図5ではコード位相誤差は正となり、図6ではコード位相誤差は負となる。前述したように、擬似距離はコード位相を用いて算出されるため、コード位相の誤差は擬似距離に誤差として重畳される。すなわち、図5のようにコード位相に正の誤差が含まれると、擬似距離誤差εは正の誤差となる。一方、図6のようにコード位相に負の誤差が含まれると、擬似距離誤差εは負の誤差となる。このように、マルチパス環境等においては、コード位相誤差の正負の符号に応じて、擬似距離誤差εは正の値にもなり得るし、負の値にもなり得る。勿論、擬似距離誤差εにはコード位相誤差以外の誤差成分も含まれる。そのため、これらの誤差成分の大きさ如何によっても、擬似距離誤差εは種々の値をとり得る。このことから、本願発明者は、例えばマルチパス環境においては、擬似距離誤差εの分布を正規分布と仮定するのではなく、複数の正規分布を混合した正規混合分布と仮定することが適切であると判断した。擬似距離誤差εに正規混合分布を仮定した確率分布モデルとして、本実施形態では2種類のモデルを例示する。図7は、擬似距離誤差GMMの一例である第1の擬似距離誤差正規混合分布モデル第1の擬似距離誤差GMMGaussianMixtureModelの説明図である。第1の擬似距離誤差GMMfεは、第1の正規分布モデルp1εと、第2の正規分布モデルp2εとの2種類の正規分布モデルを合成した関数として定義される。図7において、横軸は擬似距離誤差ε単位はメートルであり、縦軸はその確率密度fεである。また、第1の正規分布モデルp1εを点線で示し、第2の正規分布モデルp2εを一点鎖線で示し、第1の擬似距離誤差GMMfεを太実線でそれぞれ示す。第1の正規分布モデルp1εは、GPS受信機が直接波信号を受信することを想定した擬似距離誤差εの確率分布モデルである。間接波信号の影響がなければ、擬似距離誤差εは0メートルの近傍の値となることが期待される。そのため、擬似距離誤差期待値を0メートルμε=0とし、擬似距離誤差標準偏差を比較的小さな値である10メートルσε=10とする正規分布関数N0,102によって、第1の正規分布モデルp1εをモデル化した。第2の正規分布モデルp2εは、GPS受信機が直接波信号及び間接波信号を受信することを想定した擬似距離誤差εのモデル関数である。前述したように、コード位相誤差の正負の符号により、擬似距離誤差εは正負の何れの値もとり得る。そのため、第2の正規分布モデルp2εの擬似距離誤差期待値を、正負の中心である0メートルμε=0とした。また、コード位相誤差の大きさ如何によっては、擬似距離誤差εに非常に大きな誤差が含まれ得る。そのため、擬似距離誤差標準偏差を大きめに見積もって100メートルσε=100とした。すなわち、第2の正規分布モデルp2εを正規分布関数N0,1002によってモデル化した。第1の正規分布モデルp1εと第2の正規分布モデルp2εとを混合する割合合成比率は、次のように決定した。GPS受信機がマルチパスの影響を受ける割合は、一般的にはそれほど高くないと考えられる。そこで、GPS受信機が1割の確率でマルチパスの影響を受けると仮定して、第1の正規分布モデルp1εの重みを9/109割とし、第2の正規分布モデルp1εの重みを1/101割とした。従って、第1の正規分布モデルp1εを9/10倍した関数と、第2の正規分布モデルp2εを1/10倍した関数とを混合合成した関数として、第1の擬似距離誤差GMMfεをモデル化した。図8は、擬似距離誤差GMMの別例である第2の擬似距離誤差正規混合分布モデル第2の擬似距離誤差GMMの説明図である。第2の擬似距離誤差GMMfεは、図7の第1の擬似距離誤差GMMと同様に、第1の正規分布モデルp1εと第2の正規分布モデルp2εとの2種類の正規分布モデルを混合合成した確率密度関数である。第1の正規分布モデルp1εは、GPS受信機が直接波信号を受信することを想定した擬似距離誤差εのモデル関数として、正規分布関数N0,102でモデル化した。また、第2の正規分布モデルp2εは、マルチパス環境においてGPS受信機が間接波信号を受信することを想定した擬似距離誤差εのモデル関数として、正規分布関数N100,202でモデル化した。図7の第1の擬似距離誤差GMMfεと異なるのは、第2の正規分布モデルp2εを、直接波信号及び間接波信号を受信することを想定したモデル関数として定義するのではなく、間接波信号単体を受信することを想定したモデル関数として定義した点にある。これは、例えば、GPS受信機の近傍に高層ビル等の障害物が存在する環境において起こり得る。GPS受信機の天空が完全に開けていないために、GPS衛星から送出された信号は障害物に反射し、間接波信号としてGPS受信機に到達する。しかし、障害物に遮られて、直接波信号がGPS受信機に直接到達しない環境である。アーバンキャニオン環境が典型例である。本願発明者が実験を行った結果、アーバンキャニオン環境では、擬似距離誤差εとして100~150メートル程度の正の誤差が含まれる傾向があることがわかった。そこで、第2の擬似距離誤差GMMfεでは、第2の正規分布モデル関数p2εの擬似距離誤差期待値を100メートルμε=100とし、擬似距離誤差標準偏差を20メートルσε=20とした。すなわち、p2ε=N100,202としてモデル化した。第1及び第2の正規分布モデルの合成比率は、例えばGPS受信機が2割の確率でマルチパスの影響を受けると仮定して決定した。すなわち、第1の正規分布モデルp1εの重みを8/108割とし、第2の正規分布モデルp2εの重みを2/102割とした。従って、第2の正規分布モデルp1εを8/10倍した関数と、第2の正規分布モデルp2εを2/10倍した関数とを混合合成した関数として、第2の擬似距離誤差GMMfεをモデル化した。上記のように、擬似距離誤差εを確率変数とする確率分布モデルとして、正規分布モデル及び正規混合分布モデルの2種類のモデルを定義した。位置算出を行う際には、測位環境に応じて、擬似距離誤差正規分布モデルと、擬似距離誤差正規混合分布モデルとを使い分けて位置算出演算を行う。以下、それぞれの確率分布モデルを適用する場合の具体的な位置算出方法について、数式を用いて詳細に説明する。1-1.擬似距離誤差正規分布モデル擬似距離誤差GMを適用する場合1第1の位置算出方法第1の位置算出方法は、擬似距離誤差正規分布モデルを適用し、状態変化量ベクトルδXの全ての算出パラメーターを確率変数として位置算出する方法である。この方法では、状態変化量ベクトルδXの4つの算出パラメーターδx,δy,δz,δb全てを確率変数とする。この場合、式8の多次元状態変動正規分布モデルpδXδXは、4次元状態変動正規分布モデルpδXδXとなる。最尤推定法は、尤度と呼ばれる指標値に基づいて、所与の観測データを利用して推定対象とするパラメーターの最尤解を推定する統計学的手法である。尤度は、ある観測データについて、適用する確率モデルがどの程度当てはまっているかを確率で表す尺度である。本実施形態において、尤度Lは、式6の擬似距離誤差ベクトルEのK衛星の擬似距離誤差εkk=1,2,,Kの同時発生確率と、状態変化量ベクトルδXの4つの算出パラメーターδXi=δx,δy,δz,δbi=1,2,3,4の同時発生確率との積で表される。つまり、尤度Lは、次式9のように表される。式9において、pεkεはk番目の捕捉衛星の擬似距離誤差正規分布モデルを示す。なお、図2で定義した擬似距離誤差正規分布モデルpεεは、全ての捕捉衛星についてパラメーター値を共通化した共通のモデルとしてもよいし、捕捉衛星毎にパラメーター値を個別化した異なるモデルとしてもよい。また、式9において、pδXiδXiはi番目の算出パラメーターの状態変動正規分布モデルを示す。1行目から2行目に変形する際に、4つの算出パラメーターの同時発生確率を与える4次元状態変動正規分布モデルpδXδXに書き換えた。式9の尤度Lを最大化する状態変化量ベクトルδXを最尤推定解として求める。尤度Lをそのまま用いて計算するのは困難であるため、本実施形態では、その対数をとった対数尤度logLを用いて計算する。対数尤度logLは、次式10のように計算される。式10において、constは定数項である。また、μEは、K衛星の擬似距離誤差εの期待値μεk=με1,με2,,μεKを成分とする擬似距離誤差期待値ベクトルであり、次式11で与えられる。WEは、擬似距離誤差逆共分散行列であり、K衛星分の擬似距離誤差εkの共分散でなる擬似距離誤差共分散行列の逆行列として計算される。具体的には、次式12で与えられる。μδXは、式8で登場した状態変化量期待値ベクトルであり、次式13で与えられる。但し、μδx,μδy,μδz,μδbは、それぞれδx,δy,δz,δbの期待値である。また、SδXは、式8で登場した状態変化量共分散行列であり、次式14で与えられる。但し、4×4の各成分は、δx,δy,δz,δbの4成分間の共分散である。また、WδXは、状態変化量逆共分散行列であり、状態変化量共分散行列SδXの逆行列として計算される。具体的には、次式15で与えられる。式10を最大化するδXを求めるためには、dlogL/dδX=0を解けばよい。途中の計算は省略するが、最終的に、最尤推定解δXは次式16のように求められる。但し、δρは式4の擬似距離変化量ベクトルであり、Gは式5の幾何行列である。各演算時刻それぞれについて、式16に従って状態変化量ベクトルδXの最尤推定解を求める。そして、求めた状態変化量ベクトルδXの最尤推定解を用いて、状態ベクトルXを算出する。また、本実施形態の位置算出演算では、上記の最尤推定解の算出と併せて、過去の位置算出演算の結果に基づいて状態変動正規分布モデルを更新する処理を行う。ここで、演算時刻tにおける状態変化量ベクトルδXをδXtと表記する。また、演算時刻tにおける状態変化量期待値ベクトルμδX及び状態変化量共分散行列SδXを、それぞれμδXt及びSδXtと表記する。この場合、状態変化量期待値ベクトルμδXt及び状態変化量共分散行列SδXtを、それぞれ次式17及び18に従って更新する。式17は、1つ前の演算時刻に計算された状態変化量ベクトルδXtで、今回の演算時刻における状態変化量期待値ベクトルμδXt1を更新することを意味する。これは、過去の位置算出演算の結果に基づいて、演算結果状態変化量ベクトルδXの変動の基準状態変化量期待値ベクトルμδXtを更新することに相当する。また、式18は、1つ前の演算時刻に計算された状態変化量共分散行列SδXt等の諸量を用いて、今回の演算時刻における状態変化量共分散行列SδXt1を更新することを意味する。状態変化量期待値ベクトルμδXは、状態変動正規分布モデルの期待値に係るベクトルである。また、状態変化量共分散行列SδXは、状態変動正規分布モデルの標準偏差共分散に係る行列である。よって、式17及び18は、過去の位置算出演算の結果に基づいて状態変動正規分布モデルを更新する式であると言える。2第2の位置算出方法第2の位置算出方法は、擬似距離誤差εに正規分布モデルを適用し、状態変化量ベクトルδXの一部の算出パラメーターを確率変数として位置算出する方法である。この方法では、状態変化量ベクトルδXを構成する4つの算出パラメーターδx,δy,δz,δbのうちの一部の算出パラメーターを確率変数とし、残余の算出パラメーターを非確率変数として取り扱う。確率変数及び非確率変数とする算出パラメーターは、適用するシステムに応じて適宜選択可能である。ここでは、ユーザーの位置変化量ベクトルδPを構成するδx,δy,δzの3つの算出パラメーターを確率変数とし、クロックバイアス変化量δbを非確率変数とする場合を例に挙げて説明する。この場合、式8の多次元状態変動正規分布モデルpδXδXは、3次元状態変動正規分布モデルpδXδXとなる。この場合、クロックバイアス変化量δbについては正規分布モデルを仮定しないため、式10に相当する対数尤度logLは、次式19のように変更される。また、式13で定義した状態変化量期待値ベクトルμδXは、次式20のように変更される。但し、クロックバイアス変化量δbを非確率変数とするため、それに対応する4番目の成分をゼロとしている。また、式14で定義した状態変化量共分散行列SδXは、次式21のように変更される。但し、クロックバイアス変化量δbを非確率変数とするため、4×4の行列の4行目及び4列目の成分を全てゼロとしている。この場合も、最尤推定解δXを式16のように求めることができる。また、式20で定義した状態変化量期待値ベクトルμδXと、式21で定義した状態変化量共分散行列SδXとを、それぞれ式17及び式18に従って更新する。この場合、クロックバイアス変化量δbに対応する成分はゼロとして更新しない。また、非確率変数とするクロックバイアス変化量δbについては、受信信号を用いた所定の収束演算を行って解を求める。具体的には、式3の観測方程式に対して擬似距離誤差正規分布モデルを適用して、例えばニュートン法を利用した収束演算を行って、クロックバイアス変化量δbの近似解を求めることになる。収束演算によって算出されたクロックバイアス変化量δbは、次回の演算時刻での位置算出演算において擬似距離を近似的に求めるために使用する。1-2.擬似距離誤差正規混合分布モデル擬似距離誤差GMMを適用する場合3第3の位置算出方法第3の位置算出方法は、擬似距離誤差εに正規混合分布モデルを適用し、状態変化量ベクトルδXの全ての算出パラメーターを確率変数として位置算出する方法である。考え方は第1の位置算出方法と同じであるが、擬似距離誤差の確率分布モデルとして、正規分布モデルではなく、正規混合分布モデルを適用する点が異なる。この場合、擬似距離誤差GMMfεと、4次元状態変動正規分布モデルpδXδXとを用いて、対数尤度logLを次式22のように定式化する。式22において、fkεkは、k番目の捕捉衛星の擬似距離誤差GMMを示し、次式23のように一般化される。式23において、lは、混合合成する正規分布モデルの番号であり、L個l=1,2,,Lの正規分布モデルを混合する場合の式として記載している。αlは、l番目の正規分布モデルNlの重み合成割合である。μkは、k番目の捕捉衛星に係るl番目の正規分布モデルの期待値であり、σkは、k番目の捕捉衛星に係るl番目の正規分布モデルの標準偏差である。なお、上記の擬似距離誤差GMMはfεは、全ての捕捉衛星についてパラメーター値を共通化した共通のモデルとしてもよいし、捕捉衛星毎にパラメーター値を個別化した異なるモデルとしてもよい。式22の対数尤度logLを最大化するδXを求めるためには、dlogL/dδX=0を解けばよい。しかし、正規混合分布GMMでは、対数の加算など複雑な項が存在するため、解析的に解を求めることが困難である。そこで、本実施形態では、最適化問題に対するアルゴリズムの一種であるEMExpectationMaximumアルゴリズムを用いて最尤推定解を求める。EMアルゴリズムは、確率分布モデルのパラメーターを最尤推定法に基づいて推定する手法の一種である。観測可能なデータサンプルの他に、観測不可能なデータ隠しパラメーターが存在する場合であって、当該隠しパラメーターが確率分布モデルに依存する場合に用いられる手法である。EMアルゴリズムは反復法の一種であり、尤度関数の条件付き確率に関する期待値を計算するEExpectationステップと、期待値を最大化する解を求めるMMaximumステップとの2つのステップを反復して最尤推定解を求める。途中の計算は省略するが、EMアルゴリズムを適用すると、状態変化量ベクトルδXの最尤推定解は、次式24のように求まる。但し、δX0は各演算時刻における状態変化量ベクトルδXの初期値である。式24におけるMlδX0は、次式25で与えられる。また、式25におけるhklδX0は、次式26で与えられる。また、式24におけるμlは、次式27で与えられる。各演算時刻において、式24~式27に従って最尤推定解δXを求める。また、擬似距離誤差正規分布モデルを適用した場合と同様に、式13で定義した状態変化量期待値ベクトルμδXと、式14で定義した状態変化量共分散行列SδXとを、演算時刻毎に更新する。具体的には、状態変化量期待値ベクトルμδXt及び状態変化量共分散行列SδXtを、それぞれ次式28及び29に従って更新する。式29におけるHδX及びFδXは、それぞれ次式30及び31で与えられる。また、式31におけるRlは、次式32で与えられる。4第4の位置算出方法第4の位置算出方法は、擬似距離誤差εに正規混合分布モデルを適用し、状態変化量ベクトルδXの一部の算出パラメーターを確率変数として位置算出する方法である。この方法では、第2の位置算出方法と同様に、状態変化量ベクトルδXを構成する4つの算出パラメーターδx,δy,δz,δbのうちの一部の算出パラメーターを確率変数とし、残余の算出パラメーターを非確率変数として取り扱う。擬似距離誤差の確率分布モデルとして、正規分布モデルではなく、正規混合分布モデルを適用する点が第2の位置算出方法と異なる。この場合も、第3の位置算出方法と同様の数式を用いて、状態変化量ベクトルδXの最尤推定解を求めることができる。第3の位置算出方法と異なるのは、確率変数として取り扱う算出パラメーターについてのみ正規分布モデルを適用して計算する点である。つまり、第2の位置算出方法と同様に、状態変化量ベクトルδXの成分のうち、非確率変数として取り扱う算出パラメーターに対応する成分をゼロとして計算する。例えば、クロックバイアス変化量δbを非確率変数とする場合は、クロックバイアス変化量δbに対応する成分をゼロとして計算する。そして、クロックバイアス変化量δbについては、受信信号を用いた所定の収束演算を行って解を算出する。具体的には、式3の観測方程式に対して擬似距離誤差正規混合分布モデル擬似距離誤差GMMを適用して、例えばEMアルゴリズムを利用した収束演算を行って、クロックバイアス変化量δbの近似解を求めることになる。この場合、収束演算によって算出されたクロックバイアス変化量δbは、次回の演算時刻での位置算出演算において擬似距離を近似的に求めるために使用する。2.実験結果本実施形態の位置算出方法で位置算出演算を行った実験結果について説明する。従来の手法と本実施形態の手法とのそれぞれを用いて、位置算出演算を行って実際に位置を算出する実験を行った。ユーザーに位置算出装置を携行させ、予め定められた周回経路に沿って歩行させながら位置を算出した。図9は、位置算出結果の一例を示す実験結果である。グラフの横方向は東西方向を示し、縦方向は南北方向を示す。単位はそれぞれメートル[m]である。ここで示す結果は、状態変動の確率分布モデルとして正規分布を適用し、擬似距離誤差の確率分布モデルとして正規混合分布GMMを適用した結果である。従来の位置算出方法で得られた軌跡を白丸のプロットで示し、本実施形態の位置算出方法で得られた軌跡を黒丸のプロットで示す。この結果を見ると、従来の位置算出方法を適用した場合は、所々に位置飛びが発生しており、それほど正確な位置は得られていないことがわかる。しかし、本実施形態の位置算出方法を適用した場合は、周回経路に沿った滑らかな軌跡が得られていることがわかる。位置算出の演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行ったことで、算出位置の連続性が飛躍的に向上している。図10は、上記の歩行試験において位置演算精度を計測したグラフである。グラフの横方向は東西方向の位置誤差を示し、縦方向は南北方向の位置誤差を示す。単位はそれぞれメートル[m]である。グラフの中心が位置誤差0[メートル]に相当し、プロットがグラフの中心に近いほど位置算出精度が高いことを意味する。このグラフを見ると、従来の位置算出方法を適用した場合は、プロットが全体的に広くばらついており、位置算出精度が低くなっていることがわかる。それに対し、本実施形態の位置算出方法を適用した場合は、プロットが全体的に中心部に集中しており、従来の手法と比べて位置算出精度が格段に向上していることがわかる。図11は、位置算出精度を示す他の実験結果の一例である。ここで示す結果は、状態変動の確率分布モデル及び擬似距離誤差の確率分布モデルとして、共に正規分布モデルを適用した結果である。図の見方は図10と同じである。この結果を見ると、擬似距離誤差の確率分布モデルとして正規分布モデルを適用した場合も、位置誤差を示すプロットがグラフの中心部に集中しており、位置算出精度が格段に向上していることがわかる。上記の実験結果により、本実施形態の位置算出方法を用いることで、位置算出の正確性が向上することが実証された。3.実施例次に、上記の位置算出方法を用いて位置算出を行う位置算出装置の実施例について説明する。ここでは、位置算出装置を具備した電子機器の一例である携帯型電話機の実施例として説明する。なお、本発明を適用可能な実施例が以下説明する実施例に限定されるわけでないことは勿論である。3-1.携帯型電話機の構成図12は、携帯型電話機1の機能構成の一例を示すブロック図である。携帯型電話機1は、GPSアンテナ5と、GPS受信部10と、ホスト処理部30と、操作部40と、表示部50と、携帯電話用アンテナ60と、携帯電話用無線通信回路部70と、記憶部80と、時計部90とを備えて構成される。GPSアンテナ5は、GPS衛星から発信されているGPS衛星信号を含むRFRadioFrequency信号を受信するアンテナであり、受信信号をGPS受信部10に出力する。GPS受信部10は、GPSアンテナ5から出力された信号に基づいて携帯型電話機1の位置を算出する回路或いは装置であり、GPS受信機に相当する機能ブロックである。本実施例では、このGPS受信部10が位置算出装置に相当する。GPS受信部10は、RF受信回路部11と、ベースバンド処理回路部20とを備えて構成される。GPS受信部10は、GPS衛星信号を受信する受信部に相当する。なお、RF受信回路部11と、ベースバンド処理回路部20とは、それぞれ別のLSILargeScaleIntegrationとして製造することも、1チップとして製造することも可能である。RF受信回路部11は、RF信号の受信回路である。回路構成としては、例えば、GPSアンテナ5から出力されたRF信号をA/D変換器でデジタル信号に変換し、デジタル信号を処理する受信回路を構成してもよい。また、GPSアンテナ5から出力されたRF信号をアナログ信号のまま信号処理し、最終的にA/D変換することでデジタル信号をベースバンド処理回路部20に出力する構成としてもよい。後者の場合には、例えば、次のようにRF受信回路部11を構成することができる。すなわち、所定の発振信号を分周或いは逓倍することで、RF信号乗算用の発振信号を生成する。そして、生成した発振信号を、GPSアンテナ5から出力されたRF信号に乗算することで、RF信号を中間周波数の信号以下、IFIntermediateFrequency信号と称す。にダウンコンバートし、IF信号を増幅等した後、A/D変換器でデジタル信号に変換して、ベースバンド処理回路部20に出力する。ベースバンド処理回路部20は、RF受信回路部11から出力された受信信号に対して、キャリア除去や相関演算等を行ってGPS衛星信号を捕捉する。そして、捕捉したGPS衛星信号から抽出した時刻情報や衛星軌道情報等を利用して、携帯型電話機1の位置及び時計誤差を算出する。ホスト処理部30は、記憶部80に記憶されているシステムプログラム等の各種プログラムに従って携帯型電話機1の各部を統括的に制御するプロセッサーであり、CPUCentralProcessingUnit等のプロセッサーを有して構成される。ホスト処理部30は、ベースバンド処理回路部20から取得した位置座標をもとに、表示部50に現在位置を指し示した地図を表示させたり、その位置座標を各種のアプリケーション処理に利用する。操作部40は、例えばタッチパネルやボタンスイッチ等により構成される入力装置であり、押下されたキーやボタンの信号をホスト処理部30に出力する。この操作部40の操作により、通話要求やメール送受信要求、位置算出要求等の各種指示入力がなされる。表示部50は、LCDLiquidCrystalDisplay等により構成され、ホスト処理部30から入力される表示信号に基づいた各種表示を行う表示装置である。表示部50には、位置表示画面や時刻情報等が表示される。携帯電話用アンテナ60は、携帯型電話機1の通信サービス事業者が設置した無線基地局との間で携帯電話用無線信号の送受信を行うアンテナである。携帯電話用無線通信回路部70は、RF変換回路、ベースバンド処理回路等によって構成される携帯電話の通信回路部であり、携帯電話用無線信号の変調復調等を行うことで、通話やメールの送受信等を実現する。記憶部80は、ROMReadOnlyMemoryやフラッシュROM、RAMRandomAccessMemory等の記憶装置を有して構成され、ホスト処理部30が携帯型電話機1を制御するためのシステムプログラムや、各種アプリケーション処理を実行するための各種プログラムやデータ等を記憶する。時計部90は、携帯型電話機1の内部時計であり、水晶振動子及び発振回路でなる水晶発振器等を有して構成される。時計部90の計時時刻は、ベースバンド処理回路部20及びホスト処理部30に随時出力される。時計部90は、ベースバンド処理回路部20により算出された時計誤差を用いて較正される。3-2.ベースバンド処理回路部の構成図13は、ベースバンド処理回路部20の回路構成及び記憶部23のデータ構成の説明図である。ベースバンド処理回路部20は、主要な機能構成として、処理部21と、記憶部23とを備える。処理部21は、ベースバンド処理回路部20の各機能部を統括的に制御する制御装置及び演算装置であり、CPUやDSPDigitalSignalProcessor等のプロセッサーを有して構成される。処理部21は、衛星捕捉部211と、位置算出部213とを機能部として有する。衛星捕捉部211は、GPS衛星の捕捉を行う機能部である。具体的には、RF受信回路部11から出力されるデジタル化された受信信号に対して、キャリア除去及び相関演算のデジタル信号処理を実行して、GPS衛星を捕捉する。そして、相関演算結果に対するピーク判定を行い、受信キャリア信号のドップラー周波数や受信周波数、受信C/Aコードのコード位相といった情報をメジャメント情報として取得する。また、衛星捕捉部211は、相関演算結果に基づいて航法データを復号する。受信キャリア信号の位相キャリア位相と、受信C/Aコードの位相コード位相とが検出され、相関がとれた状態になると、相関値の時間変化をもとに、航法データを構成する各ビット値をデコードすることができる。この位相同期は、例えば位相ロックループとして知られるPLLPhaseLockedLoop処理により実現される。位置算出部213は、衛星捕捉部211により捕捉された各GPS衛星に係るメジャメント情報や航法データ、時刻情報、衛星情報といった諸量を用いて、上記の原理に従って位置算出演算を行って、携帯型電話機1の位置位置座標及び時計誤差クロックバイアスを算出する。そして、算出した位置をホスト処理部30に出力するとともに、算出した時計誤差で時計部90を較正する。位置算出部213は、GPS受信部10による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部に相当する。記憶部23は、ベースバンド処理回路部20のシステムプログラムや、衛星捕捉機能、位置算出機能といった各種機能を実現するための各種プログラム、データ等を記憶する。また、各種処理の処理中データ、処理結果などを一時的に記憶するワークエリアを有する。記憶部23には、プログラムとして、処理部21により読み出され、ベースバンド処理図14参照として実行される第1のベースバンド処理プログラム231が記憶されている。また、ベースバンド処理プログラム231は、各種の位置算出処理図15~図18参照として実行される位置算出プログラム2311をサブルーチンとして含む。原理で説明した第1~第4の位置算出方法の何れかに対応するプログラムが位置算出プログラム2311として記憶されている。また、記憶部23には、主要なデータとして、衛星軌道データ232と、捕捉した各GPS衛星に係る衛星別データ233と、状態変動正規分布モデル234と、演算結果データ235とが記憶される。衛星軌道データ232は、全てのGPS衛星の概略の衛星軌道情報を記憶したアルマナックや、各GPS衛星それぞれについて詳細な衛星軌道情報を記憶したエフェメリス等のデータである。この衛星軌道データ232は、GPS衛星から受信したGPS衛星信号をデコードすることで取得する他、例えば携帯型電話機1の基地局やアシストサーバーからアシストデータとして取得する。衛星別データ233は、メジャメント情報233Aと、衛星情報233Bと、測位計算用諸量233Cとを含む。メジャメント情報233Aは、当該捕捉衛星について測定したGPS衛星信号のコード位相や受信周波数、ドップラー周波数といった情報である。衛星情報233Bは、当該捕捉衛星の位置や速度、移動方向といった情報である。測位計算用諸量233Cは、上記の原理に従って測位計算を行う際に用いる各種の諸量である。状態変動正規分布モデル234は、図1で図示説明した状態変動正規分布モデルであり、モデルパラメーター値として、例えば、状態変化量期待値234Aと、状態変化量標準偏差234Bとを含む。状態変動正規分布モデルは、演算時刻毎に更新される。演算結果データ235は、位置算出演算の結果として得られたデータであり、携帯型電話機1の位置ベクトルPやクロックバイアスbがこれに含まれる。3-3.処理の流れ図14は、記憶部23に記憶されているベースバンド処理プログラム231が処理部21により読み出されることで、ベースバンド処理回路部20において実行されるベースバンド処理の流れを示すフローチャートである。最初に、処理部21は、記憶部23にサブルーチンとして記憶されている位置算出プログラム2311に従って位置算出処理を開始するステップA1。位置算出処理は、図15~図18の第1~第4の位置算出処理の何れかの処理である。次いで、処理部21は、演算結果の出力タイミングであるか否かを判定するステップA3。演算結果の出力タイミングは、任意のタイミングを設定可能である。例えば、位置算出演算の時間間隔と同じ時間間隔のタイミングとしてもよいし、位置算出演算の時間間隔よりも長い時間間隔のタイミングとしてもよい。また、ユーザーから演算結果の出力が指示されたタイミングを出力タイミングとしてもよい。出力タイミングであると判定した場合はステップA3;Yes、処理部21は、最新の演算結果をホスト処理部30に出力するステップA5。ステップA3において出力タイミングではないと判定した場合ステップA3;No、又は、ステップA5の後、処理部21は、ベースバンド処理を終了するか否かを判定するステップA7。ベースバンド処理を終了しないと判定した場合はステップA7;No、処理部21は、ステップA3に戻る。また、終了すると判定した場合はステップA7;Yes、処理部21は、ベースバンド処理を終了する。図15は、第1の位置算出方法に対応する第1の位置算出処理の流れを示すフローチャートである。最初に、処理部21は、状態ベクトルXの初期値を設定するステップB1。具体的には、サーバーアシストによって携帯型電話機1の基地局から取得した位置基地局位置や、予め定められた位置例えば固定値を、初期位置として設定する。また、クロックバイアスbにも所定の初期値例えば固定値を設定する。次いで、処理部21は、捕捉対象衛星選定処理を行うステップB3。具体的には、時計部90で計時されている現在時刻において、初期位置の天空に位置するGPS衛星を、記憶部23に記憶されたアルマナックやエフェメリス等の衛星軌道データ232を用いて判定して、捕捉対象衛星に選定する。次いで、処理部21は、状態変化量推定処理を行うステップB5~B17。具体的には、処理部21は、各捕捉対象衛星についてループAを実行するステップB5~B13。ループAでは、処理部21は、当該捕捉対象衛星からGPS衛星信号を受信して擬似距離を測定し、その値を擬似距離測定値ρcとするステップB7。擬似距離測定値ρcは、位相サーチの結果として得られるコード位相を用いて算出することができる。また、処理部21は、擬似距離近似値ρaを算出するステップB9。具体的には、衛星軌道データ232から求まる当該捕捉対象衛星の位置ベクトルPkと、携帯型電話機1の最新の位置ベクトルP及び最新のクロックバイアスbとを用いて、擬似距離近似値をρa=∥Pk-P∥bとして算出する。但し、この式では、クロックバイアスbの単位を距離に換算している。その後、処理部21は、擬似距離変化量δρを算出するステップB11。具体的には、ステップB7で測定した擬似距離測定値ρcと、ステップB9で算出した擬似距離近似値ρaとの差を算出して、擬似距離変化量δρとする。そして、処理部21は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップB7~B11を行った後、処理部21は、ループAを終了するステップB13。次いで、処理部21は、幾何行列Gを算出するステップB15。具体的には、衛星軌道データ232から求まる各捕捉対象衛星の衛星位置ベクトルと、携帯型電話機1の最新の位置ベクトルとを用いて、携帯型電話機1から各捕捉対象衛星に向かう視線ベクトルlを計算する。そして、式5の幾何行列Gを求める。その後、処理部21は、擬似距離誤差正規分布モデル及び状態変動正規分布モデルを適用して、状態変化量δXの最尤推定解を求めるステップB17。具体的には、式16に従って状態変化量δXを計算する。擬似距離誤差εに正規分布を仮定するため、ここで得られる最尤推定解は最小二乗解に相当する。そして、処理部21は、状態変化量推定処理を終了する。次いで、処理部21は、ステップB17で最尤推定解として推定した状態変化量ベクトルδXを用いて状態ベクトルXを修正するステップB19。具体的には、最尤推定解として求めた状態変化量ベクトルδXを最新の状態ベクトルXに加算することで、今回の演算時刻における状態ベクトルXを新たに算出する。その後、処理部21は、位置算出を終了するか否かを判定しステップB23、まだ終了しないと判定した場合はステップB23;No、状態変動正規分布モデル更新処理を行うステップB25。具体的には、状態変化量期待値ベクトルμδX及び状態変化量共分散行列SδXを、それぞれ式17及び式18に従って更新する。状態変動正規分布モデル更新処理を行った後、処理部21は、ステップB3に戻り、次の演算時刻での位置算出演算へと移行する。また、ステップB23において位置算出を終了すると判定した場合はステップB23;Yes、処理部21は、第1の位置算出処理を終了する。ステップB3~B25までの一連の処理が、一の演算時刻での位置算出演算に相当する。第1の位置算出処理において特徴的であるのは、各演算時刻での位置算出演算において収束演算が省略されている点である。従来の位置算出方法では、最小二乗法を利用して解を求める際に、例えばニュートン法を利用した収束演算を行って近似解を求めていた。しかし、本実施形態では、演算結果の変動を確率変数とし、その上、演算時刻毎に確率分布モデルを更新することにしているため、収束演算を不要とする。図16は、第2の位置算出方法に対応する第2の位置算出処理の流れを示すフローチャートである。状態変化量ベクトルδXのうち、位置変化量ベクトルδPの成分δx,δy,δzを確率変数とし、クロックバイアス変化量δbを非確率変数として位置算出を行うフローチャートである。なお、図15の第1の位置算出処理と同一のステップについては同一の符号を付して説明を省略する。処理部21は、ステップB3において捕捉対象衛星選定処理を行った後、状態変化量推定処理を行うステップC5。具体的には、処理部21は、第1の位置算出処理の状態変化量推定処理ステップB5~B17と同様の処理を行う。但し、ここではクロックバイアスδbを非確率変数とするため、状態変化量ベクトルδXのクロックバイアス変化量δbに対応する成分をゼロとして計算する。つまり、状態変化量推定処理は、位置変化量ベクトルδPを推定する位置変化量推定処理となる。また、状態変化量推定処理では、1回前の位置算出演算で収束演算によって算出したクロックバイアスbの値を用いて、擬似距離近似値ρaを算出するステップB9。具体的には、衛星軌道データ232から求まる当該捕捉対象衛星の位置ベクトルPkと、前回の位置算出演算で算出した位置ベクトルPと、前回の位置算出演算で収束演算によって算出したクロックバイアスbとを用いて、擬似距離近似値ρa=∥Pk-P∥bを算出する。これは、過去に収束演算によって算出された所定のパラメーターの値を用いて位置算出演算を行うことに相当する。次いで、処理部21は、クロックバイアス変化量推定処理を行うステップC7~C19。クロックバイアス変化量推定処理は、GPS衛星信号の受信信号を用いた所定の収束演算を行って、非確率変数とするクロックバイアス変化量δbの近似解を求める処理である。詳細には、処理部21は、各捕捉対象衛星についてループBを行うステップC7~C13。ループBでは、処理部21は、擬似距離近似値ρaを算出するステップC9。そして、処理部21は、擬似距離変化量δρを算出するステップC11。具体的には、状態変化量推定処理のステップB7で測定した擬似距離測定値ρcと、ステップC9で算出した擬似距離近似値ρaとの差を算出して、擬似距離変化量δρとする。その後、処理部21は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップC9及びC11を行った後、処理部21は、ループBを終了するステップC13。次いで、処理部21は、擬似距離誤差正規分布モデルを適用して、クロックバイアス変化量δbの最小二乗解を求めるステップC17。そして、処理部21は、求めたクロックバイアス変化量δbについて、所定の収束条件が成立したか否かを判定するステップC18。具体的には、例えば、今回の反復ステップで求めたクロックバイアス変化量δbと、1回前の反復ステップで求めたクロックバイアス変化量δbとの差の絶対値が所定の閾値よりも小さくなった場合に、収束条件が成立したと判定する。収束条件が成立していないと判定した場合はステップC18;No、処理部21は、ステップC17で求めたクロックバイアス変化量δbを用いて、クロックバイアスbを修正して更新するステップC19。そして、処理部21は、ステップC7に戻る。また、収束条件が成立したと判定した場合はステップC18;Yes、処理部21は、クロックバイアス変化量推定処理を終了する。その後、処理部21は、状態変化量推定処理において推定した位置変化量ベクトルδPと、クロックバイアス変化量推定処理において推定したクロックバイアス変化量δbとを状態変化量ベクトルδXとし、これを用いて状態ベクトルXを修正するステップC21。つまり、最新の状態ベクトルXに状態変化量ベクトルδXを加算して、今回の演算時刻における状態ベクトルXを新たに算出する。その後、処理部21は、位置算出を終了するか否かを判定しステップB23、まだ終了しないと判定した場合はステップB23;No、状態変動正規分布モデル更新処理を行うステップB25。その後、処理部21は、ステップB3に戻り、次の演算時刻での位置算出演算へと移行する。また、ステップB23において位置算出を終了すると判定した場合はステップB23;Yes、処理部21は、第2の位置算出処理を終了する。この場合も、ステップB3~B25までの一連の処理が、一の演算時刻での位置算出演算に相当する。また、ステップC7~C19の繰り返し処理が、受信信号を用いた所定の収束演算に相当する。第2の位置算出処理において特徴的であるのは、確率変数として取り扱う位置変化量ベクトルδPについては収束演算を省略しているが、非確率変数として取り扱うクロックバイアス変化量δbについては収束演算を省略していない点である。クロックバイアス変化量δbは確率変数とせず、全くの未知数として取り扱うため、収束演算を行わなければ最適解が求まらないためである。図17は、第3の位置算出方法に対応する第3の位置算出処理の流れを示すフローチャートである。第3の位置算出処理の基本的な流れは、図15の第1の位置算出処理と同じである。しかし、第1の位置算出処理のステップB17がステップD17に置き換わる点が大きく異なる。ステップD17では、処理部21は、擬似距離誤差正規混合分布モデル及び状態変動正規分布モデルを適用して、状態変化量ベクトルδXの最尤推定解を求める。具体的には、処理部21は、式24~式27に従って、状態変化量ベクトルδXの最尤推定解を計算する。このステップは、EMアルゴリズムを適用した計算ステップである。また、第3の位置算出処理では、ステップB25の状態変動正規分布モデル更新処理において、処理部21は、状態変化量期待値ベクトルμδXを式28に従って更新する。また、処理部21は、状態変化量共分散行列SδXを式29~式32に従って更新する。図18は、第4の位置算出方法に対応する第4の位置算出処理の流れを示すフローチャートである。第4の位置算出処理の基本的な流れは、図16の第2の位置算出処理と同じである。しかし、第2の位置算出処理のステップC5がステップE5に置き換わり、ステップC17がステップE17に置き換わる点が大きく異なる。第4の位置算出処理では、ステップE5の状態変化量推定処理において、クロックバイアスδbを非確率変数として取り扱い、状態変化量ベクトルδXのクロックバイアス変化量δbに対応する成分をゼロとして計算を行う。つまり、この場合の状態変化量推定処理は、位置変化量ベクトルδPを推定する位置変化量推定処理となる。また、第2の位置算出処理と同様に、状態変化量推定処理では、1回前の位置算出演算で収束演算によって算出されたクロックバイアスbの値を用いて、擬似距離近似値ρaを算出するステップB9。つまり、過去に収束演算によって算出された所定のパラメーターの値を用いて位置算出演算を行う。また、第4の位置算出処理では、ステップC7~ステップC19のクロックバイアス変化量推定処理において、擬似距離誤差正規混合分布モデルを適用して、例えばEMアルゴリズムを利用した収束演算を行ってステップC18;No→ステップC19→ステップC7、クロックバイアス変化量δbの最尤推定解を求めるステップE17。そして、処理部21は、クロックバイアス変化量δbについて所定の収束条件が成立したか否かを判定しステップC18、成立していないと判定した場合はステップC18;No、ステップE17で求めたクロックバイアス変化量δbを用いてクロックバイアスbを修正更新した後ステップC19、ステップC7に戻る。また、収束条件が成立したと判定した場合はステップC18;Yes、処理部21は、ステップC21へと移行する。4.作用効果本実施形態によれば、測位用衛星信号の一種であるGPS衛星信号を受信する。そして、GPS衛星信号の受信信号を用いた演算であって、GPS受信機の位置及びクロックバイアスの変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う。具体的には、GPS受信機位置算出装置の位置及びクロックバイアスの変動を確率変数とする状態変動正規分布モデルを設定する。また、擬似距離に含まれる誤差の変動を確率変数とする確率分布モデルとして、擬似距離誤差正規分布モデル又は擬似距離誤差正規混合分布モデルを設定する。そして、これらの確率モデルを用いて、最尤推定法に基づく演算を行って、位置及びクロックバイアスの変化量の最尤推定解を求める。これにより、演算時刻毎の演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行うことが可能となる。この場合、測定した擬似距離に瞬間的に大きな誤差が混入したとしても、位置算出精度が低下することを防止できる。また、本実施形態では、過去の位置算出演算の結果に基づいて演算結果の変動の基準を更新し、更新した基準を用いて演算結果の変動を算出する。詳細には、1つ前の演算時刻に算出された位置及びクロックバイアスの変化量を用いて、今回の演算時刻での位置及びクロックバイアスの変化量の基準を更新する。これは、状態変動正規分布モデルの期待値を更新することに相当する。また、本実施形態では、1つ前の演算時刻での位置及びクロックバイアスの変化量の共分散を用いて、今回の演算時刻での位置及びクロックバイアスの変化量の共分散を更新する。これは、状態変動正規分布モデルの標準偏差を更新することに相当する。つまり、演算結果の変動に係る正規分布モデルを固定的に設定するのではなく、1つ前の演算時刻の演算結果に基づいて当該正規分布モデルを更新する。これにより、位置及びクロックバイアスの時間変化に俊敏に対応することができるようになり、移動体の位置をより正確に求めることが可能となる。5.変形例5-1.ドップラー測位上記の実施形態では、擬似距離を用いた位置算出方法を例に挙げて説明したが、ドップラーを利用した位置算出方法に対しても、本発明を実質的に同一に適用可能である。ドップラーを利用した位置算出方法はドップラー測位と呼ばれる。ドップラーに関して、次式33のモデル式が成立する。但し、単位を全て速度の次元[m/s]として式を記載している。式33において、hkは、k番目の衛星について測定した周波数のずれ以下、周波数ずれと称す。である。つまり、GPS衛星信号の受信信号の周波数受信周波数と規定された搬送波周波数1.57542[GHz]との差として計算される。受信周波数は、周波数方向の相関演算を行うことで取得することができる。Vkは、k番目の捕捉衛星の速度ベクトルであり、Vは、GPS受信機の速度ベクトルである。従って、VR=Vk-Vは、GPS受信機とk番目の捕捉衛星との相対速度ベクトルである。また、lkは、GPS受信機からk番目の捕捉衛星に向かう視線ベクトルである。相対速度ベクトルVR=Vk-Vを視線ベクトルlkの方向に投影することで、視線方向の相対速度ベクトルドップラーシフトが求まる。また、式33において、dは、GPS受信機の動作クロックのドリフトクロックドリフトである。また、εkhは、k番目の衛星の周波数ずれに内在する観測誤差である。式33に対してニュートン法を適用することを考える。GPS受信機の位置ベクトルP、クロックバイアスb、速度ベクトルV及びクロックドリフトdを成分とする状態ベクトルX=P,b,V,dTを定義する。つまり、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトを状態ベクトルXの成分とする。ここで、式33の右辺-左辺で表される次式34の関数gP,Vを定義する。関数gは、速度ベクトルV及び視線ベクトルlを含む。視線ベクトルlは、位置ベクトルPを変数とする。そのため、関数gは、位置ベクトルP及び速度ベクトルVを変数とする関数である。関数gP,Vの位置ベクトルP、速度ベクトルV、クロックバイアスb及びクロックドリフトdでの偏微分を求めると、それぞれ次式35~38のようになる。ここで、捕捉したK個の衛星について関数gP,Vを纏めた関数GP,Vを便宜的に観測関数と呼び、次式39のように定義する。この場合、次式40の観測方程式が成立する。式40の左辺は、位置ベクトルの初期値P0と、クロックバイアスの初期値b0と、速度ベクトルの初期値V0と、クロックドリフトの初期値d0とにより定まる観測関数の値GP0,b0,V0,d0である。また、式40の右辺において、ベクトルδP,δb,δV,δdTを状態ベクトルXの変動を示す状態変化量ベクトルδXと定義する。また、便宜的に、観測関数Gの各変数での偏微分で定義される行列を偏微分行列Yと定義する。状態変化量ベクトルδXには、各変数の初期値により定まる偏微分行列Yの値が乗算される。このように、式33から、状態変化量ベクトルδX=δP,δb,δV,δdTを未知数とする式40の観測方程式が導出される。この観測方程式を利用して、上記の実施形態と同様に、状態変化量ベクトルδXの最尤推定解を求める。上記の実施形態と異なるのは、位置算出演算の基礎となる観測方程式が、式3から式40に置き換わる点である。この場合、状態変化量ベクトルδX=δP,δb,δV,δdTに含まれる算出パラメーターδx,δy,δz,δb,δvx,δvy,δvz,δdの全部又は一部を確率変数とする確率分布モデルを定義する。位置情報の他に、クロックバイアス、速度情報及びクロックドリフトの変動を確率変数として取り扱う。8つの算出パラメーターの全てを確率変数とする場合は、式8の多次元状態変動正規分布モデルを、8次元状態変動正規分布モデルpδXδXとして定義する。また、8つの算出パラメーターのうち、例えばクロックバイアスb及びクロックドリフトdを除いた6つの算出パラメーターを確率変数とするのであれば、式8の多次元状態変動正規分布モデルを、6次元状態変動正規分布モデルpδXδXとして定義する。なお、式40では未知数が8個であるため、解を求めるためには、少なくとも8個の衛星を捕捉する必要がある。そこで、未知数を減らす工夫をすると実用的である。例えば、速度ベクトルVを固定値例えばゼロとし、位置ベクトルP、クロックバイアスb及びクロックドリフトdを未知数とする。この場合、式40の観測方程式は、次式41のように書き換えることができる。式41では、状態ベクトルXの変動を示す状態変化量ベクトルをδX=δP,δb,δdTと定義している。速度ベクトルV=vx,vy,vzを固定値とすることで、未知数を8個から5個まで減らすことができる。これにより、少なくとも5個の衛星を捕捉することができれば、式41の解を求めることができる。図19は、上記の実施例における処理部21が実行するドップラー測位処理の流れを示すフローチャートである。このドップラー測位処理は、図14のベースバンド処理のステップA1において位置算出処理に置き換えて適用可能な処理である。また、このドップラー測位処理は、式40の観測方程式に基づく処理であり、状態変化量ベクトルδXの全ての算出パラメーターを確率変数とする図15の第1の位置算出処理に対応する処理である。最初に、処理部21は、状態ベクトルXの初期値を設定するステップF1。つまり、位置ベクトルP、クロックバイアスb、速度ベクトルV及びクロックドリフトdに所定の初期値を設定する。次いで、処理部21は、捕捉対象衛星選定処理を行うステップF3。そして、処理部21は、状態変化量推定処理を行うステップF5~F19。具体的には、処理部21は、各捕捉対象衛星についてループCを実行するステップF5~F13。ループCでは、処理部21は、当該捕捉対象衛星からGPS衛星信号を受信して周波数ずれhを算出するステップF7。周波数ずれhは、周波数サーチを行うことで取得した受信信号の周波数受信周波数と、規定された搬送波の周波数規定搬送波周波数との差として算出される。この周波数ずれhは、ドップラー測位における観測量となる。次いで、処理部21は、位置算出装置と当該捕捉対象衛星との相対速度ベクトルVRを算出するステップF9。相対速度ベクトルVRは、位置算出装置の最新の速度ベクトルVと、当該捕捉対象衛星の最新の速度ベクトルVkとを用いて算出される。初回測位の場合は、速度ベクトルの初期値V0を最新の速度ベクトルとし、2回目以降の測位の場合は、前回の測位で算出した速度ベクトルを最新の速度ベクトルとして計算する。その後、処理部21は、位置算出装置から当該捕捉対象衛星に向かう視線ベクトルlを算出するステップF11。視線ベクトルlは、位置算出装置の最新の位置ベクトルPと、当該捕捉対象衛星の最新の位置ベクトルPkとを用いて算出される。初回測位の場合は、位置ベクトルの初期値P0を最新の位置ベクトルとし、2回目以降の測位の場合は、前回の測位で算出した位置ベクトルを最新の位置ベクトルとして計算する。その後、処理部21は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップF7~F11を行った後、処理部21は、ループCを終了するステップF13。ループCの後、処理部21は、最新の位置ベクトルP、クロックバイアスb、速度ベクトルV及びクロックドリフトdを用いて、式40の左辺の観測関数Gの値を算出するステップF15。また、処理部21は、最新の位置ベクトルP、クロックバイアスb、速度ベクトルV及びクロックドリフトdを用いて、式40の右辺の偏微分行列Yの値を算出するステップF17。次いで、処理部21は、式40の観測方程式に従って、観測誤差正規分布モデル及び状態変動正規分布モデルを適用して、状態変化量δXの最尤推定解を求めるステップF19。観測誤差正規分布モデルは、観測量である周波数ずれhに含まれる観測誤差εhを確率変数とし、その分布として正規分布を仮定したモデルである。ステップF19を行った後、処理部21は、状態変化量推定処理を終了する。以降のステップは、第1の位置算出処理と同様である。なお、ここではフローチャートの図示及び説明を省略するが、図16~図18の第2~第4の位置算出処理に対応するドップラー測位処理についても、それぞれ同様に構成することが可能である。5-2.確率変数上記の実施形態では、位置算出演算の演算結果の変動を確率変数とする確率分布モデルを定義したが、位置算出演算の演算結果そのものを確率変数とする確率分布モデルを定義して位置算出演算を行うことも可能である。位置ベクトルPに着目して考えてみる。位置ベクトルPは、位置ベクトルの初期値P0に位置変化量ベクトルδPを加算することで求められる。位置ベクトルの初期値P0は固定値であるため、位置変化量ベクトルδP=P-P0を確率変数として取り扱うことと、位置ベクトルPそのものを確率変数として取り扱うこととは同値である。従って、位置変化量ベクトルδPの代わりに位置ベクトルPそのものを確率変数として位置算出演算を行うことも可能である。クロックバイアスb、速度ベクトルV及びクロックドリフトdについても同様である。この場合は、例えばGPS受信機の状態ベクトルXそのものを確率変数とする状態正規分布モデルを定義する。そして、当該状態正規分布モデルを用いて、上記の実施形態と同様に位置算出演算を行う。状態正規分布モデルの更新では、過去に算出された位置ベクトルPやクロックバイアスbを用いて、状態正規分布モデルの期待値や標準偏差を更新する。5-3.状態変動正規分布モデルの更新上記の実施形態では、1時刻前の位置算出演算の結果に基づいて状態変動正規分布モデルを更新したが、あくまでもこれは一例である。例えば、2時刻前や3時刻前の位置算出演算の結果に基づいて状態変動正規分布モデルを更新してもよい。また、過去所定期間分例えば5時刻分や10時刻分の位置算出演算の結果を平均処理し、その結果に基づいて状態変動正規分布モデルを更新してもよい。平均処理としては、単純な算術平均や幾何平均を適用してもよいし、加重平均を適用してもよい。加重平均を適用する場合は、現在の演算時刻により近い時刻での位置算出演算の結果ほど、重みをより高く設定して加重平均すると効果的である。また、外部から取得可能な参照情報を用いて状態変動正規分布モデルを更新してもよい。例えば、自動車等の移動体に本発明の位置算出装置を設定する場合において、車速検出システム例えば車速パルスや磁気センサー等の方位センサーから、移動体の速度や移動方向つまり速度ベクトルを参照情報として取得可能な場合は、当該参照情報を用いて状態変動正規分布モデルを更新するように構成することも可能である。具体的には、最新の移動体の位置ベクトルPに速度ベクトルVを加味して、現在の演算時刻における移動体の位置ベクトルPを算出する。このようにして位置ベクトルPが算出可能であるため、演算時刻毎の位置ベクトルの変化量δPを算出することができる。従って、速度ベクトルVを利用して位置変化量ベクトルδPを算出し、これを用いて状態変化量期待値ベクトルμδXを更新すればよい。また、状態変動正規分布モデルの標準偏差も同様に更新することが可能である。移動体の進行方向がわかれば、当該進行方向以外の方向に対する速度成分はほとんど無いため、当該方向に対する位置や速度の誤差はゼロに近くなると考えられる。そこで、状態変化量共分散行列SδXのうち、移動体の進行方向以外の方向に対応する共分散の成分を小さな値例えばゼロとするように状態変化量共分散行列SδXを更新すると効果的である。5-4.非確率変数上記の実施形態では、位置算出演算によって算出可能な算出パラメーターのうち、位置ベクトルPの各成分を確率変数とし、クロックバイアスbを非確率変数とする場合を例に挙げて説明した。しかし、これはあくまでも一例であり、確率変数及び非確率変数とする算出パラメーターは適宜に選択可能である。例えば、3次元の位置成分のうちの一部の成分の変動を確率変数とするといったことも可能である。具体的には、位置変化量ベクトルδPのうちの高度方向の位置変化量δzを非確率変数とし、δzについては受信信号を利用した所定の収束演算を行って近似解を求めることとしてもよい。また、高度方向の位置変化量δz及びクロックバイアス変化量δbの2種類の算出パラメーターを非確率変数とし、δz及びδbについては収束演算を行って近似解を求めるといったことも当然に可能である。5-5.最尤推定アルゴリズム上記の実施形態では、最尤推定のためのアルゴリズムとしてEMアルゴリズムを適用した場合を例に挙げて説明したが、他のアルゴリズムを適用してもよいことは勿論である。例えば、山登り法として知られる勾配法を適用してもよい。勾配法は、最適化問題に対するアルゴリズムの一種であり、所定の評価関数の勾配を用いてパラメーターの最適化を行う手法である。5-6.電子機器上記の実施例では、GPS受信機を具備する電子機器として携帯型電話機を例に挙げて説明したが、本発明を適用可能な電子機器はこれに限られるわけではない。例えば、カーナビゲーション装置や携帯型ナビゲーション装置、パソコン、PDAPersonalDigitalAssistant、腕時計といった他の電子機器についても本発明を同様に適用可能である。5-7.衛星測位システムまた、上記の実施形態では、衛星測位システムとしてGPSを例に挙げて説明したが、WAASWideAreaAugmentationSystem、QZSSQuasiZenithSatelliteSystem、GLONASSGLObalNAvigationSatelliteSystem、GALILEO等の他の衛星測位システムであってもよいことは勿論である。5-8.処理の主体上記の実施例では、GPS受信部の内部に設けられた処理部が位置算出処理を行うものとして説明した。つまり、上記の実施例では、GPS受信部ベースバンド処理回路部が位置算出装置として機能するものとして説明した。しかし、位置算出の主体を、電子機器のプロセッサーであるホスト処理部としてもよい。この場合は、例えば、ベースバンド回路部の処理部がGPS衛星信号を捕捉してメジャメント情報を算出取得し、ホスト処理部に出力する。そして、ホスト処理部は、処理部から入力したメジャメント情報を用いて、上記の原理に基づく位置算出演算を行う。この場合は、電子機器ホスト処理部が位置算出装置として機能することになる。"