Нормализация бинов гистограммы в gnuplot

В gnuplot 4.4 функции имеют другое свойство: они могут выполнять несколько последовательных команд, а затем возвращать значение (см. уловки gnuplot). Это означает, что вы действительно можете рассчитать количество точек n в файле gnuplot, даже не зная об этом заранее. Этот код работает для файла out.dat, содержащего один столбец: список из n образцов из нормального распределения:

binwidth = 0.1
set boxwidth binwidth
sum = 0

s(x)          = ((sum=sum+1), 0)
bin(x, width) = width*floor(x/width) + binwidth/2.0

plot "out.dat" u ($1):(s($1))
plot "out.dat" u (bin($1, binwidth)):(1.0/(binwidth*sum)) smooth freq w boxes

Первый оператор построения графика читает файл данных и увеличивает сумму один раз для каждой точки, выводя на график ноль.

Второй оператор графика фактически использует значение суммы для нормализации гистограммы.

В gnuplot 4.6 вы можете подсчитать количество очков с помощью команды stats, которая быстрее, чем plot. На самом деле, вам не нужен такой трюк s(x)=((sum=sum+1),0), а просто посчитайте число по переменной STATS_records после запуска stats 'out.dat' u 1.

Вот как я поступил бы с n = 500 случайными гауссовскими переменными, сгенерированными из R с помощью следующей команды:

Rscript -e 'cat(rnorm(500), sep="\\n")' > rnd.dat

Я использую ту же идею, что и ваша, для определения нормализованной гистограммы, где y определяется как 1 / (binwidth * n), за исключением того, что я использую int вместо floor, и я не повторно центрировал значение bin. Короче говоря, это быстрая адаптация демонстрационного скрипта smooth.dem и аналогичный подход описан в учебнике Джанерта Gnuplot в действии (Глава 13, стр. 257, в свободном доступе). Вы можете заменить мой образец файла данных на random-points, который находится в папке demo, поставляемой с Gnuplot. Обратите внимание, что нам нужно указать количество точек как Gnuplot, поскольку нет средств подсчета для записей в файле.

bw1=0.1
bw2=0.3
n=500
bin(x,width)=width*int(x/width)
set xrange [-3:3]
set yrange [0:1]
tstr(n)=sprintf("Binwidth = %1.1f\n", n) 
set multiplot layout 1,2
set boxwidth bw1
plot 'rnd.dat' using (bin($1,bw1)):(1./(bw1*n)) smooth frequency with boxes t tstr(bw1)
set boxwidth bw2
plot 'rnd.dat' using (bin($1,bw2)):(1./(bw2*n)) smooth frequency with boxes t tstr(bw2)

Вот результат с двумя ячейками шириной

Кроме того, это действительно грубый подход к гистограмме, и более сложные решения доступны в R. Действительно, проблема в том, как определить подходящую ширину бина, и этот вопрос уже обсуждался на stats.stackexchange.com: использование правила объединения Freedman-Diaconis не должно быть слишком сложным для реализации, хотя вам нужно будет вычислить межквартильный диапазон.

Вот как R будет действовать с тем же набором данных, с опцией по умолчанию (правило Стерджеса, потому что в этом конкретном случае это не будет иметь значения) и с одинаковым интервалом, как те, которые использовались выше.

Используемый код R приведен ниже:

par(mfrow=c(1,2), las=1)
hist(rnd, main="Sturges", xlab="", ylab="", prob=TRUE)
hist(rnd, breaks=seq(-3.5,3.5,by=.1), main="Binwidth = 0.1", 
     xlab="", ylab="", prob=TRUE)

Вы даже можете посмотреть, как R выполняет свою работу, проверив значения, возвращаемые при вызове hist():

> str(hist(rnd, plot=FALSE))
List of 7
 $ breaks     : num [1:14] -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 ...
 $ counts     : int [1:13] 1 1 12 20 49 79 108 87 71 43 ...
 $ intensities: num [1:13] 0.004 0.004 0.048 0.08 0.196 0.316 0.432 0.348 0.284 0.172 ...
 $ density    : num [1:13] 0.004 0.004 0.048 0.08 0.196 0.316 0.432 0.348 0.284 0.172 ...
 $ mids       : num [1:13] -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25 0.75 1.25 ...
 $ xname      : chr "rnd"
 $ equidist   : logi TRUE
 - attr(*, "class")= chr "histogram"

Все это говорит о том, что вы можете использовать результаты R для обработки ваших данных с помощью Gnuplot, если хотите (хотя я бы рекомендовал использовать R напрямую :-).

Другой способ подсчета количества точек данных в файле - использование системной команды. Это оказывается полезным, если вы строите несколько файлов и заранее не знаете количество точек. Я использовал:

countpoints(file) = system( sprintf("grep -v '^#' %s| wc -l", file) )
file1count = countpoints (file1)
file2count = countpoints (file2)
file3count = countpoints (file3)
...

Функции countpoints избегают подсчета строк, начинающихся с символа "#". Затем вы должны использовать уже упомянутые функции для построения нормализованной гистограммы.

Вот полный пример:

n=100
xmin=-50.
xmax=50.
binwidth=(xmax-xmin)/n

bin(x,width)=width*floor(x/width)+width/2.0
countpoints(file) = system( sprintf("grep -v '^#' %s| wc -l", file) )

file1count = countpoints (file1)
file2count = countpoints (file2)
file3count = countpoints (file3)

plot file1 using (bin(($1),binwidth)):(1.0/(binwidth*file1count)) smooth freq with boxes,\
     file2 using (bin(($1),binwidth)):(1.0/(binwidth*file2count)) smooth freq with boxes,\
     file3 using (bin(($1),binwidth)):(1.0/(binwidth*file3count)) smooth freq with boxes
...

Просто

plot 'file' using (bin($2, binwidth)):($4/$4) smooth freq with boxes