[ 岡山大学 | 理学部 | 地球科学科 | 地球および惑星大気科学研究室 ]
- 全部の問題を解かなければならないというものではありません．
- とばしてもよいと思う問題には「＃」を付けています．
- 難しいと思う問題には「％」を付けています．
- かなり難しい問題も混ざっています．
- できる人はいろいろやってみてください．
- それ以外の人は，できる範囲でやってみてください．
Skymonitor の稼働率を見る
- 以下の３つのファイルには，Skymonitorが稼働し始めた2016年10月から2021年9月末までに撮影した写真について，各年に撮影した枚数，各月に撮影した枚数，各時間帯(0時台，1時台，．．．23時台)に撮影した枚数，が書かれている
/home/atmos/ipesc/sample/gnuplot/year.dat
/home/atmos/ipesc/sample/gnuplot/month.dat
/home/atmos/ipesc/sample/gnuplot/hour.dat
- 各年に撮影された枚数を図にして見る
- 各月に撮影された枚数を図にして見る
- 各時間帯に撮影された枚数を図にして見る
- まずはファイルの中身をみて，データの書式(データは何列あって，各列に何の数値が入っているのか)を確認する
- 図は見やすくなるように工夫する
- 点の種類を変えてみる，例えば pointtype を使う
- 点の大きさを変えてみる，例えば pointsize を使う
- 年ごとのデータ
- 1440 で割った値を表示してみる
- まる１日稼働すると1440枚の写真が撮影される，すなわち1440で割ると稼働日数になる
- 2017年と2020年は稼働率が高い，2021年も9月末までであることを考えると高い稼働率である
- 月ごとのデータ
- 1440 で割った値を表示してみる
- まる１日稼働すると1440枚の写真が撮影される，すなわち1440で割ると稼働日数になる
- 1440 で割ったものをさらに 5 で割る
- 5 はスカイモニターが動き始めてからの(おおよそ)年数
- スカイモニターが定常運用を開始したのは2017年1月で，2016年のデータは少ない(ほとんどない)
- (おおよその)年数で割れば，(おおよその)１年あたりの稼働日数になる
- ２月と８月は欠測が少ない
- ８月は31日，２月はうるう年が１回あるので28.2日
- ５月と６月の欠測が多いのは，2019年の５月と６月にスカイモニターの運用を停止していたため
- スカイモニターのカメラを他の用途で使っていた
- ５月と６月を 5 でなく 4 で割れば，2017-2018 と 2020 の３年間の稼働日数(１年あたり)になる
- 時間帯ごとのデータ
- 60 で割った値を表示してみる
- まる１時間稼働すると60枚の写真が撮影される，すなわち60で割ると稼働時間数になる
- 各時間帯は１日につき１時間なので，稼働時間数＝稼働日数である
- 60 で割ったものをさらに 5 で割る
- (おおよその)年数で割れば，(おおよその)１年あたりの稼働日数になる
- 最も稼働率が低いのは15時台(15:00-15:59)，最も稼働率が高いのは21時台(21:00-21:59)
- 時間帯によって稼働率が異なる理由ははっきりとはわからないが，おそらく以下の２つで説明可能ではないかと思う
- スカイモニターが不具合を起こして停止する時間帯は特に決まっていない
- 不具合を修正して再稼働を始める時間は夕方であることが多い
文字の出現頻度を比べる
- 以下の６つのファイルには，４つのテキストにおいて各文字が使用された回数が書かれている
/home/atmos/ipesc/sample/text/char01.dat
/home/atmos/ipesc/sample/text/char02.dat
/home/atmos/ipesc/sample/text/char03.dat
/home/atmos/ipesc/sample/text/char04.dat
/home/atmos/ipesc/sample/text/char05.dat
/home/atmos/ipesc/sample/text/char06.dat
- ６つのテキストにおける文字の出現頻度を比較する図を作る
- まず６つのファイルの中身を見る
- ファイルの中身を見るには cat とか less とかを使う
- emacs で開いてもよい
- 図を作るとき，使いにくいものは無理に使わなくてもよい
- 1列目は数字じゃないので使いにくい
- 図にするだけなら0列目を使えばよい
- ６つのテキストは文字数が異なるので，使用回数そのままを図にしても，わかりやすい図にならない．文字数で規格化して図にするのがよい．
- ファイル datafile の2列目にある数字の和は，以下で計算することができる
$ awk '{print $2}' datafile | xargs | tr ' ' '+' | bc
- 文字の出現頻度からテキストを推定する
- ６つのテキストは
/home/atmos/ipesc/sample/text/Close_to_you.txt
/home/atmos/ipesc/sample/text/Hotel_Mauna_Kea.txt
/home/atmos/ipesc/sample/text/Lorem_ipsum.txt
/home/atmos/ipesc/sample/text/My_favorite_things.txt
/home/atmos/ipesc/sample/text/Potential_Vorticity.txt
/home/atmos/ipesc/sample/text/Pangram.txt
Lorenzアトラクタを描いてみる
- Lorenz方程式を数値積分した結果 /home/atmos/ipesc/sample/chaos/lorenz.dat
- データは 4 列で，左から t, x, y, z
- x-y 平面，x-z 平面，y-z 平面，などで図を描いてみる
- t-x 平面，t-y 平面，t-z 平面，などで図を描いてみる
- 初期値をちょっとだけ変えてLorenz方程式を積分した結果 /home/atmos/ipesc/sample/chaos/lorenz-dx.dat と比較してみる
- 点が大きいと重なってしまって見えにくくなるので，サイズを小さくしてみる．たとえば pointsize 0.1
- あるいは点で描くのをやめて線で描いてもいいかも
Edward N. Lorenzが大気の流れを表す方程式を簡略化して導いた連立方程式．
dx/dt = P ( y - x )
dy/dt = - y - x z + R x
dz/dt = x y - b z
決定論的な常微分方程式であるが，カオス的な振る舞いをする．
- 大雑把に分けて，x>0 でぐるぐるまわる場合と，x<0 でぐるぐるまわる場合，２通りあるのがわかる．
- x>0 でぐるぐるまわるのと，x<0 でぐるぐるまわる，この間の遷移に規則性があるように見えない(非周期性)．
- 初期値がほんのちょっと違う場合，しばらく経つとまったく違う状態になってしまう(初期値鋭敏性)．
- lorenz.dat と lorenz-dx.dat の１行目を比べれば，初期値にどれくらいの差があるか確認できる．
- 予測不可能性とも呼ばれる．ある期間より先の天気予報ができない原因．
gnuplot> splot '/home/atmos/ipesc/sample/chaos/lorenz.dat' using 2:3:4 with lines
とすると，xyz を３次元でプロットする．
図の上にカーソルを持っていって，クリックしたままカーソルを動かすと，視点を変えることができる．
CO2濃度
- マウナロアで測定されたCO2濃度のデータ（ソース）
/home/atmos/ipesc/sample/co2/insitu.txt
- まずはデータを見る
- シェルで cat や less を使う
- 先頭から151行目まではデータについての説明が書いてある
- データの本体があるのは152行目から
- gnuplot は行頭が # で始まる行を注釈とみなして無視する
- 先頭から始まるデータの説明が書かれた部分について特別なことをしないでもよい
- CO2濃度は9列目，観測時刻を年で表した数値が8列目にあるので，時系列を表示するなら
gnuplot> plot '/home/atmos/ipesc/sample/co2/insitu.txt' using 8:9 with lines
- 欠測したもの(-999.99)は表示しないようにする
gnuplot> set datafile missing '-999.99'
gnuplot> plot '/home/atmos/ipesc/sample/co2/insitu.txt' using 8:9 with lines
円周率に出てくる数字の出現頻度を見る．
- 円周率100000桁に出てくる各数字の出現頻度は「シェルスクリプト練習問題２−５」で数えている
- 数えていないという人は，出現頻度を数えた結果が以下のファイルにあるので，それを使う
/home/atmos/ipesc/sample/number/pi.dat
- ランダムであるなら，各数字の出現確率は等しく 1 / 10 = 0.1
- 各数字の出現頻度の期待値は 100000 x 0.1 = 10000
出現確率が 0.1 であった場合の出現頻度の95％信頼区間を重ね描きする．
- n 回試行したとき，確率 p の事象が生じる回数の期待値は
n p
- n 回の試行で np が観測される，というわけではない
- サイコロを6回振って1の目が1回出ることもあるが，それ以外(0回だったり，2回だったり．．．)のこともある
- n が大きいとき，観測される回数の分布は np を中心とした正規分布になる
[ n p - 1.96 ( n p ( 1 - p ) )^0.5 , n p + 1.96 ( n p ( 1 - p ) )^0.5 ]
- 横軸に平行な直線は，数字そのままでプロットできる．例えば y=2 であれば
gnuplot> plot 2
- もし数字がランダムであるなら，p=0.1 であり，n=100000 であるなら，観測される回数は95％の確率で以下の範囲にはいる
[ n p - 1.96 ( n p ( 1 - p ) )^0.5 , n p + 1.96 ( n p ( 1 - p ) )^0.5 ]
- いずれの数字も95％の範囲に入っているので，各数字の出現頻度はランダムである場合に予想されるばらつきの範囲内にあると言える
- 予想されるばらつきの範囲内にあるからといって，「ランダムである」と結論することはできない
- もし仮に95％信頼区間外の頻度が観測されていたとしても，それをもって「ランダムではない」ということにはならない
- とりあえずWikipedia
- 信頼区間と信頼度
- 1.64 sigma 90%
- 1.96 sigma 95%
- 2.58 sigma 99%
- 3.29 sigma 99.9%
各数字の出現確率を見る．
- ２列目を100000で割ったらよい，簡単にやるなら
using 1:($2/100000.0)
各数字の出現確率を区間推定し，出現確率と95％信頼区間を表示する．
- N個の中にn個あったなら，その出現確率の点推定は n/N
- 95%信頼区間は，p = n/N として，
[ p - 1.96 ( p ( 1 - p ) / N )^0.5 , p + 1.96 ( p ( 1 - p ) / N )^0.5 ]
- using でデータを３列分指定し，with errorbars を使うと３列目を縦方向の誤差の大きさとして表示してくれる
using 1:2:3 with errorbars
- 区間の大きさを計算してデータの３列目として書き加えてもよいけど，２列目から計算するのが簡単だと思う
using 1:($2/100000.0):(1.96*(($2/100000.0)*(1.0-$2/100000.0)/100000.0)**(0.5)) with errorbars
N=100000.0
using 1:($2/N):(1.96*(($2/N)*(1.0-$2/N)/N)**(0.5)) with errorbars
- 区間推定の意味は，95％の確率でこの区間に入る，ということ
- 区間推定された範囲外にある確率は5％で，大きい可能性と小さい可能性がそれぞれ2.5％ずつ
- 信頼区間の幅は N の平方根に反比例する
- N を大きくすれば，出現確率をより狭い範囲で推定することができる
- 100桁だけとか1000桁だけとかで同じことをやってみて比べてみるといいかも
Last Updated: 2021/11/11, Since: 2019/11/07.
This page is generated by Makefile.rd2html.