両対数グラフの上で等間隔な点を打つ

細かい話ですが、何度も計算し直すのが面倒になったのでここでまとめておきます。

設定

あるパラメータ$\theta \in [\theta_{start} , \theta_{end}]$を動かして$\theta$の値ごとにシミュレーションを行い、その結果得られた観測値$f(\theta)$をグラフ上にプロットしたいものとします。このとき$f(\theta) \propto (\theta-\theta_c)^{\alpha}$のような依存性があり、この$\alpha$を求めたいものとします。$\theta_c$はあらかじめ求まっている、あるいは何かしらの値を一時的に仮定するものとします。具体的には臨界指数などを求めたい状況です。
このとき$f(\theta)$と$\theta-\theta_c$を両対数グラフ上にプロットすると、その傾きが求めたい指数$\alpha$を表しています。そこで、素朴に$\theta_{start}$から$\theta_{end}$まで等間隔に$N$個の$\theta$の値を選び、そのグラフをプロットしてみましょう。

このように、確かに直線が得られているため$f(\theta) \propto (\theta-\theta_c)^{\alpha}$という関係は成り立っていそうです。ただ、$\theta-\theta_c$が等間隔であるために、両対数グラフではグラフの一部分が密集してしまいます。これを両対数上でも等間隔になるように$\theta$をサンプルしたいというのがモチベーションです。特に臨界指数を計算するにあたっては$\theta-\theta_c$が十分近いときの情報が重要になるため、グラフの左端の部分をより細かく刻みたいというのは自然な発想です。

パラメータリストの決定

いま求める$N$個のパラメータのリスト$\theta_i$が満たすべき条件は以下のようになります。
まず$\theta_0 = \theta_{start}$と$\theta_{N-1} = \theta_{end}$です。また両対数グラフ上で$\theta_i-\theta_c$が等間隔になるため、ある値$r$が存在して以下の式が成立します。
$$
\theta_{i+1} - \theta_c = r (\theta_{i} - \theta_c)
$$
よってこの式を繰り返し適用すると$\theta_{N-1} - \theta_c = r^{N-1} (\theta_0 - \theta_c)$が得られるため、両端の条件を考慮すると$r$の満たすべき式は
$$
r = \left( \dfrac{\theta_{end}-\theta_c}{\theta_{start}-\theta_c} \right)^{1/(N-1)}
$$
となります。
以上から、パラメータ$\theta$の選び方は次のような手順で示すことが出来ます。

1. 考えているモデルから$\theta_c$を適当に決める(既知の場合は理論値を用いる、わからないときは一時的に適当に決め、べき則に従うようなものを採用する)。
2. パラメータの数$N$と$\theta$の範囲$\theta_{start}$と$\theta_{end}$を決める。
3. $r$を上の式に基づいて求める。
4. $\theta_i = \theta_c + r^i (\theta_{start}-\theta_c)$を順に計算する。

この方法に基づいて先ほどと同じグラフを計算したものが次のグラフです。確かに両対数グラフ上で等間隔になっていることが確認できます。

一応データの生成用のコードとgnuplot用のスクリプトを示します。gnuplotでは$\theta_c$を定義しておくことで、簡単に$\theta-\theta_c$を横軸にすることができます。

data_test.py

import numpy as np
theta_c = 2.0
theta_start = theta_c + 0.001
theta_end = 5.0

def f(x):
	return 2.0 * np.sqrt(x-theta_c)

N = 50

#linspace
theta_list = np.linspace(theta_start,theta_end,N)

data = [[] for i in range(N)]
for i in range(N):
	data[i] = [theta_list[i], f(theta_list[i])]

data = np.array(data)
np.savetxt("data_linspace.txt",data)

#power
r = ( (theta_end-theta_c)/(theta_start-theta_c) )**(1.0/(N-1))
theta_list = [theta_c + r**float(i)*(theta_start-theta_c) for i in range(N)]

data = [[] for i in range(N)]
for i in range(N):
	data[i] = [theta_list[i],f(theta_list[i])]
data = np.array(data)
np.savetxt("data_power.txt",data)

test_plt.txt

theta_c = 2.0
set xlabel "{/symbol q}-{/symbol q}_c"
set ylabel "f({/symbol q})"
set logscale
#plot "data_linspace.txt" u (($1)-theta_c):2 title "data"
plot "data_power.txt" u (($1)-theta_c):2 title "data"
set terminal jpeg
#set output "data_linspace.jpg"
set output "data_power.jpg"
replot
q