# fitting a distribution to skewed data with negative values

I am trying to model data about altruistic behavior in a simple lab experiment. I have one value for each participant in the sample (N=479), describing how altruistic that person was. As you can see in the graph below, the data are positively skewed, with most people having an altruism score of 0. However, there are also a few people with a negative altruism score - these "spiteful" people are willing to pay a cost to hurt others.

For simulation purposes, I am trying to fit this data to a parametric distribution. However, most distributions that can be used to model positively-skewed data, like gamma or lognormal, cannot take negative values.

One solution to this problem would be to transform my data by adding a constant so that all negative values become positive, and then fit the data to, e.g., a gamma distribution. However, this does not feel entirely satisfactory: while a gamma distribution has a strict lower bound (at 0), the variable I am trying to model doesn't (there is no hard theoretical limit on how spiteful someone could be).

Is there a probability distribution that allows for skew and also can take negative values?

edit: here is the data

  0.00  0.00  0.15  0.74  0.15  0.00  0.35  0.00  0.55  0.00  0.15  0.15 1.55  0.15  0.55  0.00  0.15  0.00 0.15 -0.45  0.15  0.75  0.00  0.35  0.00  0.00  0.00  0.00  0.15  0.55  0.35  0.00  0.35  0.00  0.55  0.00 0.15  1.55  0.35  1.55  0.00  0.15  0.00  0.15  1.35  0.15  0.15  0.00  0.00 -0.26  0.00  0.00  0.55  0.15 0.55  0.15  0.00 -0.45  1.35  0.55  0.00  0.35  0.15  0.35  0.95  0.15  0.00  0.00  0.15  0.00  0.15  0.00 0.00  0.00  0.75  0.00  0.00  0.00  0.00  0.15  0.55  0.74 -0.26  1.55  0.15 -0.26  0.00  0.00  0.15  0.00  0.00  0.15  0.15  0.15  0.15  0.00  0.15  0.00 -0.26  0.35  0.00  0.35  1.55  0.00  1.55  0.15  0.00  0.00  0.00  0.00  0.55  0.00  0.15 -0.45  0.94 -0.26  0.15  0.15  0.00  0.00  0.00 -0.45  0.00  0.00  0.00  0.35  1.55 -0.26  0.55  0.55  0.35  0.35  0.35  0.00  0.00  0.55  0.74  0.00  0.00  0.00  0.15  0.15  0.74  0.00  0.15  0.15  0.55  0.15  0.15  0.00  0.00  0.00  0.75  0.15  0.55  0.00  0.00  0.15  0.15  0.55  0.35  0.00  0.00  0.00  1.35  0.00  0.00  0.00  0.75  0.15  0.00 -0.26  0.00 -0.45  0.15  0.15  0.00  0.00  0.35  0.00  0.55  0.35  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.94  1.55  0.15  0.00  0.55  0.00  0.00  0.00  1.55  0.15  0.15 -0.26  0.35  0.00  0.15  0.15  0.00  0.15  0.00 -0.45  0.00  0.00  0.00  0.15  0.00  0.00  0.94  0.00  0.35  0.00  0.00  0.55  0.55  0.15  0.35  0.00  0.00 -0.26  0.00  0.15  0.00  0.15  0.00  0.35  0.00  0.15  0.00  0.15  1.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15 0.00 -0.45  0.00  0.00  0.74  0.74  0.00  0.35  0.15  0.55  0.00  0.15  0.15  0.15  0.55  1.55  0.00  0.00  0.00 -0.45  0.15  1.15  0.00  0.15  0.15  0.15  0.15  0.35  0.15  0.55  0.94  0.00  0.15  0.00  0.00  0.35  0.94  0.75  0.00  0.00  0.15  0.00  0.15  0.55  0.35  0.00  0.00  0.35  1.55 -0.26  0.00  0.00  0.21  0.00  0.00  0.15  0.00  1.15 -0.45  0.00  0.55  0.00  0.35  0.35  0.35  0.15  0.15  0.15  0.00  0.95 -0.26  0.94  0.00  1.35  0.15  0.35  0.00  0.55  0.65  0.00  1.15 -0.45  0.15  0.00  0.15  0.15  0.15  0.15  0.35  0.15  -0.26  0.15  0.35  0.35  0.00  0.94  0.00  0.00  0.00  0.15  0.15  0.74  0.00  0.94 -0.45  0.00  0.00  0.00 0.35  0.15  0.00  0.00  0.15  1.35  0.00  0.15  0.00  0.15  0.35  0.94  0.35  0.00  0.00  0.55  0.55  0.00  0.00  0.00  0.55  1.35  0.00  0.55  0.55  0.00  0.15  0.35  0.15  0.00  0.00  0.00  1.35  0.15  0.00  0.35  0.35  0.00  0.15  1.55  0.55  0.00  0.35  0.15  0.15  0.15  0.15  0.15  0.94  0.35  0.95  0.15  0.00  0.00  0.00  0.15  0.00  0.55  0.15  0.00  0.15  1.55  0.00  0.75 -0.45  0.00  0.00  0.35  0.00  0.94  0.00  0.15  0.00 -0.45  0.00  0.00  0.00  1.15  0.00  0.15  0.15 -0.26  0.35  0.35  0.35  0.15  0.15 -0.26  0.35 -0.45 0.00  0.55  0.00 -0.26  0.15  0.55 -0.26  0.00  0.15  0.15  0.94  0.94  0.15  0.00  0.00  0.00  0.35  0.55  0.15  0.00  0.00  0.00  0.15  0.00  0.00  0.00  0.00  0.55  0.00

