What probability distribution best describes my data? I have a variable whose values are summated scale scores constrained within the range 3 to 21, inclusive. I considered this variable as a scale variable (though with restricted range). In order to model the relation of this variable to a number of explanatory variables, I first wanted to know what the distribution of this variable is. For this purpose I used Cox’s dpplot command in Stata to fit various distributions to the variable. The plots are given below.

Based on the plot, what may I conclude about the distribution of the variable? It seems to be close to all distributions except to  the exponential distribution.
 A: Presumably by 'restricted range' you mean that there's both an upper and lower bound to the possible values the data can take.
Several of the distributions are reasonably consistent with the general shape of your data, but since your data are (i) discrete and (ii) bounded above, your data cannot actually come from any of those distributions.
As a way of choosing a distributional model, this activity strikes me as potentially a form of data dredging.

Note that if you're looking to use regression, then your immediate assumption is that you have not one distribution, but a different distribution (at least in respect of location) at each set of $x$ values.
As such, looking at the marginal distribution of $y$ (in terms of trying to identify a single distribution to describe it) is of little use - it doesn't relate to the regression assumptions, which involve the conditional distribution, not the marginal.
Even when the conditional distributions are not normal, you may still be able to use least squares - most forms of inference will still be okay in large samples, and even in small samples it may be okay as long as you adapt any inferential procedures. Of greater concern than distribution shape will be the assumption of linearity and homoskedasticity.
A: Let $X$ denote your variable. Assuming that the scoring system is such that 3 is the smallest possible value that could be observed, $X-3$ has integer values between 0 and 18. If the theoretical range of values is different, than you might have to do a different adjustment (or no adjustment at all). I would explore overdispersed generalizations of the binomial distribution such as the beta-binomial distribution to model your data. Depending on what you want to do, it is also possible that you don't need a parametric distribution.
