gnuplotで正規分布を表示

gnuplotで正規分布のグラフをterminal上に表示する。環境はLinux(CentOS 6)

[手順]

1. perlで正規分布に従う乱数のリストをファイルに出力
2. gnuplotでterminalにグラフ出力

事前準備(perl, gnulpotのinstall)

$ sudo yum -y install gnuplot
$ sudo yum -y install perl

正規分布に従う乱数を出力するスクリプト
ndist.pl

#!/usr/bin/perl

$Pi = 3.14159265359;

$input = 1.0;
$sigma = 0.01;
$count = 10000;

while ($count > 0) {
    printf("%lf\n", RandBoxMuller($input, $sigma));
    $count--;
}

sub RandBoxMuller {
    my ($m, $sigma) = @_;
    my ($r1, $r2) = (rand(), rand());

    while ($r1 == 0) { $r1 = rand(); }

    return ($sigma * sqrt(-2 * log($r1)) * sin(2 * $Pi * $r2)) + $m;
}

gnuplot用スクリプト
plothist.gp

reset

# set parameters
interval=100
min=0.950
max=1.050
width=(max-min)/interval

#function used to map a value to the intervals
hist(x,width)=width*floor(x/width)+width/2.0
set terminal dumb
set xrange [min:max]
set yrange [0:500]
set xtics min,(max-min)/5,max
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
set tics out nomirror
set xlabel "x"
set ylabel "Frequency"

#calculate and plot
plot "test.dat" u (hist($1,width)):(1.0) smooth freq w boxes

実行用シェルスクリプト
show_histgram.sh

#!/bin/sh

/usr/bin/perl ./ndist.pl > test.dat
/usr/bin/gnuplot plothist.gp

実行結果

$ ./show_histgram.sh



  Frequency
    500 ++------------------------------------------------------------------+
         |                       "test.dat" u (hist($1,width)):(1.0) ****** |
         |                                                                  |
         |                               ***                                |
    400 ++                               *****                              |
         |                             ********                             |
         |                            *********                             |
    300 ++                           ************                           |
         |                          *************                           |
         |                          ***************                         |
         |                         ****************                         |
    200 ++                        ******************                        |
         |                       *********************                      |
         |                      ***********************                     |
    100 ++                     ************************                     |
         |                    ***************************                   |
         |                  *******************************                 |
         |               ************************************               |
      0 ++--------*-************************************************-*------+
         +            +             +            +             +            +
        0.95         0.97          0.99         1.01          1.03         1.05
                                          x
$