Here is how I would do it, with n = 500 random Gaussian variations generated from R using the following command:
Rscript -e 'cat(rnorm(500), sep="\\n")' > rnd.dat
I use the same idea as yours to define a normalized histogram where y is defined as 1 / (bin width * n), except that I use int instead of floor , and I did not return to the bin value. In short, this is a quick adaptation from smooth.dem demo script, and a similar approach is described in the Janert, Gnuplot tutorial in action ( Chapter 13 , p. 257, freely available). You can replace the random-points sample data file, which is available in the demo folder that comes with Gnuplot. Please note that we need to specify the number of points in the form of Gnuplot as there are no counters for entries in the file.
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)
Here is the result, with two hopper widths
Moreover, this is a really rough approach to the histogram, and more complex solutions are easily accessible in R. Indeed, the problem is how to determine a good bean width, and this problem has already been discussed at stats.stackexchange.com : using the Freedman-Diaconis rule is not should be too complicated to implement, although you will need to calculate the range between quartiles.
This is how R will work with the same dataset, with a default parameter (Sturges rule, because in this particular case it will not make any difference) and an equally distributed bin, like the ones used above.
The R code used is shown below:
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)
You can even see how R does its work by checking the values returned by 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"
All that can be said is that you can use the results of R to process your data with Gnuplot if you want (although I would recommend using R: - directly).