物理学者(ポスドク)による日々の研究生活のメモ書きです ( python/emacs/html/Japascript/シェルスクリプト/TeX/Mac/C言語/Linux/git/tmux/R/ポケモンGO)
あるデータの1列目の累積分布を書きたいときは以下のような感じでOK
nはデータ点数、これが間違ってると累積分布が[0:1]の範囲に収まらない
n=100
plot "exp.dat" using 1:(1/n) smooth cumulative
plot "exp.dat" using 1:(1/n) smooth cumulative
他に必要であろう設定を加えたのが以下の通り
n=100
set size square
set key bottom right
set xlabel "p-value"
set ylabel "cumulative distribution"
plot "p-value.txt" using 1:(1/n) smooth cumul title "simulation", x w l lw 2 title "uniform"
set size square
set key bottom right
set xlabel "p-value"
set ylabel "cumulative distribution"
plot "p-value.txt" using 1:(1/n) smooth cumul title "simulation", x w l lw 2 title "uniform"
PR
昔は、linetypeでなんかやった覚えがあるようなないような・・・
pdfでは点線にならなくて、そのときtestというコマンドを覚えた気がします
調べたら8年前でした(mjd)↓
■ 過去記事 : 【Gnuplot】plotするときの線の太さ、色、点線など
gnuplot5では点線のやり方が変わったらしい
■ 参考 : 破線の使用(dashtype)
その後の(10, 5)は白い空白と実線の間隔を決めるコマンド
これを変えれば、いろんなパターンの点線を書ける
めんどくさかったら dt 1とかでOK
注意点としては、事前に今使っているターミナルで点線が使えるか確認しておくと良い
↓ set terminal pngでは点線は使えない
↓ set terminal pngcairoでは点線が使える
ツイート
pdfでは点線にならなくて、そのときtestというコマンドを覚えた気がします
調べたら8年前でした(mjd)↓
■ 過去記事 : 【Gnuplot】plotするときの線の太さ、色、点線など
gnuplot5では点線のやり方が変わったらしい
■ 参考 : 破線の使用(dashtype)
plot sin(x) with lines dt (10,5)
という感じでdtを使えばOKその後の(10, 5)は白い空白と実線の間隔を決めるコマンド
これを変えれば、いろんなパターンの点線を書ける
めんどくさかったら dt 1とかでOK
plot sin(x) with lines dt (10, 10, 3, 10)
とか注意点としては、事前に今使っているターミナルで点線が使えるか確認しておくと良い
↓ set terminal pngでは点線は使えない
↓ set terminal pngcairoでは点線が使える
だいぶ前に使って、ブックマークに入れたまんまだったので代わりにここに貼り付けとく
■ 参考 : gnuplotでファイルの一行目をtitleにする
ファイルの1行目になんか文字があって
file = "hoge.txt"
time = system("head -1 " . file)
set title sprintf("%2.2f hrs", time / 3600)
time = system("head -1 " . file)
set title sprintf("%2.2f hrs", time / 3600)
ファイルの1行目を#から始まるようにしておかないとこのファイルをplotとかしたときに抵触する気がする
または文字が書いてあるだけのファイルを別に用意しておけばいい?
もう一つ別の例
■ 参考 : gnuplot
これはplotのtitleにファイルから読み込んだ文字を入力したい場合
図の凡例にデータファイルから読み込んだデータを入れたいときには, title column (line number)とすればいい。(参考URL)
p "filename" u ($0):2:xtic(1) title column(2) w p
p "filename" u ($0):2:xtic(1) title column(2) w p
これ↓
■ 参考 : Gnuplotting/gnuplot-palettes
ここに乗ってるサンプルの配色例↓
使い方は
load 'hoge.pal'
hoge.palを自分が使いたいものに置き換えればOK例えばbrbg.palの中身はこんな感じ
これを直接 ~/.gnuplot に書き込んでおいてもOK
# line styles for ColorBrewer BrBG
# for use with divering data
# provides 8 colors with brown low, white middle, and blue-green high
# compatible with gnuplot >=4.2
# author: Anna Schneider
# line styles
set style line 1 lt 1 lc rgb '#8C510A' # dark brown
set style line 2 lt 1 lc rgb '#BF812D' # medium brown
set style line 3 lt 1 lc rgb '#DFC27D' #
set style line 4 lt 1 lc rgb '#F6E8C3' # pale brown
set style line 5 lt 1 lc rgb '#C7EAE5' # pale blue-green
set style line 6 lt 1 lc rgb '#80CDC1' #
set style line 7 lt 1 lc rgb '#35978F' # medium blue-green
set style line 8 lt 1 lc rgb '#01665E' # dark blue-green
# palette
set palette defined ( 0 '#8C510A',\
1 '#BF812D',\
2 '#DFC27D',\
3 '#F6E8C3',\
4 '#C7EAE5',\
5 '#80CDC1',\
6 '#35978F',\
7 '#01665E' )
# for use with divering data
# provides 8 colors with brown low, white middle, and blue-green high
# compatible with gnuplot >=4.2
# author: Anna Schneider
# line styles
set style line 1 lt 1 lc rgb '#8C510A' # dark brown
set style line 2 lt 1 lc rgb '#BF812D' # medium brown
set style line 3 lt 1 lc rgb '#DFC27D' #
set style line 4 lt 1 lc rgb '#F6E8C3' # pale brown
set style line 5 lt 1 lc rgb '#C7EAE5' # pale blue-green
set style line 6 lt 1 lc rgb '#80CDC1' #
set style line 7 lt 1 lc rgb '#35978F' # medium blue-green
set style line 8 lt 1 lc rgb '#01665E' # dark blue-green
# palette
set palette defined ( 0 '#8C510A',\
1 '#BF812D',\
2 '#DFC27D',\
3 '#F6E8C3',\
4 '#C7EAE5',\
5 '#80CDC1',\
6 '#35978F',\
7 '#01665E' )
high sierraにアップデートしたのでgnuplotのversionもv5.2になりました
デフォルトのterminalがaquaだったのでそのまましばらく使ってました
なんかグラフを書いたときに、aqua termが初めて立ち上がったときはその画面がアクティブになるものの、
replotしたときは画面が一番うしろにいったっきりで、最前面に出てこない・・・
なんとかできないものかなぁ〜とあれこれ調べたけど、解決策は見つからず
しょうがないので、XQuartzでx11 terminalを使ってます・・・
aquaのが描画はきれいだから、そっちを使いたいんだけどなぁ〜
いい解決方法を見つけたら追記します
数ヶ月前におっぱい関数というのがあることを初めて知った・・・
■ 参考記事 : gnuplot でおっぱい
その記事は以下のツイートに触発されて書かれてるっぽい
元のツイートではMathematicaで書かれてるけど、それをgnuplotに置き換えたものらしい
さらに遡ってみると、2011年にtwitter上でおっぱいの方程式についてツイートを発見した
さらに最近見つけたのはこの動画
これらのおっぱい関数について私の感想
ツイート
■ 参考記事 : gnuplot でおっぱい
その記事は以下のツイートに触発されて書かれてるっぽい
元のツイートではMathematicaで書かれてるけど、それをgnuplotに置き換えたものらしい
某工大の叡智の結集 pic.twitter.com/TBhRfTsNrn
— T s u z u k i (@ttzux19) 2014年11月4日
さらに遡ってみると、2011年にtwitter上でおっぱいの方程式についてツイートを発見した
いつの時代も天才はいるものだ・・・おっぱいの方程式を発見したよ! pic.twitter.com/nM2rtxXC
— しもまっしぐま (@shimoMathSiGMA) 2011年11月6日
さらに最近見つけたのはこの動画
これらのおっぱい関数について私の感想
ちょい前からgnuplotのメーリングリストに参加している
開発の最新情報や「これわからんのやけどどーしたらいいんや?」という疑問などが1日に1メール以上の割合でやり取りされているっぽい
質問に関しては、質問する側の説明不足がやや見受けられる気もするが、答える側も辛抱強く答えてくれてるのを端から見てて感じた
その中で二重ループ(double loop)のやり方について質問してた人がいた
そのうち使うかもしれないので、コードだけメモ・・・
title(i, j)=sprintf("i=%d, j=%d", i, j)
set xrange [0:pi]
plot for [i=1:3] for [j=1:3] sin(i*j*x) title title(i, j)
set xrange [0:pi]
plot for [i=1:3] for [j=1:3] sin(i*j*x) title title(i, j)
今は必要じゃないので調べないけど気になった疑問
・for [i=1:3]でforループの開始になるのはいいとして、終了場所はどこなのか?
・もし複数行に渡ったループを書きたいときはどうしたらいいのか? {}とかでくくればいける?
gnuplotでエラーバーをつける方法は以前まとめた
■ 過去記事 : 【gnuplot5】エラーバーを付けたい
今回はエラーバーじゃなくて、エラーの領域を少し違う色で薄く塗りつぶしたグラフを書きたい
検索したらあった
■ 参考 : gnuplot demo script: errorbars.dem
これの一番下の例
ただ、v5.2なのでもしかしたらけっこう最近追加された機能の可能性があるなぁ・・・
実際v5.0のdemoにはなかったので・・・(記事タイトルを修正)
■ 参考 : Gnuplot smooth confidence interval lines as opposed to error bars
こっちの例は簡単で with filledcurvesを使って
エラー領域の下限値〜上限値の範囲を塗りつぶすプロット + センター値のプロットで2回プロットしてる
やるならこっちか
ツイート
■ 過去記事 : 【gnuplot5】エラーバーを付けたい
今回はエラーバーじゃなくて、エラーの領域を少し違う色で薄く塗りつぶしたグラフを書きたい
検索したらあった
■ 参考 : gnuplot demo script: errorbars.dem
これの一番下の例
ただ、v5.2なのでもしかしたらけっこう最近追加された機能の可能性があるなぁ・・・
実際v5.0のdemoにはなかったので・・・(記事タイトルを修正)
■ 参考 : Gnuplot smooth confidence interval lines as opposed to error bars
こっちの例は簡単で with filledcurvesを使って
エラー領域の下限値〜上限値の範囲を塗りつぶすプロット + センター値のプロットで2回プロットしてる
やるならこっちか
twitter で見かけたこれをやってみたい
■ 参考 : gnuplotでハートの形を描く
陰関数を使う方法は、
■ 参考 : 3.7 どうやったら陰関数のグラフが書けますか
を参考に
注意としてはset tableしたあとでunset tableしないとその後のコマンドがきちんと動かない
で、できたのがこれ
コマンドは
with filledcurve
でやると、塗りつぶしになっていい感じの図になるかと思ったけど、
できたのは下の図
その値〜x軸までの範囲を塗りつぶしてるような・・・
ツイート
■ 参考 : gnuplotでハートの形を描く
陰関数を使う方法は、
■ 参考 : 3.7 どうやったら陰関数のグラフが書けますか
を参考に
注意としてはset tableしたあとでunset tableしないとその後のコマンドがきちんと動かない
で、できたのがこれ
コマンドは
len=1.5
set xrange [-len:len]
set yrange [-len:len]
set zrange [-len:len]
f(x,y)= (x**2+y**2-1)**3-x**2*y**3
set contour base
set cntrparam levels discrete 0.0
unset surface
set isosample 400, 400
set table 'curve.txt'
splot f(x, y)
unset table
set term png
set output "hoge.png"
#plot [-len:len][-1.0:len] "curve.txt" with filledcurve lc rgb 'pink' title "heart shape"
plot [-len:len][-1.0:len] "curve.txt" with line lw 4 lc rgb 'pink' title "heart shape"
最後の部分をset xrange [-len:len]
set yrange [-len:len]
set zrange [-len:len]
f(x,y)= (x**2+y**2-1)**3-x**2*y**3
set contour base
set cntrparam levels discrete 0.0
unset surface
set isosample 400, 400
set table 'curve.txt'
splot f(x, y)
unset table
set term png
set output "hoge.png"
#plot [-len:len][-1.0:len] "curve.txt" with filledcurve lc rgb 'pink' title "heart shape"
plot [-len:len][-1.0:len] "curve.txt" with line lw 4 lc rgb 'pink' title "heart shape"
with filledcurve
でやると、塗りつぶしになっていい感じの図になるかと思ったけど、
できたのは下の図
その値〜x軸までの範囲を塗りつぶしてるような・・・
f(x)=1/x
set samples 11
plot [0:10][0:1] f(x) w lp lw 2 pt 6 ps 3
set term png
set output "hoge.png"
rep
set samples 11
plot [0:10][0:1] f(x) w lp lw 2 pt 6 ps 3
set term png
set output "hoge.png"
rep
最初
set samples 10
にしてて、なんか数字が合わないなぁ〜・・・と悩んでたけど、1~10を10区切るならまだしも
0~10を10個に区切ったら半端になるわなぁ〜・・・
■ 参考 : 87.34 凡例 (key)
set key left top
とかではtitleの表示する場所が左上になってしまうtitleの文字を左揃えにしたいときは
set key Left
また、titleとサンプル直線の場所を逆にしたいときはreverseを追加する
set key Left reverse
gnuplotでfittingをする方法は以前まとめた
■ 過去記事 : 【gnuplot】を使ってヒストグラムをRayleigh分布でfittingしてみる
このfittingの出力結果を取り出してさらに別のところで使いたい・・・・
fittingした結果を毎回、手でメモって代入する方法では5000個とかfittingするときに骨が折れるので・・・・
■ 参考 : 変数や関数のファイルへの書き出し(save)
これを参考にしてgnuplotで、
あとはシェルスクリプトの中でこの値を取り出したいときは、
参考までに、
もう一度gnuplotを起動して
先ほどの変数たち a, b, cとかfittingのエラーとかがそれぞれの変数に代入されてるはず
ツイート
■ 過去記事 : 【gnuplot】を使ってヒストグラムをRayleigh分布でfittingしてみる
このfittingの出力結果を取り出してさらに別のところで使いたい・・・・
fittingした結果を毎回、手でメモって代入する方法では5000個とかfittingするときに骨が折れるので・・・・
■ 参考 : 変数や関数のファイルへの書き出し(save)
これを参考にしてgnuplotで、
f(x)=a*exp(-(x-b)*(x-b)/2.0/c/c); a=10000; b=0.001; c=1
fit f(x) "hist.txt" u 1:3 via a, b, c
plot [:][0.001:] f(x), "hist.txt" u 1:3
save variables "test.var"
とすると、以下のようなファイルが生成されるはずfit f(x) "hist.txt" u 1:3 via a, b, c
plot [:][0.001:] f(x), "hist.txt" u 1:3
save variables "test.var"
#
#
# G N U P L O T
# Version 5.0 patchlevel 3 last modified 2016-02-21
#
# Copyright (C) 1986-1993, 1998, 2004, 2007-2016
# Thomas Williams, Colin Kelley and many others
#
# gnuplot home: http://www.gnuplot.info
# faq, bugs, etc: type "help FAQ"
# immediate help: type "help" (plot window: hit 'h')
GNUTERM = "qt"
a = 23470.5784871032
b = -0.788993648242623
c = 17.8279442010367
GPFUN_f = "f(x)=a*exp(-(x-b)*(x-b)/2.0/c/c)"
x = 0.0
FIT_CONVERGED = 1
FIT_NDF = 97
FIT_STDFIT = 74.9163111531704
FIT_WSSR = 544408.006649469
FIT_P = 0.0
FIT_NITER = 6
a_err = 20.6412240133938
b_err = 0.0181088448166268
c_err = 0.0181088509499784
# EOF
#
# G N U P L O T
# Version 5.0 patchlevel 3 last modified 2016-02-21
#
# Copyright (C) 1986-1993, 1998, 2004, 2007-2016
# Thomas Williams, Colin Kelley and many others
#
# gnuplot home: http://www.gnuplot.info
# faq, bugs, etc: type "help FAQ"
# immediate help: type "help" (plot window: hit 'h')
GNUTERM = "qt"
a = 23470.5784871032
b = -0.788993648242623
c = 17.8279442010367
GPFUN_f = "f(x)=a*exp(-(x-b)*(x-b)/2.0/c/c)"
x = 0.0
FIT_CONVERGED = 1
FIT_NDF = 97
FIT_STDFIT = 74.9163111531704
FIT_WSSR = 544408.006649469
FIT_P = 0.0
FIT_NITER = 6
a_err = 20.6412240133938
b_err = 0.0181088448166268
c_err = 0.0181088509499784
# EOF
あとはシェルスクリプトの中でこの値を取り出したいときは、
var_a=`grep "^a =" test.var | awk '{print $3}'`
var_b=`grep "^b =" test.var | awk '{print $3}'`
var_c=`grep "^c =" test.var | awk '{print $3}'`
みたいな感じで取り出せるはずvar_b=`grep "^b =" test.var | awk '{print $3}'`
var_c=`grep "^c =" test.var | awk '{print $3}'`
参考までに、
もう一度gnuplotを起動して
load "test.var"
とすると、先ほどの変数たち a, b, cとかfittingのエラーとかがそれぞれの変数に代入されてるはず
それぞれgamma(x)とigamma(a, x)とかで呼び出せる
ヘルプは以下のような感じ
gnuplot> help gamma
The `gamma(x)` function returns the gamma function of the real part of its
argument. For integer n, gamma(n+1) = n!. If the argument is a complex
value, the imaginary component is ignored.
gnuplot> help igamma
The `igamma(a,x)` function returns the normalized incomplete gamma
function of the real parts of its arguments, where a > 0 and x >= 0.
The standard notation is P(a,x), e.g. Abramowitz and Stegun (6.5.1),
with limiting value of 1 as x approaches infinity. If the arguments
are complex, the imaginary components are ignored.
The `gamma(x)` function returns the gamma function of the real part of its
argument. For integer n, gamma(n+1) = n!. If the argument is a complex
value, the imaginary component is ignored.
gnuplot> help igamma
The `igamma(a,x)` function returns the normalized incomplete gamma
function of the real parts of its arguments, where a > 0 and x >= 0.
The standard notation is P(a,x), e.g. Abramowitz and Stegun (6.5.1),
with limiting value of 1 as x approaches infinity. If the arguments
are complex, the imaginary components are ignored.
■ 参考 : x軸が共通のplotを複数並べる場合(stacked plots by multiplot w/ set origin, set size)
↑
ここに何例か、実際のコード付きで載っている
■ 参考 : Gnuplot で複数のグラフを並べる方法
↑
こっちは横に2枚のグラフを並べるとき
■ 参考 : 2次元プロットのあれこれ (その6)
2つめのページでは、2枚を並べるときの余白の調整をmarginでやるとずれるのでやめた方がいいって書いてある
ただ、これについてはmarginがないとlabelとか軸に付いてる数字がグラフの範囲外にはみ出ることがあるので
やはり微調整をset bmarginとかset tmarginでやって方がいいんじゃないかなァ〜と自分は思う・・・
そのうち、実例とともに記事にしたいけど、今回はブックマーク整理のためのメモ記事で
ホント感動した・・・
こんな図がgnuplotでも作れる
(正確にはシェルコマンドで綺麗に整形したデータをプロットしてるだけだけど・・・そこは野暮か・・・)
■ 参考 : Plotting raster data from Natural Earth
プロットするためのコマンドは、
■ 参考 : プロットコマンド
プロットに使うデータは、
■ 参考 : プロットに必要なデータ
このデータは convert_natural_earth.sh という(おそらく)自作のスクリプトで作られている
使い方も書いてあるけど、元の世界地図のTIFデータ(?)がないと使えない・・・・
■ 参考 : convert_natural_earth
上記のデータよりももうちょい解像度が高いデータはgithubにあった
■ 参考 : もうちょい解像度が高いデータ
これで書いたのが以下、
関係ありそうな記事リンク
■ 参考 : Plotting the world
■ 過去記事 : 【gnuplot】塗りつぶしプロット
■ 過去記事 : 【gnuplot】world.datを活用したい
(2017/05/22 追記)
TIFファイルという地図データは
■ 過去記事 : 1:50m Raster Data - Natural Earth
からダウンロードできるっぽい
上にリンクを貼ったデータはそれぞれ
・Cross-blended Hypsometric Tints
・Natural Earth 2
からダウンロードされたデータを加工したもの
加工するためのシェルスクリプト convert_natural_earth.sh をダウンロードしてきたTIFファイルにかけても、エラーが出ておかしなデータが出てきてしまう・・・
原因はMacのsedとLINUXのsedが違う挙動をすることで、
これはMacのsedで-iとするときはバックアップのファイルのための拡張子を指定しないといけないらしい
なので以下ではsed -i '.bak'と置き換えた
あとはpm3dのために、適時空行をいれないといけないが、元々のソースコードにあったsedではうまく動かなかったのでawkで代用
途中で257で割っているのは、RGBの値がきちんと取得できていないので、よくわからんが257で割っている
時間をかけたらうまくできそうだけど、めんどくさいからいいや・・・・
n=2000 でプロットしたのが以下の図
この図の上に色々なものを書いていきたい
ツイート
こんな図がgnuplotでも作れる
(正確にはシェルコマンドで綺麗に整形したデータをプロットしてるだけだけど・・・そこは野暮か・・・)
■ 参考 : Plotting raster data from Natural Earth
プロットするためのコマンドは、
■ 参考 : プロットコマンド
プロットに使うデータは、
■ 参考 : プロットに必要なデータ
このデータは convert_natural_earth.sh という(おそらく)自作のスクリプトで作られている
使い方も書いてあるけど、元の世界地図のTIFデータ(?)がないと使えない・・・・
■ 参考 : convert_natural_earth
上記のデータよりももうちょい解像度が高いデータはgithubにあった
■ 参考 : もうちょい解像度が高いデータ
これで書いたのが以下、
関係ありそうな記事リンク
■ 参考 : Plotting the world
■ 過去記事 : 【gnuplot】塗りつぶしプロット
■ 過去記事 : 【gnuplot】world.datを活用したい
(2017/05/22 追記)
TIFファイルという地図データは
■ 過去記事 : 1:50m Raster Data - Natural Earth
からダウンロードできるっぽい
上にリンクを貼ったデータはそれぞれ
・Cross-blended Hypsometric Tints
・Natural Earth 2
からダウンロードされたデータを加工したもの
加工するためのシェルスクリプト convert_natural_earth.sh をダウンロードしてきたTIFファイルにかけても、エラーが出ておかしなデータが出てきてしまう・・・
原因はMacのsedとLINUXのsedが違う挙動をすることで、
extra characters at the end of N command
みたいなエラーが出てくるが、これはMacのsedで-iとするときはバックアップのファイルのための拡張子を指定しないといけないらしい
なので以下ではsed -i '.bak'と置き換えた
あとはpm3dのために、適時空行をいれないといけないが、元々のソースコードにあったsedではうまく動かなかったのでawkで代用
途中で257で割っているのは、RGBの値がきちんと取得できていないので、よくわからんが257で割っている
時間をかけたらうまくできそうだけど、めんどくさいからいいや・・・・
#!/bin/sh
# data : http://www.naturalearthdata.com/downloads/50m-raster-data/
if [ $# -ne 2 ]; then
echo "usage : `basename $0` [resolution] [inputfile]"
exit
fi
RES="$1"
IN="$2"
SCALED="${IN%.tif}_${RES}px.tif"
OUT="${SCALED%.tif}.txt"
# resize tiff image to desired resolution
convert "$IN" -resize $RES "$SCALED"
# convert to txt file
convert "$SCALED" "$OUT"
# remove unneeded header
sed -i '.bak' '1d' "$OUT"
# replace ": (" with ","
sed -i '.bak' 's/: (/\,/' "$OUT"
# remove all unneeded entries at end of line
sed -i '.bak' 's/).*$//' "$OUT"
# convert longitude and latitude to the right values
# 0..RES => -180..180
# 0..RES/2 => -90..89
awk -F, -v res=${RES} '{$1=$1*360/res-180; $2=-1*($2*360/res-90); }1' OFS=, "$OUT" > tmp.txt
awk -F, '{print $1, $2, $3/257, $4/257, $5/257}' OFS=, tmp.txt > tmp2.txt
mv tmp2.txt ${OUT}
rm tmp.txt
# add blank lines for pm3d
awk -v n=${RES} '{print $0} NR>2 && 0==NR % n{printf "\n"}' ${OUT} > tmp.txt
mv tmp.txt ${OUT}
# add new header
awk 'BEGIN{print "# long, lat, r, g, b"} {print $0}' ${OUT} > tmp.txt
mv tmp.txt ${OUT}
rm *.bak
これならMacでも動くはず・・・・# data : http://www.naturalearthdata.com/downloads/50m-raster-data/
if [ $# -ne 2 ]; then
echo "usage : `basename $0` [resolution] [inputfile]"
exit
fi
RES="$1"
IN="$2"
SCALED="${IN%.tif}_${RES}px.tif"
OUT="${SCALED%.tif}.txt"
# resize tiff image to desired resolution
convert "$IN" -resize $RES "$SCALED"
# convert to txt file
convert "$SCALED" "$OUT"
# remove unneeded header
sed -i '.bak' '1d' "$OUT"
# replace ": (" with ","
sed -i '.bak' 's/: (/\,/' "$OUT"
# remove all unneeded entries at end of line
sed -i '.bak' 's/).*$//' "$OUT"
# convert longitude and latitude to the right values
# 0..RES => -180..180
# 0..RES/2 => -90..89
awk -F, -v res=${RES} '{$1=$1*360/res-180; $2=-1*($2*360/res-90); }1' OFS=, "$OUT" > tmp.txt
awk -F, '{print $1, $2, $3/257, $4/257, $5/257}' OFS=, tmp.txt > tmp2.txt
mv tmp2.txt ${OUT}
rm tmp.txt
# add blank lines for pm3d
awk -v n=${RES} '{print $0} NR>2 && 0==NR % n{printf "\n"}' ${OUT} > tmp.txt
mv tmp.txt ${OUT}
# add new header
awk 'BEGIN{print "# long, lat, r, g, b"} {print $0}' ${OUT} > tmp.txt
mv tmp.txt ${OUT}
rm *.bak
n=2000 でプロットしたのが以下の図
この図の上に色々なものを書いていきたい
ただのメモ
■ 参考 : Bending the arrows - "delaying" the plot
普通にset arrow 1とかで矢印を書くと、直線の矢印しか書けない
ので、ここでは定義したf(x)って関数の一部だけをプロットしたものに、擬似的に矢印の先っぽをくっつけて曲がった矢印を再現してるんだと思う・・・(たぶん)
うちの環境ではリンク先と同じプロットはできなかった・・・
なぜか矢印の矢じりの部分が表示されない・・・・
ツイート
■ 参考 : Bending the arrows - "delaying" the plot
eps = 0.001
set style arrow 1 head filled size screen 0.03, 15, 45 lt -1
cut(x,x1,x3) = ((x >= x1 && x <= x1 + (1.0-eps)*(x3-x1)) ? 1.0 : 1/0)
f(x) = 0.5+(x-1)*(x-1.2)*(x-1.4)
x1 = 0.5
x3 = 1.95
new_x = x1 + (1.0-eps)*(x3-x1)
set arrow from new_x, f(new_x) to x3,f(x3) as 1
plot [0:3] sin(x) with point ps 1 pt 6, f(x)*cut(x,x1,new_x) with line lt -1
set style arrow 1 head filled size screen 0.03, 15, 45 lt -1
cut(x,x1,x3) = ((x >= x1 && x <= x1 + (1.0-eps)*(x3-x1)) ? 1.0 : 1/0)
f(x) = 0.5+(x-1)*(x-1.2)*(x-1.4)
x1 = 0.5
x3 = 1.95
new_x = x1 + (1.0-eps)*(x3-x1)
set arrow from new_x, f(new_x) to x3,f(x3) as 1
plot [0:3] sin(x) with point ps 1 pt 6, f(x)*cut(x,x1,new_x) with line lt -1
普通にset arrow 1とかで矢印を書くと、直線の矢印しか書けない
ので、ここでは定義したf(x)って関数の一部だけをプロットしたものに、擬似的に矢印の先っぽをくっつけて曲がった矢印を再現してるんだと思う・・・(たぶん)
うちの環境ではリンク先と同じプロットはできなかった・・・
なぜか矢印の矢じりの部分が表示されない・・・・
前回は塗りつぶしのプロットを書いた
■ 過去記事 : 【gnuplot 6】gsl-histogramの結果を塗りつぶしてプロットしたい
(今気づいたけどgnuplot6になってますよ・・・おいおい・・・)
今回は、with stepsでこれをプロットしたい
今まで単一のデータをプロットするときは、塗りつぶしプロットでも用は足りていたが、2種類のデータを塗りつぶしのwith boxesで重ねると
多少は透過させているとはいえ、見にくい・・・・
単純にwith stepsで2つのデータをプロットすると、一見綺麗に書けるが
set log yとしてみると、データ点が入っていないbin(y軸の値が0になるbin)の処理がlogスケールなので、想像と違うことになってしまう・・・
↑
カウントが0のbinの端の処理が嫌・・・・
with boxesでもいいかなぁ〜と思っていたが解決方法に気づいたので、ここに書いとく・・・
まずはデータの用意
gsl-histogramで0から50の間を30個のbinで分けてヒストグラムのデータを生成する
特に問題なさそう・・・・
けど、もっとデータ点が増えてくると、y軸をlogスケールにしたくなる・・・
矢印を書いたあたりの側面の処理がなんか嫌・・・
with stepsは(x_1, y_1)の点から(x_2, y_2)の点に向かって垂直な直線を引くようなプロットらしい
それでそのうちのいずれかがy=0だとy軸がログスケールだと、直線が引かれない仕様らしい
思いついた解決方法は、三項演算子でy=0になっている点はy=0.0001のようなグラフに影響が出ないような小さい数字を代わりに当てはめてプロットする方法
ツイート
■ 過去記事 : 【gnuplot 6】gsl-histogramの結果を塗りつぶしてプロットしたい
(今気づいたけどgnuplot6になってますよ・・・おいおい・・・)
今回は、with stepsでこれをプロットしたい
今まで単一のデータをプロットするときは、塗りつぶしプロットでも用は足りていたが、2種類のデータを塗りつぶしのwith boxesで重ねると
多少は透過させているとはいえ、見にくい・・・・
単純にwith stepsで2つのデータをプロットすると、一見綺麗に書けるが
set log yとしてみると、データ点が入っていないbin(y軸の値が0になるbin)の処理がlogスケールなので、想像と違うことになってしまう・・・
↑
カウントが0のbinの端の処理が嫌・・・・
with boxesでもいいかなぁ〜と思っていたが解決方法に気づいたので、ここに書いとく・・・
まずはデータの用意
gsl-randist 0 1000 rayleigh 10 | gsl-histogram 0 50 30 > hoge.txt
gsl-randistで1000個の乱数データを生成して、gsl-histogramで0から50の間を30個のbinで分けてヒストグラムのデータを生成する
plot [:][0.1:] "hoge.txt" u 1:3 w steps lw 3 lc 2
でこれ特に問題なさそう・・・・
けど、もっとデータ点が増えてくると、y軸をlogスケールにしたくなる・・・
set log y
plot [:][0.1:] "hoge.txt" u 1:3 w steps lw 3 lc 2
plot [:][0.1:] "hoge.txt" u 1:3 w steps lw 3 lc 2
矢印を書いたあたりの側面の処理がなんか嫌・・・
with stepsは(x_1, y_1)の点から(x_2, y_2)の点に向かって垂直な直線を引くようなプロットらしい
それでそのうちのいずれかがy=0だとy軸がログスケールだと、直線が引かれない仕様らしい
思いついた解決方法は、三項演算子でy=0になっている点はy=0.0001のようなグラフに影響が出ないような小さい数字を代わりに当てはめてプロットする方法
set log y
plot [:][0.1:] "hoge.txt" u 1:($3 == 0 ? 0.001 : $3) w steps lw 3 lc 2
plot [:][0.1:] "hoge.txt" u 1:($3 == 0 ? 0.001 : $3) w steps lw 3 lc 2
plot "hoge.txt" every n
で最初の1点を打った後、次の(n-1)点を飛ばして、その次はn番目のデータをプロットする実際に使うときはnじゃなくて(n=3と事前に打っておくならnでもいいけど)3とか整数値に置き換えて使う
plot "hoge.txt" every n u 1:3 w lp
のようにプロットタイプはeveryの後に打つっぽいevery nの前に打つと、エラーが出たのでダメっぽい
このeveryはfor文のようにもっと複雑なことができるけど、そこまでは・・・・いいや・・・
昔書いた気もするけど、グラフにエラーバーを付けたことがあるのって4年生以来だから一体何年舞なんだよ・・・って感じだ・・・・
とはいえ、実際に測定したものにエラーバーを付けるのではなくて
自分の場合ほとんどがモンテカルロシミュレーションを繰り返しやって、それにエラーバーを付ける程度の使い方だけど・・・・
どういうエラーバーを付けられるかは、次のページを見てもらうのが一番早い
■ 参考 : gnuplot demo script: mgr.dem
X軸方向のエラーバー、Y軸方向のエラーバー、X軸とY軸両方に伸びてるエラーバー全部できる
ただし、ちょっとめんどくさい制約が・・・・
詳細は次のページに書かれてる
■ 参考 : 34.2 Errorbars
まずはデータ作り
u 1:2:3だけではダメで、それにw e(with errorbarsの略)を付ければOK
データ点を線で結びながら、エラーバーを付ける場合は2回プロットしないといけなくて、
ここまで調べてわかったのが w eはwith yerrorbarsの略称らしい
他には
w xe : x軸方向のエラーバー
w ye : y軸方向のエラーバー
w xye : x軸とy軸方向両方にエラーバー
があるらしい
w xyeを使う場合はu 1:2:3:4のように4つ使うデータ列を指定しないといけない
例えば、
他の形のエラーバーを使いたかったら、以下を参照
■ 参考 : gnuplot demo script: mgr.dem
ツイート
とはいえ、実際に測定したものにエラーバーを付けるのではなくて
自分の場合ほとんどがモンテカルロシミュレーションを繰り返しやって、それにエラーバーを付ける程度の使い方だけど・・・・
どういうエラーバーを付けられるかは、次のページを見てもらうのが一番早い
■ 参考 : gnuplot demo script: mgr.dem
X軸方向のエラーバー、Y軸方向のエラーバー、X軸とY軸両方に伸びてるエラーバー全部できる
ただし、ちょっとめんどくさい制約が・・・・
詳細は次のページに書かれてる
■ 参考 : 34.2 Errorbars
まずはデータ作り
% seq 1 20 > hoge.txt
% awk '{x=rand*10; print $1, x, sqrt(x)}' hoge.txt > foo.txt
% cat foo.txt
1 8.40188 2.8986
2 3.94383 1.98591
3 7.83099 2.79839
4 7.9844 2.82567
5 9.11647 3.01935
6 1.97551 1.40553
7 3.35223 1.83091
8 7.6823 2.7717
9 2.77775 1.66666
10 5.5397 2.35366
11 4.77397 2.18494
12 6.28871 2.50773
13 3.64784 1.90993
14 5.13401 2.26584
15 9.5223 3.08582
16 9.16195 3.02687
17 6.35712 2.52133
18 7.17297 2.67824
19 1.41603 1.18997
20 6.06969 2.46367
% awk '{x=rand*10; print $1, x, sqrt(x)}' hoge.txt > foo.txt
% cat foo.txt
1 8.40188 2.8986
2 3.94383 1.98591
3 7.83099 2.79839
4 7.9844 2.82567
5 9.11647 3.01935
6 1.97551 1.40553
7 3.35223 1.83091
8 7.6823 2.7717
9 2.77775 1.66666
10 5.5397 2.35366
11 4.77397 2.18494
12 6.28871 2.50773
13 3.64784 1.90993
14 5.13401 2.26584
15 9.5223 3.08582
16 9.16195 3.02687
17 6.35712 2.52133
18 7.17297 2.67824
19 1.41603 1.18997
20 6.06969 2.46367
plot "foo.txt" u 1:2:3 w e
でエラーバー付きのグラフが書けるu 1:2:3だけではダメで、それにw e(with errorbarsの略)を付ければOK
plot "foo.txt" u 1:2:(sqrt($2)) w e
でも同じプロットが書けるはずデータ点を線で結びながら、エラーバーを付ける場合は2回プロットしないといけなくて、
plot "foo.txt" u 1:2:3 w e, "" u 1:2 w l lw 2
ここまで調べてわかったのが w eはwith yerrorbarsの略称らしい
他には
w xe : x軸方向のエラーバー
w ye : y軸方向のエラーバー
w xye : x軸とy軸方向両方にエラーバー
があるらしい
w xyeを使う場合はu 1:2:3:4のように4つ使うデータ列を指定しないといけない
例えば、
plot "foo.txt" u 1:2:($2/10.0):($2/10.0) w xye lc 7
で、他の形のエラーバーを使いたかったら、以下を参照
■ 参考 : gnuplot demo script: mgr.dem
こういうcontourにlabelを書き込むことがやりたかった
一応やり方は記事にまとめられてるけど、よくわからん
■ 参考 : Contour plots
それで他の方法を探してみたら、こんなのを見つけた
■ 参考 : Labels with white background in LaTeX terminals
■ 参考 :Label size in epslatex terminal
■ 参考 :Images within a graph
特に最後の写真にlabelを貼付けるのは職人すぎるわ・・・・
自分ならkeynoteで矢印を付けて、それを保存してるところだ・・・
■ 過去記事 :【gnuplot】を使ってリサージュ曲線のgifを作りたいとき、
以前こんな記事を書いた↑
プロットするデータを(例えば)100行ごとに空行を入れて生成すると
gnuplotでindexというオプションを使って分割プロットができる
(1つの空行でいけるはずだけど、なぜかうまくいかなかったので今回は2行入れた)
↑
(2017/06/07 追記)このトラブルに何度かハマったので赤太字にしといた
なんでだろう?
for文を使って1行で一括プロットするには
最初はi=0の場合とi!=0の場合をわざわざ分けて、replotとかを使って書いていたけど、そんなめんどくさいことはしなくてもいいっぽい
ちなみにこのfor loopの中でindex 2のデータにきちんと"2"とタイトルを付けたいときは
この書き方でそれぞれのタイトルがNo.1、No.2、、、、となるはず
■ 参考 : Loop structure inside gnuplot?
他にも調べたらeveryとかdo forとかでてきたけど、難し過ぎてよくわからん・・・・
一応リンクだけメモしておく
■ 参考 : gnuplotのコマンドまとめ(ループとかeveryとか)
■ 参考 : 80.10 Plot コマンドの for ループ (for loops in plot command)
■ 参考 : X行毎にデータをプロットする
ツイート
以前こんな記事を書いた↑
プロットするデータを(例えば)100行ごとに空行を入れて生成すると
gnuplotでindexというオプションを使って分割プロットができる
(1つの空行でいけるはずだけど、なぜかうまくいかなかったので今回は2行入れた)
↑
(2017/06/07 追記)このトラブルに何度かハマったので赤太字にしといた
なんでだろう?
for文を使って1行で一括プロットするには
plot for [i=0:40] "hoge.txt" index i
とかでOK最初はi=0の場合とi!=0の場合をわざわざ分けて、replotとかを使って書いていたけど、そんなめんどくさいことはしなくてもいいっぽい
ちなみにこのfor loopの中でindex 2のデータにきちんと"2"とタイトルを付けたいときは
plot for [i=0:40] "hoge.txt" index i title "No.".i
とするこの書き方でそれぞれのタイトルがNo.1、No.2、、、、となるはず
■ 参考 : Loop structure inside gnuplot?
他にも調べたらeveryとかdo forとかでてきたけど、難し過ぎてよくわからん・・・・
一応リンクだけメモしておく
■ 参考 : gnuplotのコマンドまとめ(ループとかeveryとか)
■ 参考 : 80.10 Plot コマンドの for ループ (for loops in plot command)
■ 参考 : X行毎にデータをプロットする
twitterで以下の記事に関するツイートを見かけた
■ 参考 : 【xkcd】pythonコードにたった一行で漫画のようなグラフを作る!
そもそも「xkcd」って漫画を知らない・・・
リンク先にあるように、
へぇ〜
で、もしかしてこれgnuplotでもできるんじゃない?
と思って「gnuplot xkcd」とかで調べてみたらあったわ・・・
■ 参考 : xkcd-gnuplot
このコードのミソはここで宣言されているjigglea関数jiggle関数
これによってsin(x)とcos(x)に揺らぎを与えて漫画調のぐらぐらした線を再現している
また漫画調のフォントは"Humor Sans"というらしい
ググったら以下のページに落ちていた
■ 参考 : xkcd-desktop/Humor-Sans.ttf
ダウンロードしたttfファイルをダブルクリックすれば、Macなら簡単にインストールできるはず
インストールしたフォントファイルにPATHを通さないとgnuplotでは使えないので注意
PATHの通し方は以下を参考に
■ 過去記事 : 【gnuplot】Could not find/open font...というエラーが出るのでフォントの指定をする
これでこんな感じのグラフが書ける
コードは以下の通り、
(オリジナルからpng出力に変更しただけ)
ツイート
■ 参考 : 【xkcd】pythonコードにたった一行で漫画のようなグラフを作る!
そもそも「xkcd」って漫画を知らない・・・
リンク先にあるように、
xkcd()
とpytyonスクリプトの最初に書いてしまうと、あとはmatplotlibがやってくれるらしいへぇ〜
で、もしかしてこれgnuplotでもできるんじゃない?
と思って「gnuplot xkcd」とかで調べてみたらあったわ・・・
■ 参考 : xkcd-gnuplot
このコードのミソはここで宣言されているjigglea関数jiggle関数
これによってsin(x)とcos(x)に揺らぎを与えて漫画調のぐらぐらした線を再現している
また漫画調のフォントは"Humor Sans"というらしい
ググったら以下のページに落ちていた
■ 参考 : xkcd-desktop/Humor-Sans.ttf
ダウンロードしたttfファイルをダブルクリックすれば、Macなら簡単にインストールできるはず
インストールしたフォントファイルにPATHを通さないとgnuplotでは使えないので注意
PATHの通し方は以下を参考に
■ 過去記事 : 【gnuplot】Could not find/open font...というエラーが出るのでフォントの指定をする
これでこんな感じのグラフが書ける
コードは以下の通り、
(オリジナルからpng出力に変更しただけ)
set term png font 'Humor Sans'
set output "hoge.png"
# Some line types
set style line 10 lt 1 lc rgbcolor "#ffffff" lw 15 #thick white
set style line 11 lt 1 lc rgbcolor "#000000" lw 4 #black
set style line 1 lt 1 lc rgbcolor "#ff0000" lw 4 #red
set style line 2 lt 1 lc rgbcolor "#0000ff" lw 4 #blue
# No border with tics
set border 0
set noxtics
set noytics
# Show the axis
set xzeroaxis ls 11
set yzeroaxis ls 11
#Arrow end to the axis
set arrow from graph 0.95, first 0 to graph 1, first 0 size 2,3 front
set arrow from first 0, graph 0.95 to first 0, graph 1 size 2,3 front
set yrange [-1.1:1.1]
set xrange [-1:15]
set key bottom
set label 'Wasted time' at 11,0.7 right
set arrow from 11.1,0.7 to 12,((12/15.0)**2) ls 11
set title 'Check this out! XKCD in Gnuplot'
#Jiggling functions
range = 2.0 #Range for the absolute jiggle
jigglea(x) = x+range*(2*(rand(0)-0.5)*0.005) #Absolute (as a fraction of the range)
jiggle(x) = x*(1+(2*rand(0)-0.5)*0.015) #Relative +-1.5% (as a fraction of the y value)
dpsin(x) = sin(x) * exp(-0.1 * (x - 5) ** 2)
dpcos(x) = - cos(x) * exp(-0.1 * (x - 5) ** 2)
plot jiggle(dpsin(x)) ls 10 t '', \
jiggle(dpsin(x)) ls 1 t 'Damped Sin',\
jiggle(dpcos(x)) ls 10 t '', \
jiggle(dpcos(x)) ls 2 t 'Damped Cos',\
jigglea((x/15)**2) ls 10 t '',\
jigglea((x/15)**2) ls 11 t ''
set output "hoge.png"
# Some line types
set style line 10 lt 1 lc rgbcolor "#ffffff" lw 15 #thick white
set style line 11 lt 1 lc rgbcolor "#000000" lw 4 #black
set style line 1 lt 1 lc rgbcolor "#ff0000" lw 4 #red
set style line 2 lt 1 lc rgbcolor "#0000ff" lw 4 #blue
# No border with tics
set border 0
set noxtics
set noytics
# Show the axis
set xzeroaxis ls 11
set yzeroaxis ls 11
#Arrow end to the axis
set arrow from graph 0.95, first 0 to graph 1, first 0 size 2,3 front
set arrow from first 0, graph 0.95 to first 0, graph 1 size 2,3 front
set yrange [-1.1:1.1]
set xrange [-1:15]
set key bottom
set label 'Wasted time' at 11,0.7 right
set arrow from 11.1,0.7 to 12,((12/15.0)**2) ls 11
set title 'Check this out! XKCD in Gnuplot'
#Jiggling functions
range = 2.0 #Range for the absolute jiggle
jigglea(x) = x+range*(2*(rand(0)-0.5)*0.005) #Absolute (as a fraction of the range)
jiggle(x) = x*(1+(2*rand(0)-0.5)*0.015) #Relative +-1.5% (as a fraction of the y value)
dpsin(x) = sin(x) * exp(-0.1 * (x - 5) ** 2)
dpcos(x) = - cos(x) * exp(-0.1 * (x - 5) ** 2)
plot jiggle(dpsin(x)) ls 10 t '', \
jiggle(dpsin(x)) ls 1 t 'Damped Sin',\
jiggle(dpcos(x)) ls 10 t '', \
jiggle(dpcos(x)) ls 2 t 'Damped Cos',\
jigglea((x/15)**2) ls 10 t '',\
jigglea((x/15)**2) ls 11 t ''
↓これを書きたい
こんな感じで書ける
ちなみにset term epsでは文字化けしてダメだった
set term postscript enhanced color linewidth 2 fontscale 2
set title '{/Symbol @\326\140}@x{/Symbol \140}@+{/Symbol \140}3'
set ylabel '{/Symbol @\326\140}@x{/Symbol \140}@+{/Symbol \140}3'
set xlabel 'a~{x}{.8\~}'
set output "foo.eps"
plot sin(x)
set title '{/Symbol @\326\140}@x{/Symbol \140}@+{/Symbol \140}3'
set ylabel '{/Symbol @\326\140}@x{/Symbol \140}@+{/Symbol \140}3'
set xlabel 'a~{x}{.8\~}'
set output "foo.eps"
plot sin(x)
単純にルート(sqrt)の記号を表示するだけならば
set xlabel "{/Symbol @\326}"
とかでOK「それだと上にのばし棒がないから冴えなくて困ってる なんか解決方法しらない?」って質問の答えが以下のページであった
■ 参考 : sqrt symbol for sqrt(x+3) and sqrt(x)+3
よー考えるわ・・・
ちなみにx11でも表示されるけど、上記のグラフほどは綺麗ではない
ちなみにpngとして出力すると盛大に文字化けする
上に書いたterminalコマンドにも付いてるけど、enhanced modeってのがある
set term postscript enhanced color
terminalを指定するときに、enhancedと付けておくとさらに拡張機能が使える
今のところこれを付けないデメリットがないのでいつでも付けておけばいいんじゃない?と思うんだがどうなんだろう・・・・
今まではtitleやlabelにアンダーバー(_)や^があると勝手に下付きや上付き文字扱いされて困っていたけど、
gnuplot5になってから
set xlabel "hoge^2" noenhanced
と適宜そのオプションを外せるようになったのでデメリットが思いつかないで、そのenhanced modeの中でチルダも付けられる
set xlabel 'a~{x}{0.7\~}'
xlabelの後はだぶるコーテーションではダメでシングルコーテーションじゃないといけないっぽい
後、aを消すと、なぜか何もlabelが表示されないっぽい
set xlabel ' ~{x}{0.7\~}'
aの代わりにスペースを入れたら動く
1つ目の「~」は重ね書きをするという意味
~の後のxに文字や記号を重ねる
{}の中にその設定などを書く
{0.7}は何倍して上に乗せるかを示す、例えばここを{3\~}とかにすると、x軸を突き抜けて上にチルダが表示される
バックスラッシュはチルダを記号として認識させるためのエスケープ?
2つ目の~は上にチルダを乗せたいからここに書く、例えばこれは-とかでもOK
その場合はTeXで言うところの\bar{x}になる
詳しい説明は以下にある
■ 参考 : gnuplotのTips
GSLをインストールすると、一緒にgsl-histogramというコマンドもインストールされる
■ 使い方 :
もしデータが1列じゃないが、3列目のみのヒストグラムを書きたい場合はawkを使って
実例を示すためにサンプルデータを作る
gnuplotを使ってプロットしてみる
問題はヒストグラムは0から5までなのにboxesの方は-0.5から値があるように描かれていること
理想的にはx=0~1の範囲が3、x=1~2の範囲が2、x=2~3の範囲が2、x=3~4の範囲が3、x=4~5の範囲が0のステップ関数になっていて欲しいので、この場合はstepsの方が求めているものに近い
ただ左端と右端(今回は0だからわからないけど)の挙動がなんか思ってるのと違う・・・
そこも縦棒を描いてほしいんだが・・・
しょうがないのでwith boxesの方にオプションを追加してみる
どういうパターンがあるかは
■ 参考 : gnuplot demo script: fillstyle.dem
■ 参考 : Boxes
同じような塗りつぶしはwith stepsでもできる(最初こっちでやっていて、両端がなんか気に入らなかったので使用をやめた)
なので今回はwith boxesのほうを使う
■ 参考 : filling under step functions in gnuplot
まったく関係ないけど、with histepsというのもある(hiは何の略だろう?)
これをwith stepsの代わりに使うとboxesみたいに半ビンだけ左にずれる
(2016/11/30 追記)
一応上記のwith boxes+囲い線のオプションを使えばこんな感じのグラフも作れる
画像検索でのHIT狙いで追記してみた(スクショだからちょっと線の太さが違う)
(2016/12/22 追記)
ずっと書こう書こうと思って、放置してたけど、ついぞ追記・・・
2種類のヒストグラムを重ね書きするときに半透明にするオプションもあります
(データ生成の方法はパス)
これで
あと、上の方のグラフで+0.5とすべきところを-0.5としてたので修正しておきました
ツイート
■ 使い方 :
Usage: gsl-histogram xmin xmax [n]
Computes a histogram of the data on stdin using n bins from xmin to xmax.
If n is unspecified then bins of integer width are used.
要するにComputes a histogram of the data on stdin using n bins from xmin to xmax.
If n is unspecified then bins of integer width are used.
gsl-histogram 0 100 50
で0から100の間を50ビンで区切ってヒストグラムを計算してくれるもしデータが1列じゃないが、3列目のみのヒストグラムを書きたい場合はawkを使って
awk '{print $3}' | gsl-histogram 0 100 50
とかでおk実例を示すためにサンプルデータを作る
# データ作り
%gsl-randist 0 10 rayleigh 2 > hoge.txt
# データの中身はこれ
% cat hoge.txt
0.0454563
3.81005
3.17951
0.65875
3.42053
2.40609
0.589601
1.53702
2.2201
1.55221
# この行でヒストグラムを取っている、0~5を5ビンで分ける
% awk '{print $1}' hoge.txt | gsl-histogram 0 5 5 > hist.txt
# ヒストグラムの中身、1列目がビンの左端の点、2列目がビンの右端の点、3列目がそのビンに入っているサンプル数
# echo 1 | gsl-histogram 0 2 2 とかで確認できる
% cat hist.txt
0 1 3
1 2 2
2 3 2
3 4 3
4 5 0
%gsl-randist 0 10 rayleigh 2 > hoge.txt
# データの中身はこれ
% cat hoge.txt
0.0454563
3.81005
3.17951
0.65875
3.42053
2.40609
0.589601
1.53702
2.2201
1.55221
# この行でヒストグラムを取っている、0~5を5ビンで分ける
% awk '{print $1}' hoge.txt | gsl-histogram 0 5 5 > hist.txt
# ヒストグラムの中身、1列目がビンの左端の点、2列目がビンの右端の点、3列目がそのビンに入っているサンプル数
# echo 1 | gsl-histogram 0 2 2 とかで確認できる
% cat hist.txt
0 1 3
1 2 2
2 3 2
3 4 3
4 5 0
gnuplotを使ってプロットしてみる
# with boxesで描いてみる
# with boxesは省略不能
plot [-1:5][0:4] "hist.txt" u 1:3 with boxes
または# with boxesは省略不能
plot [-1:5][0:4] "hist.txt" u 1:3 with boxes
# with stepsで描いてみる
# with stと省略可能
plot [-1:5][0:4] "hist.txt" u 1:3 with st
この二つを見比べるとboxesは全体を線で囲って箱にしているのに対して、stepsではそれに沿ってステップ関数を描いているだけ(というかまさに名前がそうなんだけど・・・)# with stと省略可能
plot [-1:5][0:4] "hist.txt" u 1:3 with st
問題はヒストグラムは0から5までなのにboxesの方は-0.5から値があるように描かれていること
理想的にはx=0~1の範囲が3、x=1~2の範囲が2、x=2~3の範囲が2、x=3~4の範囲が3、x=4~5の範囲が0のステップ関数になっていて欲しいので、この場合はstepsの方が求めているものに近い
ただ左端と右端(今回は0だからわからないけど)の挙動がなんか思ってるのと違う・・・
そこも縦棒を描いてほしいんだが・・・
しょうがないのでwith boxesの方にオプションを追加してみる
plot [0:5][0:4] "hist.txt" u ($1+0.5):3 with boxes fs solid 0.25 lc 6 lw 3
他にもplot [0:5][0:4] "hist.txt" u ($1+0.5):3 with boxes fill pattern 4 lc 6
とかもできるどういうパターンがあるかは
set term pnt
set output "hoge.png"
test
で確認できるset output "hoge.png"
test
■ 参考 : gnuplot demo script: fillstyle.dem
■ 参考 : Boxes
同じような塗りつぶしはwith stepsでもできる(最初こっちでやっていて、両端がなんか気に入らなかったので使用をやめた)
plot [-1:5][0:4] "hist.txt" u 1:3 with fillsteps fs solid 0.25 noborder ls 7, "" u 1:3 with steps ls 7 lw 3
やっぱり左端(見えないけど右端も)の扱いがなんか嫌・・・なので今回はwith boxesのほうを使う
■ 参考 : filling under step functions in gnuplot
まったく関係ないけど、with histepsというのもある(hiは何の略だろう?)
これをwith stepsの代わりに使うとboxesみたいに半ビンだけ左にずれる
plot [-1:5][0:4] "hist.txt" u 1:3 with fillsteps fs solid 0.25 noborder ls 7, "" u 1:3 with histeps ls 7 lw 3
(2016/11/30 追記)
一応上記のwith boxes+囲い線のオプションを使えばこんな感じのグラフも作れる
画像検索でのHIT狙いで追記してみた(スクショだからちょっと線の太さが違う)
(2016/12/22 追記)
ずっと書こう書こうと思って、放置してたけど、ついぞ追記・・・
2種類のヒストグラムを重ね書きするときに半透明にするオプションもあります
(データ生成の方法はパス)
set style fill transparent solid 0.15
set term png
set output "hoge.png"
plot [0:20]"hist_2.txt" u ($1+0.5):3 with boxes fs lc 7 lw 3 title "rayleigh 2", "hist_5.txt" u ($1+0.5):3 with boxes fs lc 6 lw 3 title "rayleigh 5"
set term png
set output "hoge.png"
plot [0:20]"hist_2.txt" u ($1+0.5):3 with boxes fs lc 7 lw 3 title "rayleigh 2", "hist_5.txt" u ($1+0.5):3 with boxes fs lc 6 lw 3 title "rayleigh 5"
これで
あと、上の方のグラフで+0.5とすべきところを-0.5としてたので修正しておきました
何がやりたいかというとこういうグラフが書きたい
左下がメインの図で、2つのデータの散布図
左上が散布図のx方向のデータのヒストグラム、右下が散布図のy軸方向のデータのヒストグラム
おまけで平均値まで付けてみた(awkさん、ありがとう)
普通にgnuplotのmultiplotで2x2の分割を使っていると左上→右上→左下→右下という順番でグラフが埋まっていく
けど今回は右上を空けておきたい
そういう白紙にしてスキップするようなコマンドがあるのかマニュアル読んだけど、
英語でわからんちんだったのでゴリ押すことにした
あとgnuplotのdemoも見たけど、multiplotのページは全部埋まってる図しかなかった
(2016/11/5 追記)
右下の図を少し回転させたのが次の図で、下のスクリプトはこの図を書くものに修正しました
赤字が少し工夫した点
■ 画面に白紙を表示させる方法
まず枠(border)やgrid、tics(グラフに付いてる目盛り)はすべて手で消さないといけない
それらを表示させない設定をした後で、白紙の画面を描画→設定を元に戻す必要がある
そのときにresetとしたくなるが3番目、4番目の図の大きさが他の2枚と一致しなくなるのでやめた方がいい
具体的には
最後の plot [:][-1:1]1/0はy軸の範囲指定をきちんとしないとエラーになるっぽい
他に出てきたのはplot sqrt(-1)とか
というかプロットする範囲を指定した上で、そこに入らないように何かをプロットするならなんでもいいと思われる
■ データの作り方
なぜdata1がないのか気になってしまうが何かしらの試行錯誤の名残なんだと思う・・・・(シャッと書いたので気にしない)
gsl-randistというコマンドはGSLをインストールしたら一緒に入ってくると思う
そのコマンドがない人はなんでもいいので1次元のデータを手で書いてください
data0というのをわざわざ作っているのは散布図を作っているから散布図を書きたいから
■ プロットのスクリプト
ツイート
左下がメインの図で、2つのデータの散布図
左上が散布図のx方向のデータのヒストグラム、右下が散布図のy軸方向のデータのヒストグラム
おまけで平均値まで付けてみた(awkさん、ありがとう)
普通にgnuplotのmultiplotで2x2の分割を使っていると左上→右上→左下→右下という順番でグラフが埋まっていく
けど今回は右上を空けておきたい
そういう白紙にしてスキップするようなコマンドがあるのかマニュアル読んだけど、
英語でわからんちんだったのでゴリ押すことにした
あとgnuplotのdemoも見たけど、multiplotのページは全部埋まってる図しかなかった
(2016/11/5 追記)
右下の図を少し回転させたのが次の図で、下のスクリプトはこの図を書くものに修正しました
赤字が少し工夫した点
■ 画面に白紙を表示させる方法
まず枠(border)やgrid、tics(グラフに付いてる目盛り)はすべて手で消さないといけない
それらを表示させない設定をした後で、白紙の画面を描画→設定を元に戻す必要がある
そのときにresetとしたくなるが3番目、4番目の図の大きさが他の2枚と一致しなくなるのでやめた方がいい
具体的には
set key off
set xlabel ""
set ylabel ""
unset ytics
unset xtics
set nomy2tics
set nomx2tics
unset border
unset grid
plot [:][-1:1]1/0
とかでおkset xlabel ""
set ylabel ""
unset ytics
unset xtics
set nomy2tics
set nomx2tics
unset border
unset grid
plot [:][-1:1]1/0
最後の plot [:][-1:1]1/0はy軸の範囲指定をきちんとしないとエラーになるっぽい
他に出てきたのはplot sqrt(-1)とか
というかプロットする範囲を指定した上で、そこに入らないように何かをプロットするならなんでもいいと思われる
■ データの作り方
data2=hoge.txt
data3=foo.txt
gsl-randist $seed $nsample rayleigh 1 > $data2
gsl-randist $seed $nsample chisq 5 > $data3
hist2=hist_2.txt
ave_2=`awk '{sum+=$1}END{print sum/NR}' ${data2} | cut -c1-4`
max_2=`awk 'NR==1 {max=$1} {if($1 > max) max = $1} END {print max}' ${data2}`
cat ${data2} | gsl-histogram $xmin ${max_2} $nbin > ${hist2}
echo "$ave_2"
hist3=hist_3.txt
ave_3=`awk '{sum+=$1}END{print sum/NR}' ${data3} | cut -c1-4`
max_3=`awk 'NR==1 {max=$1} {if($1 > max) max = $1} END {print max}' ${data3}`
cat ${data3} | gsl-histogram $xmin ${max_3} $nbin > ${hist3}
echo "$ave_3"
data0=paste.txt
paste $data2 $data3 > $data0
これでデータができるdata3=foo.txt
gsl-randist $seed $nsample rayleigh 1 > $data2
gsl-randist $seed $nsample chisq 5 > $data3
hist2=hist_2.txt
ave_2=`awk '{sum+=$1}END{print sum/NR}' ${data2} | cut -c1-4`
max_2=`awk 'NR==1 {max=$1} {if($1 > max) max = $1} END {print max}' ${data2}`
cat ${data2} | gsl-histogram $xmin ${max_2} $nbin > ${hist2}
echo "$ave_2"
hist3=hist_3.txt
ave_3=`awk '{sum+=$1}END{print sum/NR}' ${data3} | cut -c1-4`
max_3=`awk 'NR==1 {max=$1} {if($1 > max) max = $1} END {print max}' ${data3}`
cat ${data3} | gsl-histogram $xmin ${max_3} $nbin > ${hist3}
echo "$ave_3"
data0=paste.txt
paste $data2 $data3 > $data0
なぜdata1がないのか気になってしまうが何かしらの試行錯誤の名残なんだと思う・・・・(シャッと書いたので気にしない)
gsl-randistというコマンドはGSLをインストールしたら一緒に入ってくると思う
そのコマンドがない人はなんでもいいので1次元のデータを手で書いてください
data0というのをわざわざ作っているのは散布図を作っているから散布図を書きたいから
■ プロットのスクリプト
#!/bin/sh
xmin=0
xmax=10
nbin=40
nsample=10000
seed=0
data2=hoge.txt
data3=foo.txt
gsl-randist $seed $nsample rayleigh 1 > $data2
gsl-randist $seed $nsample chisq 5 > $data3
hist2=hist_2.txt
ave_2=`awk '{sum+=$1}END{print sum/NR}' ${data2} | cut -c1-4`
max_2=`awk 'NR==1 {max=$1} {if($1 > max) max = $1} END {print max}' ${data2}`
cat ${data2} | gsl-histogram $xmin ${max_2} $nbin > ${hist2}
echo "$ave_2"
hist3=hist_3.txt
ave_3=`awk '{sum+=$1}END{print sum/NR}' ${data3} | cut -c1-4`
max_3=`awk 'NR==1 {max=$1} {if($1 > max) max = $1} END {print max}' ${data3}`
cat ${data3} | gsl-histogram $xmin ${max_3} $nbin > ${hist3}
echo "$ave_3"
data0=paste.txt
paste $data2 $data3 > $data0
gnuplot <<EOF
set term png size 720, 720 fontscale 1.2 linewidth 1.5
set output "hoge.png"
set tmargin 0
set multiplot layout 2,2
set size square
set lmargin 6
set rmargin 4
set tmargin 0
set bmargin 2
## (1, 1)
set xlabel "rayleigh dist"
plot [${xmin}:${max_2}] "$hist2" u 1:3 w step lw 2 lc 1 title "mean=${ave_2}"
## (1, 2)
set key off
set xlabel ""
set ylabel ""
unset ytics
unset xtics
set nomy2tics
set nomx2tics
unset border
unset grid
plot [:][-1:1]1/0
## (2, 1)
set xtics rotate by -60
set key on
set ytics
set xtics mirror
set border
set grid
set bmargin 4
set xlabel "rayleigh dist"
set ylabel "chi^2 dist"
plot [${xmin}:${max_2}][${xmin}:${max_3}] "$data0" u 1:2 w d lc 0 notitle
## (2, 2)
set xlabel ""
set ylabel "chi^2 dist"
plot [:] [${xmin}:${max_3}] "$hist3" u 3:1 w step lw 2 lc 3 title "mean=$ave_3"
unset multiplot
EOF
rm ./*.txt
xmin=0
xmax=10
nbin=40
nsample=10000
seed=0
data2=hoge.txt
data3=foo.txt
gsl-randist $seed $nsample rayleigh 1 > $data2
gsl-randist $seed $nsample chisq 5 > $data3
hist2=hist_2.txt
ave_2=`awk '{sum+=$1}END{print sum/NR}' ${data2} | cut -c1-4`
max_2=`awk 'NR==1 {max=$1} {if($1 > max) max = $1} END {print max}' ${data2}`
cat ${data2} | gsl-histogram $xmin ${max_2} $nbin > ${hist2}
echo "$ave_2"
hist3=hist_3.txt
ave_3=`awk '{sum+=$1}END{print sum/NR}' ${data3} | cut -c1-4`
max_3=`awk 'NR==1 {max=$1} {if($1 > max) max = $1} END {print max}' ${data3}`
cat ${data3} | gsl-histogram $xmin ${max_3} $nbin > ${hist3}
echo "$ave_3"
data0=paste.txt
paste $data2 $data3 > $data0
gnuplot <<EOF
set term png size 720, 720 fontscale 1.2 linewidth 1.5
set output "hoge.png"
set tmargin 0
set multiplot layout 2,2
set size square
set lmargin 6
set rmargin 4
set tmargin 0
set bmargin 2
## (1, 1)
set xlabel "rayleigh dist"
plot [${xmin}:${max_2}] "$hist2" u 1:3 w step lw 2 lc 1 title "mean=${ave_2}"
## (1, 2)
set key off
set xlabel ""
set ylabel ""
unset ytics
unset xtics
set nomy2tics
set nomx2tics
unset border
unset grid
plot [:][-1:1]1/0
## (2, 1)
set xtics rotate by -60
set key on
set ytics
set xtics mirror
set border
set grid
set bmargin 4
set xlabel "rayleigh dist"
set ylabel "chi^2 dist"
plot [${xmin}:${max_2}][${xmin}:${max_3}] "$data0" u 1:2 w d lc 0 notitle
## (2, 2)
set xlabel ""
set ylabel "chi^2 dist"
plot [:] [${xmin}:${max_3}] "$hist3" u 3:1 w step lw 2 lc 3 title "mean=$ave_3"
unset multiplot
EOF
rm ./*.txt
プロフィール
HN:coffee
職業:物理屋(自称)
趣味:映画鑑賞、登山
出身:大阪府の南の田舎
自己紹介:
import MyProfile
import coffee_pote from TWITTER
import amazonのほしい物リスト from WISH_LIST
print "先月子供が産まれました!"
# 最終更新 2022/10/25
職業:物理屋
趣味:映画鑑賞、登山
出身:大阪府の南の田舎
自己紹介:
import MyProfile
import coffee_pote from TWITTER
import amazonのほしい物リスト from WISH_LIST
print "先月子供が産まれました!"
# 最終更新 2022/10/25
カテゴリー
ブログ内検索
リンク
相互リンク募集中です
(Twitterにてお知らせください)
Demo scripts for gnuplot version 5 (gnuplotのさまざまなデモ画像と作り方がまとめられている、眺めているだけでできるようになった気分になれる)
gnuplotスクリプトの解説 (米澤進吾さんの個人ページ、gnuplotと言えばこのかた)
gnuplot のページ (Takeno Lab、うちのブログがリンクされていたのでリンク返し)
Twitterから映画の評価が分かる & 映画の鑑賞記録が残せる coco (映画の感想をまとめられるサイト、いつもお世話になっています)
Astronomy Picture of the Day Archive (天文や宇宙関連の最新の話題について画像とともにNASAが説明しているページ)
今日のほしぞら (任意の時刻の空で見える星を表示してくれる、国立天文台が管理している)
GNUPLOTとアニメーション (応用の項目の「見せてあげよう!ラピュタの雷を!!」あたりからすごすぎる)
読書メーター (読んだ本をリストできる便利なサイト)
flickr難民の写真置き場 (20XX年、flickrは有料化の炎に包まれた。あらゆるflickr無料ユーザーは絶滅したかに見えた。 しかし、tumblr移住民は死に絶えてはいなかった。)
教授でもできるMac OS X へのLaTeX, X11, gccのインストレーションと環境設定 (阪大の山中卓さんのwebページ、タイトルにセンスが溢れている、内容は超充実してる、特にTeX関連、学振DCとかPDの申請書類作成時にはお世話になっております)
英語論文執筆用の例文検索サービス (とんでもないものを見つけてしまった・・・・ arXivに収録されている 811,761報の 英語論文から,例文を検索するための検索エンジン)
Demo scripts for gnuplot version 5 (gnuplotのさまざまなデモ画像と作り方がまとめられている、眺めているだけでできるようになった気分になれる)
gnuplotスクリプトの解説 (米澤進吾さんの個人ページ、gnuplotと言えばこのかた)
gnuplot のページ (Takeno Lab、うちのブログがリンクされていたのでリンク返し)
Twitterから映画の評価が分かる & 映画の鑑賞記録が残せる coco (映画の感想をまとめられるサイト、いつもお世話になっています)
Astronomy Picture of the Day Archive (天文や宇宙関連の最新の話題について画像とともにNASAが説明しているページ)
今日のほしぞら (任意の時刻の空で見える星を表示してくれる、国立天文台が管理している)
GNUPLOTとアニメーション (応用の項目の「見せてあげよう!ラピュタの雷を!!」あたりからすごすぎる)
読書メーター (読んだ本をリストできる便利なサイト)
flickr難民の写真置き場 (20XX年、flickrは有料化の炎に包まれた。あらゆるflickr無料ユーザーは絶滅したかに見えた。 しかし、tumblr移住民は死に絶えてはいなかった。)
教授でもできるMac OS X へのLaTeX, X11, gccのインストレーションと環境設定 (阪大の山中卓さんのwebページ、タイトルにセンスが溢れている、内容は超充実してる、特にTeX関連、学振DCとかPDの申請書類作成時にはお世話になっております)
英語論文執筆用の例文検索サービス (とんでもないものを見つけてしまった・・・・ arXivに収録されている 811,761報の 英語論文から,例文を検索するための検索エンジン)
最新記事
(11/20)
(03/05)
(02/29)
(02/21)
(02/21)
(02/21)
(02/21)
(01/13)
(01/05)
(01/05)