# Choosing center of histogram bins for fitting

I have a bimodal distribution, and if plotted with Mathematica it looks like this:

Now, the lowest value from the actual data is 8196 and 690720, but as seen in the plot, Mathematica lets the data range go from 0 to 744572. Is Mathematica choosing a bad histogram data range?

What is in general a good choice for defining the center of the bins and data range, so that I can fit a distribution through the histogram?

My approach would be:

(1) Calculate (bin width) = (Max-Min)/(number of bins)
[I'm aware that there are different rules how to choose the optimal number of
bins depending on the underlying distribution, let's just assume this is 12]

(2) then I have 12 equal bins, starting from 8196 and ending at 690720, each having a width of 56877

(3) The first bin goes then from 8196 to 65073=8196+56877 and so on

(4) As the center of the bin I define the middle between 8196 and 65073 which is 36624.5 and I position my first bin there.

(5) Then I get 12 data pairs of bin center position and number of observations and I can fit a bimodal distribution through it


Am I making a mistake if I do that, or what is the reasoning behind Mathematica's choice of the histogram range exceeding the actual data range?

Edit: I've uploaded the raw data here: raw data

Edit2: For clarity explaining the mysterious frequency of 17 which was pointed out by Stephan: The data is a confocal photoluminescence map where a laser scans an emitter and it looks like this: The laser scans row by row so in the middle of each row the emitter lights up which explains the frequency of 17 when the raw data is plotted as it originates from a single list.

• Fit a distribution to the raw data, not bin summaries. From that point of view, histogram binning is irrelevant. It's hard to know what named distribution would make sense here, but perhaps you plan to fit a mixture. The Mathematica default is one of many possible reasonable choices. At a wild guess you have ~600 data points, so many researchers would favour or tolerate many more bins. Most advice is centred on the number of bins, but it can make as much or more sense to think about what bin width makes sense scientifically or practically. Jun 29 '18 at 8:04
• Questions on (1) how to fit a (bimodal) distribution (2) good defaults or criteria for histogram binning are in my view pretty much unrelated. Jun 29 '18 at 8:06
• I've added the raw data (see edit). If I use the raw data (binning with a bin width of 1, then the highest occurrences are 2. Or what exactly did you mean? There are 799 data points and I know I need more bins, the number of 12 was just chosen for simplicity.
– user213095
Jun 29 '18 at 8:19
• The plan is to fit a mixture of 2 normal distributions (see also en.wikipedia.org/wiki/…), I'm interested in the center position of the second peak which I want to extract, as well as the width of the second peak.
– user213095
Jun 29 '18 at 8:25
• A mixture of two normals is wrong in principle, unless you add truncation. No normal distribution starts with finite density at zero. I'd use kernel density estimation to get a better idea of your distribution. Jun 29 '18 at 8:50

As @Nick Cox says, fit your distribution directly to the data. Do not first bin the data into a histogram. Why would you want to do so?

Instead, fit a standard kernel density. I'll use R, because I know it better, but I assume Mathematica has similar functionalities. (If it doesn't, I recommend you learn R.) Below is code that will fit such a density to your data and extract the $x$ value for the second peak.

For added enlightenment, we can assess how sure we are of this coordinate by bootstrapping it. I'm also plotting a bootstrapped 95% quantile. Notice how this is slightly asymmetrical.

dataset <- unlist(read.table("https://files.fm/down.php?i=qucxqxgw"))

foo <- density(dataset)
max.index <- which(foo$x>4e5)[which.max(foo$y[foo$x>4e5])] plot(foo) points(foo$x[max.index],foo$y[max.index],pch=19,col="red") text(foo$x[max.index],foo$y[max.index],round(foo$x[max.index]),pos=3,col="red")

library(boot)
bootstrap <- boot(dataset,statistic=function(dataset,index){
foo <- density(dataset[index])
max.index <- which(foo$x>4e5)[which.max(foo$y[foo$x>4e5])] foo$x[max.index]
}, R=1e3)

lines(quantile(bootstrap$t,c(0.025,0.975)),rep(foo$y[max.index],2),col="red",lwd=2)


If you want the width of the second peak, you can extract it from the density (and bootstrap it) after you have decided how you define a peak (anything more than 95% of the peak value, or a fixed offset, or something else).

(Yes, in principle we could correct for the fact that your data seem to be all nonnegative, whereas the density estimate goes negative. In practice, since you are only interested in the second peak, I don't really see the point.)

# However...

Here is a plot of your original data:

plot(dataset,type="o")


This looks strangely regular. So, after playing around a bit with the frequency parameter, we find the following seasonplot:

library(forecast)
seasonplot(ts(dataset,frequency=17))


Unless you did some very strange sorting to your raw data, your data is actually seasonal with a period of 17. Thus, I'd question whether finding the location of your second mode in such data is really what you want to be doing at all.

• I was drafting a partly similar answer using Stata, which is now mostly redundant, so I won't post it. Good catch on the periodicity! I will add that the smoothing of the density function below zero may not make physical sense if the variable must be positive. That's more than a detail because arguably the density so smoothed below 0 should be reflected back, which may shift the estimate of the lower (main) mode. Jun 29 '18 at 9:07
• @NickCox: you are right about constraining the density to be nonnegative, but I kind of doubt it will have a major impact on the second mode, in particular given the uncertainty in bandwidth selection... Jun 29 '18 at 9:13
• The data is a 17x47 matrix. Actually, the data is a map of different intensity values generated from a confocal photoluminescence map, where a laser scans over an emitter. And this is done 17 times at different locations, hence the periodicity of 17. I'll add a picture in the original post for clarity. I think we can ignore that fact, as this sorting is just an artifact of how the data recorded.
– user213095
Jun 29 '18 at 9:21
• Yes. I would like to isolate the bright intensities of the 2D map, but instead of using a simple threshold a better solution would be getting the second peak position in the histogram and use this as a threshold value.
– user213095
Jun 29 '18 at 9:33
• If you just want all the medium-bright intensities, then you can proceed along the answer I provide. (Getting the very brightest ones would mean looking at the highest mode.) Alternatively, you may want to look at two-dimensional kernel density estimates. Jun 29 '18 at 9:37

You can use Mathematica to draw a SmoothHistogram. I use here the data you provided in a flattened list.

data = Flatten[Import["raw.txt", "TSV"]]

Histogram[data]

SmoothHistogram[data]

Now assume you want to fit to a mixture of two distributions, and obtain the parameters and also the mixture, first define a distribution mixture:

distMix = MixtureDistribution[{p, 1 - p}, {NormalDistribution[a, b], NormalDistribution[c, d]}]

And then, obtain the parameters:

params = FindDistributionParameters[data, distMix]

{p -> 0.382944, a -> 516841., b -> 101764., c -> 124503., d -> 81260.4}

• Also a very good answer. But the parameters found are not the peak positions. Does this mean that this mixture of normals is bad and different distributions should be used? Or how can I extract the peak positions without fitting some distribution to SmoothHistogram?
– user213095
Jul 1 '18 at 8:07
• If I let Mathematica itself find the distribution the result is better, but still not close to the actual result
– user213095
Jul 1 '18 at 8:12
• I assume you are using FindDistribution to obtain the MixtureDistribution automatically. The parameters depend on the data and choice of distributions in the mixture. I am not familiar with your data, and not sure what the "actual result" means - that should be whatever the model gives, but I suggest you use different mixtures, and then you can also use DistributionFitTest to check these for a fit.
– GIM
Jul 5 '18 at 16:42