IT pass HikiWiki - [itbase2024]Fortran 実習 最小二乗法 Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
{{toc_here}}
= 練習問題 1
地球の温度の時間変化のデータを用いて, 移動平均と最小二乗法でのフィッティングを行うプログラムを作りなさい.
((<<これ>|URL:http://itpass.scitec.kobe-u.ac.jp/~itbase/exp/fy2024/gnuplot_data/Earth_TempPRCP.txt>))は, 地球全球の年平均温度偏差と北半球の年平均温度偏差, 南半球の年平均温度偏差のデータである (gnuplot の練習問題に載せているものと同じものです).
ファイルの内容 (数値の意味) は, ファイルの先頭 (ヘッダー) に書いてある.
== 注意
ただし, Fortran ではヘッダーの文字が邪魔になる. 下のような方法で回避することができる.
* ファイルを読み込む部分の最初に,
read( <装置番号>, * )
をヘッダーの行の回数だけ繰り返せばヘッダー部分を読み飛ばすことができる.
* あるいは, Fortran プログラムでファイルの数値を読む前に, エディタを使ってファイルの中のヘッダーを削除してしまうと簡単である.
== 練習問題 1.1
上のファイルを読み込み, 地球全球の年平均温度偏差に 5 年間の移動平均を適用したデータを作りなさい.
また元データと移動平均したデータの両方を gnuplot でプロットし, どのような操作がなされたのか確認しなさい.
=== 移動平均
移動平均とは, 区切られた期間ごとの平均のことである.
例えば, 一年ごとの温度データを, T1, T2, T3, T4, T5, T6, T7, T8, ... とする.
このデータの 3 年間の移動平均とは,
(T1+T2+T3)/3, (T2+T3+T4)/3, (T3+T4+T5)/3, ...
のデータとなる.
こうすることで, 長い時間スケールの温度変化が見やすくなる.
== 練習問題 1.2
上のファイルを読み込み, 地球全球の年平均温度偏差に 10 年間の移動平均を適用したデータを作りなさい.
また元データと 5 年間, 10 年間の移動平均を適用したデータの両方を gnuplot でプロットし, どのような操作がなされたのか確認しなさい.
== 練習問題 1.3
元データを使って (移動平均したデータを使わないこと), 最小二乗法に基づいて地球全球の年平均温度偏差の時間変化に最もよく合う一次関数を求めなさい.
そして, 元データと求めた一次関数を gnuplot でプロットし, 一次関数が合っているように見えるかどうか確認しなさい.
#なお, 上に指示したファイルに保存されている個々の温度の値には, 標準偏差の値は含まれていない. 今回はすべての温度の値の標準偏差は 1 (K) として最小二乗法を使いなさい.
= 練習問題 2
二酸化炭素濃度の時間変化のデータを用いて, 平均的な濃度の時間変化率 (ppm/年; 1 年あたりの濃度変化率) を最小二乗法に基づいて求めなさい.
なお, ppm は parts per million の略であり, 百万分の一を表す.
((<<これ>|URL:http://itpass.scitec.kobe-u.ac.jp/~itbase/exp/fy2024/average_dispersion/co2_monthave_ryo.txt>))は, 綾里 (岩手県) で観測された月平均二酸化炭素濃度のデータである.
ファイルの内容 (数値の意味) は, ファイルの先頭 (ヘッダー) に書いてある.
なお, 濃度の値の中には負の値が含まれている. 負の値は計算からは除くこと.
#また, 上に指示したファイルに保存されている個々の濃度の値には, 標準偏差の値は含まれていない. 最小二乗法を用いる際には, 標準偏差を 1 (ppm) とするか, 各月ごとの平均値と標準偏差を計算し, それらを用いて最小二乗法に用いると良い.
== 注意
ただし, Fortran ではヘッダーの文字が邪魔になる. 下のような方法で回避することができる.
* ファイルを読み込む部分の最初に,
read( <装置番号>, * )
をヘッダーの行の回数だけ繰り返せばヘッダー部分を読み飛ばすことができる.
* あるいは, Fortran プログラムでファイルの数値を読む前に, エディタを使ってファイルの中のヘッダーを削除してしまうと簡単である.
= 練習問題 3
神戸気象台における日平均温度のデータを用いて下のことを行いなさい.
((<<これ>|URL:http://itpass.scitec.kobe-u.ac.jp/~itbase/exp/fy2024/average_dispersion/Kobe_Temp.txt>))は, 神戸気象台で観測された 10 年間の日平均温度のデータである.
ファイルの内容 (数値の意味) は, ファイルの先頭 (ヘッダー) に書いてある.
== 注意
ただし, Fortran ではヘッダーの文字が邪魔になる. 下のような方法で回避することができる.
* ファイルを読み込む部分の最初に,
read( <装置番号>, * )
をヘッダーの行の回数だけ繰り返せばヘッダー部分を読み飛ばすことができる.
* あるいは, Fortran プログラムでファイルの数値を読む前に, エディタを使ってファイルの中のヘッダーを削除してしまうと簡単である.
上のデータにおける「日付」には, "/" (スラッシュ) が含まれるため, 文字型変数を使って読むと良いだろう.
== 練習問題 3.1
上のファイルを読み込み, 月平均温度の時間変化のデータを作りなさい.
また元データと移動平均したデータの両方を gnuplot でプロットし, どのような操作がなされたのか確認しなさい.
== 練習問題 3.2
最小二乗法に基づいて, 3.1 で作成した月平均データに合う, 1 年周期の余弦関数 (cosine) を求めなさい. それにより, 1 年周期の温度変化の振幅と位相 (例えば温度が最大となる日付) を求めなさい.
=== 補足 1
余弦関数は, 一般には,
f(t) = A + B cos( omega * t - Phi )
と書ける. ここで, A, B, Phi は定数であり, omega は周波数である.
(一年周期であれば, omega = 2 pi / (365 日).)
A, B, Phi を最小二乗法を用いて決定すればよい.
しかし, 上の形の f(t) は Phi に関する非線型関数であり, 本実習で説明した方法
のみでは解くことができない.
そこで,
f(t) = A + B cos( omega * t - Phi )
= A + B cos( omega * t ) * cos( Phi ) - B sin( omega * t ) * sin( Phi )
と変形し,
C = B cos( Phi )
D = B sin( Phi )
と置くことで,
f(t) = A + C cos( omega * t ) - D sin( omega * t )
と書けば, f(t) は A, C, D に対する線形関数となるため, 容易に最小二乗法を適用できる.
これにより求めた C, D から, B, Phi を求めれば, 振幅と位相が求められる.
なお, この説明では,
f(t) = A + B cos( omega * t - Phi )
と書いたが, これは,
f(t) = A cos( 0 * t - Phi0 ) + B cos( omega * t - Phi )
とも書けるため, 定数 A は周波数ゼロの成分の振幅である.
=== 補足 2
振幅と位相を求める別の方法はフーリエ変換である.
本実習では数値積分には触れないため, (離散)フーリエ変換の方法は説明しない.
興味があれば調べてみると良いだろう.
なお, 最小二乗法と(離散)フーリエ変換は実はやっていることはほとんど同じである.
= 練習問題 1
地球の温度の時間変化のデータを用いて, 移動平均と最小二乗法でのフィッティングを行うプログラムを作りなさい.
((<<これ>|URL:http://itpass.scitec.kobe-u.ac.jp/~itbase/exp/fy2024/gnuplot_data/Earth_TempPRCP.txt>))は, 地球全球の年平均温度偏差と北半球の年平均温度偏差, 南半球の年平均温度偏差のデータである (gnuplot の練習問題に載せているものと同じものです).
ファイルの内容 (数値の意味) は, ファイルの先頭 (ヘッダー) に書いてある.
== 注意
ただし, Fortran ではヘッダーの文字が邪魔になる. 下のような方法で回避することができる.
* ファイルを読み込む部分の最初に,
read( <装置番号>, * )
をヘッダーの行の回数だけ繰り返せばヘッダー部分を読み飛ばすことができる.
* あるいは, Fortran プログラムでファイルの数値を読む前に, エディタを使ってファイルの中のヘッダーを削除してしまうと簡単である.
== 練習問題 1.1
上のファイルを読み込み, 地球全球の年平均温度偏差に 5 年間の移動平均を適用したデータを作りなさい.
また元データと移動平均したデータの両方を gnuplot でプロットし, どのような操作がなされたのか確認しなさい.
=== 移動平均
移動平均とは, 区切られた期間ごとの平均のことである.
例えば, 一年ごとの温度データを, T1, T2, T3, T4, T5, T6, T7, T8, ... とする.
このデータの 3 年間の移動平均とは,
(T1+T2+T3)/3, (T2+T3+T4)/3, (T3+T4+T5)/3, ...
のデータとなる.
こうすることで, 長い時間スケールの温度変化が見やすくなる.
== 練習問題 1.2
上のファイルを読み込み, 地球全球の年平均温度偏差に 10 年間の移動平均を適用したデータを作りなさい.
また元データと 5 年間, 10 年間の移動平均を適用したデータの両方を gnuplot でプロットし, どのような操作がなされたのか確認しなさい.
== 練習問題 1.3
元データを使って (移動平均したデータを使わないこと), 最小二乗法に基づいて地球全球の年平均温度偏差の時間変化に最もよく合う一次関数を求めなさい.
そして, 元データと求めた一次関数を gnuplot でプロットし, 一次関数が合っているように見えるかどうか確認しなさい.
#なお, 上に指示したファイルに保存されている個々の温度の値には, 標準偏差の値は含まれていない. 今回はすべての温度の値の標準偏差は 1 (K) として最小二乗法を使いなさい.
= 練習問題 2
二酸化炭素濃度の時間変化のデータを用いて, 平均的な濃度の時間変化率 (ppm/年; 1 年あたりの濃度変化率) を最小二乗法に基づいて求めなさい.
なお, ppm は parts per million の略であり, 百万分の一を表す.
((<<これ>|URL:http://itpass.scitec.kobe-u.ac.jp/~itbase/exp/fy2024/average_dispersion/co2_monthave_ryo.txt>))は, 綾里 (岩手県) で観測された月平均二酸化炭素濃度のデータである.
ファイルの内容 (数値の意味) は, ファイルの先頭 (ヘッダー) に書いてある.
なお, 濃度の値の中には負の値が含まれている. 負の値は計算からは除くこと.
#また, 上に指示したファイルに保存されている個々の濃度の値には, 標準偏差の値は含まれていない. 最小二乗法を用いる際には, 標準偏差を 1 (ppm) とするか, 各月ごとの平均値と標準偏差を計算し, それらを用いて最小二乗法に用いると良い.
== 注意
ただし, Fortran ではヘッダーの文字が邪魔になる. 下のような方法で回避することができる.
* ファイルを読み込む部分の最初に,
read( <装置番号>, * )
をヘッダーの行の回数だけ繰り返せばヘッダー部分を読み飛ばすことができる.
* あるいは, Fortran プログラムでファイルの数値を読む前に, エディタを使ってファイルの中のヘッダーを削除してしまうと簡単である.
= 練習問題 3
神戸気象台における日平均温度のデータを用いて下のことを行いなさい.
((<<これ>|URL:http://itpass.scitec.kobe-u.ac.jp/~itbase/exp/fy2024/average_dispersion/Kobe_Temp.txt>))は, 神戸気象台で観測された 10 年間の日平均温度のデータである.
ファイルの内容 (数値の意味) は, ファイルの先頭 (ヘッダー) に書いてある.
== 注意
ただし, Fortran ではヘッダーの文字が邪魔になる. 下のような方法で回避することができる.
* ファイルを読み込む部分の最初に,
read( <装置番号>, * )
をヘッダーの行の回数だけ繰り返せばヘッダー部分を読み飛ばすことができる.
* あるいは, Fortran プログラムでファイルの数値を読む前に, エディタを使ってファイルの中のヘッダーを削除してしまうと簡単である.
上のデータにおける「日付」には, "/" (スラッシュ) が含まれるため, 文字型変数を使って読むと良いだろう.
== 練習問題 3.1
上のファイルを読み込み, 月平均温度の時間変化のデータを作りなさい.
また元データと移動平均したデータの両方を gnuplot でプロットし, どのような操作がなされたのか確認しなさい.
== 練習問題 3.2
最小二乗法に基づいて, 3.1 で作成した月平均データに合う, 1 年周期の余弦関数 (cosine) を求めなさい. それにより, 1 年周期の温度変化の振幅と位相 (例えば温度が最大となる日付) を求めなさい.
=== 補足 1
余弦関数は, 一般には,
f(t) = A + B cos( omega * t - Phi )
と書ける. ここで, A, B, Phi は定数であり, omega は周波数である.
(一年周期であれば, omega = 2 pi / (365 日).)
A, B, Phi を最小二乗法を用いて決定すればよい.
しかし, 上の形の f(t) は Phi に関する非線型関数であり, 本実習で説明した方法
のみでは解くことができない.
そこで,
f(t) = A + B cos( omega * t - Phi )
= A + B cos( omega * t ) * cos( Phi ) - B sin( omega * t ) * sin( Phi )
と変形し,
C = B cos( Phi )
D = B sin( Phi )
と置くことで,
f(t) = A + C cos( omega * t ) - D sin( omega * t )
と書けば, f(t) は A, C, D に対する線形関数となるため, 容易に最小二乗法を適用できる.
これにより求めた C, D から, B, Phi を求めれば, 振幅と位相が求められる.
なお, この説明では,
f(t) = A + B cos( omega * t - Phi )
と書いたが, これは,
f(t) = A cos( 0 * t - Phi0 ) + B cos( omega * t - Phi )
とも書けるため, 定数 A は周波数ゼロの成分の振幅である.
=== 補足 2
振幅と位相を求める別の方法はフーリエ変換である.
本実習では数値積分には触れないため, (離散)フーリエ変換の方法は説明しない.
興味があれば調べてみると良いだろう.
なお, 最小二乗法と(離散)フーリエ変換は実はやっていることはほとんど同じである.