• Of course there are. Any positive valued continuous function defined on an interval where the lower endpoint is negative is not symmetric and integrates to 1 will work. One such example is the skew normal distribution which you can look up on Wikipedia. – Michael R. Chernick Apr 5 '18 at 22:30
• Do you mind sharing the data? – Georg M. Goerg Apr 5 '18 at 22:58
• With a sample size of 479, why do you feel the need to fit a parametric probability distribution? Why not just get a nonparametric density estimate (such as what one gets with the density function in R)? This is not to suggest that you don't have a need to fit a parametric probability distribution but only that you might want to express why you think that is necessary (assuming that it is necessary). – JimB Apr 5 '18 at 23:23
• I have two reasons in mind: 1)Using a parametric distribution makes it easier to write code that generates data from this distribution (at least, using what I know in R). 2)Assuming I publish the outcome of this research, it makes it easier for readers of the paper to reproduce my results, just by using the parameters of the distribution given in the paper. – Choirofangels Apr 6 '18 at 0:05
• If the parametric distribution is a good fit, then for whatever it's worth I think your reasons are good. If you can't find one with a good fit, then consider stats.stackexchange.com/questions/82797/… to eliminate reason (1). – JimB Apr 6 '18 at 4:25

Not a complete answer, but people will see why this is not a comment. I get 473 values here and suggest that the ambition to fit a smooth distribution here needs to surmount marked granularity in the data. What protocol lies behind this? The small moral that binning can conceal fine structure as well as noise should need little emphasis. The square root scale (compare J.W. Tukey's "rootograms") is chosen as a matter of convenience.

      score |      freq.           %      cum. %
------------+-----------------------------------
-.45 |         14        2.96        2.96
-.26 |         16        3.38        6.34
0 |        186       39.32       45.67
.15 |        119       25.16       70.82
.21 |          1        0.21       71.04
.35 |         46        9.73       80.76
.55 |         36        7.61       88.37
.65 |          1        0.21       88.58
.74 |          7        1.48       90.06
.75 |          6        1.27       91.33
.94 |         13        2.75       94.08
.95 |          3        0.63       94.71
1.15 |          5        1.06       95.77
1.35 |          7        1.48       97.25
1.55 |         13        2.75      100.00
------------+-----------------------------------
total |        473      100.00


• The binning comes from the nature of the experiment (getting less coarse data would have required too many experimental trials). I understand that this level of granularity is less than ideal, but how problematic would you say it is (minor inconvenience or fatal problem)? – Choirofangels Apr 6 '18 at 0:58
• I am more puzzled by the small gaps (e.g. 0.74, 0.75; 0.94, 0.95) than by the large ones! Also, any story on why no-one at around $-$0.15 given that 0 is a strong mode? – Nick Cox Apr 6 '18 at 1:07
• Were there limits on possible scores? My guess is a rescaled beta distribution is a fair starting point. – Nick Cox Apr 6 '18 at 1:40
• There were limits on possible scores: a participant couldn't have a score lower than -.45 or higher than 1.55. However, these bounds are fairly arbitrary: one could have designed a similar experiment with more trials that could have enabled higher and lower scores. In other words, -.45 does not represent "maximum possible spitefulness" but rather the maximum amount of spitefulness that could be detected in this particular experiment. Therefore, ideally I would be looking for a distribution that has no bounds either. – Choirofangels Apr 6 '18 at 4:01
• As the limits can be attained, I think beta is off the table. I don't have further suggestions as I can't tell between respecting your data, respecting your protocol and respecting your desire to infer the distribution of a hypothetical latent variable. How can anyone guess whether the asymmetry is a fact of psychology or imparted by your procedure? – Nick Cox Apr 6 '18 at 8:22

Just some ideas, but yours seems to be a difficult case for distribution fitting. I will disregard the bounds for now, and show some example code in R.

text  <-  "0.00  0.00  0.15  0.74  0.15  0.00  0.35  0.00  0.55  0.00  0.15  0.15 1.55  0.15  0.55  0.00  0.15  0.00 0.15 -0.45  0.15  0.75  0.00  0.35  0.00  0.00  0.00  0.00  0.15  0.55  0.35  0.00  0.35  0.00  0.55  0.00 0.15  1.55  0.35  1.55  0.00  0.15  0.00  0.15  1.35  0.15  0.15  0.00  0.00 -0.26  0.00  0.00  0.55  0.15 0.55  0.15  0.00 -0.45  1.35  0.55  0.00  0.35  0.15  0.35  0.95  0.15  0.00  0.00  0.15  0.00  0.15  0.00 0.00  0.00  0.75  0.00  0.00  0.00  0.00  0.15  0.55  0.74 -0.26  1.55  0.15 -0.26  0.00  0.00  0.15  0.00  0.00  0.15  0.15  0.15  0.15  0.00  0.15  0.00 -0.26  0.35  0.00  0.35  1.55  0.00  1.55  0.15  0.00  0.00  0.00  0.00  0.55  0.00  0.15 -0.45  0.94 -0.26  0.15  0.15  0.00  0.00  0.00 -0.45  0.00  0.00  0.00  0.35  1.55 -0.26  0.55  0.55  0.35  0.35  0.35  0.00  0.00  0.55  0.74  0.00  0.00  0.00  0.15  0.15  0.74  0.00  0.15  0.15  0.55  0.15  0.15  0.00  0.00  0.00  0.75  0.15  0.55  0.00  0.00  0.15  0.15  0.55  0.35  0.00  0.00  0.00  1.35  0.00  0.00  0.00  0.75  0.15  0.00 -0.26  0.00 -0.45  0.15  0.15  0.00  0.00  0.35  0.00  0.55  0.35  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.94  1.55  0.15  0.00  0.55  0.00  0.00  0.00  1.55  0.15  0.15 -0.26  0.35  0.00  0.15  0.15  0.00  0.15  0.00 -0.45  0.00  0.00  0.00  0.15  0.00  0.00  0.94  0.00  0.35  0.00  0.00  0.55  0.55  0.15  0.35  0.00  0.00 -0.26  0.00  0.15  0.00  0.15  0.00  0.35  0.00  0.15  0.00  0.15  1.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15  0.15 0.00 -0.45  0.00  0.00  0.74  0.74  0.00  0.35  0.15  0.55  0.00  0.15  0.15  0.15  0.55  1.55  0.00  0.00  0.00 -0.45  0.15  1.15  0.00  0.15  0.15  0.15  0.15  0.35  0.15  0.55  0.94  0.00  0.15  0.00  0.00  0.35  0.94  0.75  0.00  0.00  0.15  0.00  0.15  0.55  0.35  0.00  0.00  0.35  1.55 -0.26  0.00  0.00  0.21  0.00  0.00  0.15  0.00  1.15 -0.45  0.00  0.55  0.00  0.35  0.35  0.35  0.15  0.15  0.15  0.00  0.95 -0.26  0.94  0.00  1.35  0.15  0.35  0.00  0.55  0.65  0.00  1.15 -0.45  0.15  0.00  0.15  0.15  0.15  0.15  0.35  0.15  -0.26  0.15  0.35  0.35  0.00  0.94  0.00  0.00  0.00  0.15  0.15  0.74  0.00  0.94 -0.45  0.00  0.00  0.00 0.35  0.15  0.00  0.00  0.15  1.35  0.00  0.15  0.00  0.15  0.35  0.94  0.35  0.00  0.00  0.55  0.55  0.00  0.00  0.00  0.55  1.35  0.00  0.55  0.55  0.00  0.15  0.35  0.15  0.00  0.00  0.00  1.35  0.15  0.00  0.35  0.35  0.00  0.15  1.55  0.55  0.00  0.35  0.15  0.15  0.15  0.15  0.15  0.94  0.35  0.95  0.15  0.00  0.00  0.00  0.15  0.00  0.55  0.15  0.00  0.15  1.55  0.00  0.75 -0.45  0.00  0.00  0.35  0.00  0.94  0.00  0.15  0.00 -0.45  0.00  0.00  0.00  1.15  0.00  0.15  0.15 -0.26  0.35  0.35  0.35  0.15  0.15 -0.26  0.35 -0.45 0.00  0.55  0.00 -0.26  0.15  0.55 -0.26  0.00  0.15  0.15  0.94  0.94  0.15  0.00  0.00  0.00  0.35  0.55  0.15  0.00  0.00  0.00  0.15  0.00  0.00  0.00  0.00  0.55  0.00"
dat  <-  scan(textConnection(text))


One idea is to try some Box-Cox transform, but the standard formulation cannot be used as there are some negative values.

library(MASS) # We will use this later
library(car) # We use boxCox from car
boxCox(dat ~ 1, family="yjPower",param="gamma")  # this produces a plot, indicating lambda around -0.5  (not shown)
boxCox(dat ~ 1, family="bcnPower")  # This gives somewhat conflicting advice, lambda close to zero but positive, again plot not shown


I will leave that kind of method only as an idea here. Now, for a different idea, we can try to fit a skew-normal distribution (for this data I believe more in this last idea than the former one):

library(sn)
skew_mod <-  selm(dat ~ 1) # selm is "skew-elliptic lm"
summary(skew_mod)
Call: selm(formula = dat ~ 1)
Number of observations: 473
Family: SN
Estimation method: MLE
Log-likelihood: -159.4612
Parameter type: CP

CP residuals:
Min       1Q   Median       3Q      Max
-0.69865 -0.24865 -0.09865  0.10135  1.30135

Regression coefficients
estimate  std.err  z-ratio Pr{>|z|}
mean  0.24865  0.01651 15.06170        0

Parameters of the SEC random component
estimate std.err
s.d.     0.3603   0.013
gamma1   0.7498   0.035


and finally we show a histogram with the estimated skew-normal density overlaid:

hist(dat,prob=TRUE,nclass="scott") # "scott" is from MASS
plot(function(x) dsn(x, dp=skew_mod@param$dp), from=-0.5, to=1.75, col="red", add=TRUE)  The probability mass outside the bounds is so low that that should not cause problems. For sampling one could simply reject simulated values outside bounds. But the fit is not too good: the skewness is captured, but the histogram has much more mass at the center. Maybe a skew-t distribution could give a better fit. I will leave that for the OP, but the fitting can be done with the software I have illustrated here. The other answer by Nick Cox address the obvious granularity in the data, I have ignored that aspect, but before doing anything you should ask yourself why is that. If it is important to maintain that granularity in simulations, what I have done here is not enough. Thanks to everyone who helped! I have finally found that a skewed laplace distribution provides an excellent fit. #package rmutil has the skewed laplace distribution library(rmutil) #package fitdistrplus is used to fit the data to a distribution library(fitdistrplus) #estimate which parameters of a skewed laplace distribution fit the data best fit.sklaplace <- fitdist(AltruismScore, "skewlaplace", start=list(0,.2,.5)) #this yields the following parameters: m=0.00, s=.23, f=.63. #plot the distribution thus obtained over the data range <- seq(-1,2.5,by=.01) histogram<-hist(AltruismScore, prob=TRUE, xlab="AltruismScore", main="Histogram of AltruismScore", col="lightgreen", xaxt="n") histogram lines(range, dskewlaplace(range, 00, .23, .63), col="darkblue", lwd=2) axis(side=1,at=histogram$mids)


edit: histograms with better bins

• Your histogram doesn't do you any favours here. Your data show a mode at 0 and that's not quite obvious in the histogram. Evidently your bins are of the form $a < x \le b$ so that when $b = 0$ the zeros are included in the bin $(-0.2, 0]$. I'd try bin limits that go $-0.7(0.2)1.7$ or $-0.65(0.1)1.65$ so that the bin that includes $0$ has $0$ as its midpoint. On the larger issue, I remain sceptical that an unbounded distribution is of that much use or interest for a bounded response. – Nick Cox Apr 23 '18 at 1:27
• I have uploaded a new histogram above following your suggestions. The weird thing is that the axis tick at -.4 is left out for some reason. Code: histogram<-hist(AltruismScore, prob=TRUE, xlab="AltruismScore", main="Histogram of AltruismScore", col="lightgreen", xaxt="n", breaks=seq(-.7,1.7,.2)) histogram lines(range, dskewlaplace(range, 00, .23, .63), col="darkblue", lwd=2) axis(side=1,at=round(histogram$mids, 2)) – Choirofangels Apr 23 '18 at 3:01 • I see a tick at$-0.6\$ but no text -0.6. I think R does that when the axis gets too crowded. – Nick Cox Apr 23 '18 at 7:03