気象観測データに1年周期のサインカーブをあてはめる (Rによる暫定版)

(2021年度秋学期の教材ページはこれからつくる予定だが、暫定的に 2020年秋学期の教材)[海陸分布がもたらす気候の特徴]の中で、 「気温の年変化の振幅と位相」をあつかった。 緯度などの条件が同じならば、海の影響をうけるほど、気温の年変化の振幅は小さくなり、 年変化の位相は遅いほうにずれるだろう。

わたしは、1998年ごろに、日本をふくむ世界の多数の地点の毎日の気温の観測値に、 1年周期のサインカーブ (正弦波) をあてはめて、振幅や位相の分布をみた。 結果を研究論文にするところまでいかなかったが、 その下書きのようなものを[気温の年周期の振幅と位相 -- 両半球中高緯度、地上気温の日データによる]のページに 置いている。 そのときは、最小2乗法によるあてはめを、Fortranでプログラムを書いて実行した。 図は GMT (Generic Mapping Tools)というソフトウェアでかいた。

今回、その計算の手順を基本的に変えないで、R でかきなおしてみた。 まだ、計算手順(とくに単位のついた数量と無次元の数量のあいだの換算など)についても、R によるプログラミングについても、 うまくないところがあるので、さらに改良したいと考えている。

対象となるデータについても、もっと汎用化することがのぞましいが、 いまの段階では、[気象庁の観測データを使う (序論)]で例としてつかった アメダス 館林 の 2018年12月から2020年11月までの毎日の日平均気温をそのまま使ってみることにする。

Rによるプログラム例

• sinfit_sjis.R (日本語文字コード Shift_JIS)
• sinfit_utf8.R (日本語文字コード UTF-8)

このプログラムを実行してつくった図

この例では、気温の平均値は 15.8 ℃、1年周期の振幅は 11.2 ℃、 位相は1周を1として0.577、 それを1年周期成分が最大となる日付になおすと 1月1日を1として211日、つまり7月30日となった。 1年周期成分の分散は、日ごとの観測値の集まりがもつ分散の 88.9 % をしめている。

(毎日の値を見ると、実際の気温の極大は7月30日よりだいぶあとに出ている。 数値をみると、最大値は2020年8月11日の32.4℃だ。 この2年間、7月はまだ梅雨の状態で気温があまりあがらず、8月になってから暑い日が多くなったのだ。 天候には年々変動があって、この特徴はかならずしも常にみられるものではないが、 1年周期のサインカーブからのこのようなずれかたは、ありがちなことだ。)

何を計算したかの説明

最小2乗法

最小2乗法の考えかたは、つぎのようなものである(と、わたしは理解している)。

• 変量 y は、変量 x の関数 f(x) に 誤差 をくわえたものである (と認識する)。
• f(x) の関数の式の形はわかっている(とする)が、その係数の値は、x と y の観測値の組に なるべくよく合うように決める。
• 「なるべくよく合う」とは、{ y_j - f(x_j) }² を すべての j について合計したもの (残差の2乗和) が最小になることだとする。
• (観測値の誤差が一様でないばあいは、誤差を考慮した重みづけをすることもあるが、 ここでは重みづけしない。観測値どうしで誤差の性質は同様であると仮定したことになる。)
• (誤差の確率分布が正規分布であるばあい、最小2乗法は最尤[さいゆう]推定法でもある。 しかし、最小2乗法は、誤差が正規分布しないばあいにも、最尤でないことを承知で、つかわれる。)
• (データに大きくはずれた値がふくまれているばあい、 最小2乗法では、観測値の組にうまくあてはまる関数にたどりつかないこともある。)

ここでの最小2乗法の利用

• 上の一般論での x にあたるのは、時間 t である。 ただし、ちがう年を区別せず、年のうちの相対的な日付だけを問題にする。 1月1日を1としたとおしの日付を、平年は365、うるう年は366でわって無次元化したものを t とする。 結果の1年周期成分の位相をしめすときだけ、365をかけて「日」単位にしてしめす。
• 上の一般論での y にあたるのは、気温 (日平均気温) である。変数名は y にしておく。 観測値が欠損とされているばあいのほか、不完全とされているもの(「品質情報」が「8」でないもの)も 計算につかわないことにする。 (曲線のあてはめかたのうちにはデータ欠損があるとうまくいかないものもあるが、 最小2乗法は欠損は欠損としてあつかえる。)
• あてはめる関数は、定数項(年平均値にあたる)と、 1年周期、(1/2)年周期、(1/3)年周期 ... の sin 成分、cos 成分からなるものとする。 最小2乗で値をきめるべき係数は、a₀ と a_k, b_k (kは1からnまで)である。 T は時間の次元をもつ定数で、1年にあたる値をもつ。 【つぎの式は、Σの下と上に文字をうまくかけなかったが、「kを1からnまで変えた総和」のつもりである。】
  • f(t) = a₀ + Σ_{k = 1..n} {a_k cos(2π k t/T) + b_k sin(2π k t/T)}
• ここでは、1年周期成分だけがほしいので、n = 1 で計算する。 【なお、1998年には、わたしは、n = 6 で計算した。つまり、(1/6)年周期成分まで計算した。 しかし、結果を議論した対象はほぼ1年周期成分だけで、すこしだけ半年周期成分にもふれた。】
  • f(t) = a₀ + a₁ cos(2π t/T) + b₁ sin(2π t/T)
• 係数をきめる問題は、連立1次方程式になる。 R には連立1次方程式をとく関数 solve があるので、それを利用した。 【1998年には、「ガウスの消去法」を、Fortran言語で、 数値計算の教科書の例文にあった Fortran サブルーチンを利用して、実行した。】

結果としてしめした数値

• 平均値 a₀
• 1年周期の振幅 sqrt(a₁² + b₁²), ただし sqrt は平方根をあらわす。
• 1年周期の位相 arctan(b₁ / a₁), arctan は tan の逆関数。ただし a₁ が負のときは πラジアン (=180度)ずらす。 これを、2πでわって、1年を1とする相対的なわりあいとしてしめす。 さらに、365をかけて、年のうちの日付としてしめす。
• 1年周期成分 a₁ cos(2π t/T) + b₁ sin(2π t/T) の分散の yのデータの分散にしめる割合。

グラフにしめしたもの

• (t, y)のデータ値。小さな記号(〇)で。
• 平均値 a₀。横線で。
• 1年周期成分が最大値をとる日付。縦線で。
• 1年周期成分の時系列。赤の曲線で。 (Rには関数のグラフをかく機能があるのだが、 わたしはまだそれをつかいこなせず、1日きざみの折れ線で近似してグラフにしている。)