Нормализация гистограммы бинов в gnuplot
Я пытаюсь построить гистограмму, чьи ячейки нормализованы по количеству элементов в ячейке.
Я использую следующее
binwidth=5
bin(x,width)=width*floor(x/width) + binwidth/2.0
plot 'file' using (bin($2, binwidth)):($4) smooth freq with boxes
чтобы получить базовую гистограмму, но я хочу, чтобы значение каждого бина было разделено на размер бина. Как я могу сделать это в gnuplot или при необходимости использовать внешние инструменты?
5 ответов
В 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
Первый оператор plot считывает файл данных и увеличивает сумму по одному разу для каждой точки, представляя ноль.
Второе утверждение графика фактически использует значение суммы для нормализации гистограммы.
В 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
Я использую ту же идею, что и вы, для определения нормализованной гистограммы, где у определяется как 1/(binwidth * n), за исключением того, что я использую int
вместо floor
и я не отказался от стоимости бина. Короче говоря, это быстрая адаптация из демонстрационного сценария smooth.dem, и аналогичный подход описан в учебнике Джанерта " Gnuplot in Action" ( глава 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- Правило биннинга Диакониса не должно быть слишком сложным для реализации, хотя вам нужно будет вычислить интервал квартилей.
Вот как R будет работать с тем же набором данных, с опцией по умолчанию (правило Sturges, потому что в данном конкретном случае это не будет иметь значения) и с одинаково разнесенным бином, как те, что использовались выше.
Код 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