Anaconda3のインストールが完了したら,Pythonのサンプルプログラム(python_test.ipynb)が動くかどうか確かめてみましょう.
プログラムは以下の手順で実行します.
上記の1~3が上手くいかない場合(図が出ずにエラーメッセージのみが出る場合),Anaconda3のインストールに失敗している可能性があります.再度,以下のホームページから
最新版(2023.07)の Anaconda3-2023.07 をインストールしてください.それでもうまく行かない場合には,以下のお問い合わせ先にご連絡ください.
上記の1~3が上手くいった場合(エラーメッセージは出ずに図が表示される場合),Anaconda3の環境構築はこれで終了です.お疲れさまでした.アメダス気象データ分析チャレンジ!(Python版)の当日にお目にかかれますことを楽しみにしています.
応募者多数により惜しくも抽選で外れてしまった場合には,大変申し訳ございません.当日の資料は後日ホームページ上で公開しますので,この環境を活用して是非,アメダス気象データ分析に挑戦してみてください.
# numpyをnpという別名でインポートします.
import numpy as np
# matplotlibをpltという別名でインポートします.
import matplotlib.pyplot as plt
# Seabornをインポートします.
import seaborn as sns
# pandasをpdという別名でインポートします.
import pandas as pd
# sklearn(scikit-learn)は機械学習関連のライブラリーです.インポートします.
from sklearn import linear_model
%precision 3
%matplotlib inline
"""
4次のルンゲ-クッタ法によるローレンツモデル計算
dx/dt = -s*x + s*y
dy/dt = -y + r*x - x*z
dz/dt = -b*z + x*y
"""
# 1つ目の方程式dx/dt = -s*x + s*yを関数として設定
def f1(t,x,y,z,s):
return -s*x + s*y
# 2つ目の方程式dy/dt = -y + r*x - x*zを関数として設定
def f2(t,x,y,z,r):
return -y + r*x - x*z
# 3つ目の方程式dz/dt = -b*z + x*yを関数として設定
def f3(t,x,y,z,b):
return -b*z + x*y
# パラメータs,r,bの設定
s = 10.0
r = 28.0
b = 8.0/3.0
# データx,y,zの初期値の設定
x = 0.1
y = 0.1
z = 0.1
# 時間刻み
h = 0.01
# t=0からt=100までh刻みで計算します.時間tをnumpy配列で設定.
tvalues = np.arange(0.0,100.0,h)
# データxの保存先をnumpy配列とする
xvalues = []
# データyの保存先をnumpy配列とする
yvalues = []
# データzの保存先をnumpy配列とする
zvalues = []
# ルンゲ-クッタ法で時間積分
for t in tvalues:
xvalues.append(x)
yvalues.append(y)
zvalues.append(z)
k1 = h*f1(t,x,y,z,s)
l1 = h*f2(t,x,y,z,r)
m1 = h*f3(t,x,y,z,b)
k2 = h*f1(t+0.5*h,x+0.5*k1,y+0.5*l1,z+0.5*m1,s)
l2 = h*f2(t+0.5*h,x+0.5*k1,y+0.5*l1,z+0.5*m1,r)
m2 = h*f3(t+0.5*h,x+0.5*k1,y+0.5*l1,z+0.5*m1,b)
k3 = h*f1(t+0.5*h,x+0.5*k2,y+0.5*l2,z+0.5*m2,s)
l3 = h*f2(t+0.5*h,x+0.5*k2,y+0.5*l2,z+0.5*m2,r)
m3 = h*f3(t+0.5*h,x+0.5*k2,y+0.5*l2,z+0.5*m2,b)
k4 = h*f1(t+h,x+k3,y+l3,z+m3,s)
l4 = h*f2(t+h,x+k3,y+l3,z+m3,r)
m4 = h*f3(t+h,x+k3,y+l3,z+m3,b)
x += (k1+2*k2+2*k3+k4)/6
y += (l1+2*l2+2*l3+l4)/6
z += (m1+2*m2+2*m3+m4)/6
# 数値積分おわり
# 後処理
# xの時系列図
print("xの時系列図です")
plt.plot (tvalues, xvalues)
plt.xlabel("t")
plt.ylabel("x(t)")
plt.show()
# yの時系列図
print("yの時系列図です")
plt.plot (tvalues, yvalues)
plt.xlabel("t")
plt.ylabel("y(t)")
plt.show()
# zの時系列図
print("zの時系列図です")
plt.plot (tvalues, zvalues)
plt.xlabel("t")
plt.ylabel("z(t)")
plt.show()
# x-z断面図
print("x-zの断面図です")
plt.plot (xvalues, zvalues)
plt.xlabel("x(t)")
plt.ylabel("z(t)")
plt.show()
# Pandasのデータフレーム形式に変換
print("データフレームに割り当てて,相関係数を計算します")
df = pd.DataFrame({'t':tvalues,'x':xvalues,'y':yvalues,'z':zvalues})
corr=df['x'].corr(df['y'])
# 線形回帰直線
print("xとyの散布図と回帰直線を引きます")
clf = linear_model.LinearRegression()
X=df.loc[:,['x']].values
Y=df.loc[:,['y']].values
clf.fit(X,Y)
plt.plot(X, Y, 'o')
plt.plot(X, clf.predict(X))
plt.xlabel("x(t)")
plt.ylabel("y(t)")
plt.show()
print('xとyの相関係数は',corr)
df
xの時系列図です
yの時系列図です
zの時系列図です
x-zの断面図です
データフレームに割り当てて,相関係数を計算します xとyの散布図と回帰直線を引きます
xとyの相関係数は 0.8867624304628634
t | x | y | z | |
---|---|---|---|---|
0 | 0.00 | 0.100000 | 0.100000 | 0.100000 |
1 | 0.01 | 0.101300 | 0.126889 | 0.097481 |
2 | 0.02 | 0.105051 | 0.154219 | 0.095059 |
3 | 0.03 | 0.111094 | 0.182642 | 0.092737 |
4 | 0.04 | 0.119344 | 0.212771 | 0.090521 |
... | ... | ... | ... | ... |
9995 | 99.95 | 4.805565 | 2.329803 | 26.451100 |
9996 | 99.96 | 4.572819 | 2.392494 | 25.864153 |
9997 | 99.97 | 4.369234 | 2.476375 | 25.290816 |
9998 | 99.98 | 4.193928 | 2.578975 | 24.731966 |
9999 | 99.99 | 4.045893 | 2.698347 | 24.188313 |
10000 rows × 4 columns
上のサンプルのプログラムは,気象学における熱対流(積乱雲も対流の1つ)の問題を表す微分方程式であり,「ローレンツモデル」と呼ばれています(Lorenz, 1963).方程式は以下の通りです.
$\frac{dx}{dt} = - \sigma x + \sigma y $
$\frac{dy}{dt} = - y + r x - xz $
$\frac{dz}{dt} = - b z + x y $
$x$は水平方向の風速の変動,$y$は水平方向の温度の変動,$z$は鉛直方向の温度の変動(のようなもの)を表しています.ルンゲクッタ法という方法で計算しています.余裕がある人はプログラムを解読してみてください.
$x$と$y$と$z$が連動しながら不規則に変化するカオス(chaos) の振る舞いを示します.$x$と$z$の断面図はチョウのような模様(ローレンツバタフライといいます)になっています.初期値(33~35行目)を変えると結果も大きく変わります.このように大気はカオス的性質を持つために,天気予報で用いられる数値予報モデルでは初期値を少しずつ変えたアンサンブル予報が行われています.
話がそれてしまいました...
難しそうな式が出てきてびっくりしたかもしれませんが,アメダス気象データ分析チャレンジ!を受講するにあたって上記を理解している必要は全くありません のでご安心ください.Pythonのプログラムのサンプルとして用意しました.