How to model a probability distribution to phone call duration? I was plotting the phone call durations, and the result seems like a right-skewed histogram. Shape-wise,  Poisson, Raylegh, Chi-2 can be candidates for modeling that.
But is there a more rigorous way of modeling? Like a real-world intuition behind one of those options? In addition, how to tune the parameters to fit the curve? Do you know a tool for that? Should I use moments of the distribution?
Any suggestion is welcomed. I have attached the histogram.

 A: 
But is there a more rigorous way of modeling?

Visually comparing a theoretical distribution to an observed histogram is an important practice toward being rigorous in our inferences.


Like a real-world intuition behind one of those options?

This would require domain-specific knowledge about phone calls that perhaps I don't have. There are multiple distributions that appear similar to yours prima facie, and I do not know which specific theoretical distribution is 'best' according to my background information. But perhaps there isn't an available way to know other than to evaluate and compare models trained on the data. Some speculations I might have about the general shape:

*

*Unimodality is suggesting there is no evidence here of clusters of phone calls.

*Hypothesis: The tiny strip near zero might be those automatic calls along with the occassional "oh s**t I just dialed the wrong number" calls.

*Most of the rest of the weight not being close to zero even on the left-side of the mode might be suggestive that for most calls there is a minimum length that is useful to humans for verbally getting information across.

*The right skew suggests that occasional calls go much longer than the typical call length. For me these longer calls are catching up with old friends that I don't see often, but it could be for other reasons in general.

*Sanity check: Phone calls should not have negative duration, and the histogram agrees with that constraint.



In addition, how to tune the parameters to fit the curve?

There are too many ways to summarize for how to fit the parameters of a distribution to data. The choices often depend on whether you are approaching this from a Bayesian or likelihood perspective. Let me mention two estimation techniques in these broad areas, though I cannot emphasize enough that this barely scratches the surface.

*

*Maximum a posteriori estimation

*Maximum likelihood estimation


Do you know a tool for that?

For the simple cases, which yours might be, I often use SciPy to fit distributions to data. By default it uses maximum likelihood estimation, but method of moments can be used by passing method="MM" into the fit method. You could use the code in this answer to train (all?) the available SciPy distributions on your data, but this seems wasteful (if not asking for spurious results).
Taking a more practical stance of training only specific distributions, let's consider how you can do this with the gamma distribution as an example.
from scipy.stats import gamma

data = <an array of your data goes here>

parameters = gamma.fit(data) # include `floc=0` if you don't want the location parameter.

A list of available statistical functions is found in this reference.
If SciPy doesn't have a distribution you are interested in, like if you wanted to construct mixture model, you can build your own using using scipy.stats.rv_continuous or scipy.stats.rv_discrete as base classes. See this Stack Overflow Q/A for an example.
In your case you can probably model you distribution as a mixture distribution between a Dirac delta distribution centered at zero with something like gamma distribution. SciPy's base class should work for that with some customized code, but for building custom mixtures I have become a fan of pymc.Mixture.


Should I use moments of the distribution?

If the best-fitting distribution has finite moments to some sufficiently large order, then you can use them to report typical moments such as the mean, variance, skewness and kurtosis. Eye-balling your histogram, my guess is that the mean and variance are finite, but there are methods for examining whether your distribution is fat-tailed. I believe videos 11-14 cover tails and their estimation in this (fairly good in my opinion) video series on quantitative risk management. Which moments you use, and how you use them, will depend on what you are trying to accomplish with the analysis.
A: Treating time as continuous, candidate distributions that come to mind are log-normal, gamma, and Weibull.  You can use a time-to-event regression package to fit various models and compare the AIC for which provides the best fit.  You can plot the Kaplan-Meier estimator of the survival function and overlay the fitted parametric curve.  If there is no censoring you could plot a histogram and overlay the fitted density.  It looks like there is a spike at time=0.  You could model only the non-zero times or try a finite mixture model that includes 0.
A: It is usually a bad idea to fit a probability density function, bot in the least because the shape and the results would be very sensitive to the choice of the bin size.
Robust and simple options for comparing empirical distribution with parametric ones include QQ-plot (mostly to get the intuitive idea) and the Kolmiogorov-Smirnov test (both based on comparing cdf rather than the pdf.)
A: From a practical point of view, whether a proposed calculation or test or assumption makes sense depends on the purpose. What are you going to do with the distribution once you have fitted it? If you want to use the distribution as input to simulation, then you could just use the empirical model. If you want to input the distribution to a queueing-theory model then some thought has to be given to the actions that might be taken based on the results. Some queueing-theory results are not particularly sensitive to service-time distribution, and call-center planning has to contend with forecasting error and agent attendance variability as well. If you want to apply queueing-theory, as opposed to simulation, then the question is not "what is the best-fit distribution", it is "what is the best fit of the distributions that are mathematically tractable", which may steer you to Erlang, mixed-Erland, Gamma, or LogNormal. Gamma-family is attractive mathematically, lognormal is often a good fit. Plus, "how much does the fitted or assumed distribution affect resulting management decisions?".
