4

I am trying to fit a beta distribution to election forecast data. The ultimate purpose is determining with what probability the election will be decided by one vote (more on this here).

My data is as follows:

0.517 0.59 0.55 0.528 0.496 0.51 0.55 0.54 0.57

I use the scipy scipy.stats.beta.fit() function, as discussed here.

For my data I get:

$\alpha = 0.7469674643946238 $

$\beta = 1.7165349362225657 $

which leads to a sloping distribution with its maximum close to zero.

For most other artificial data I have tried, I get distributions which actually hold water. Obviously, for my case this is wrong (it may not be for others), but I would like to make certain that the minima of the distribution are close to 0 and 1 respectively, and the maximum is approximately in the range of my data points.

I had thought of fitting a normal distribution to my data, and a beta distribution to that - but that would be an approximation of an approximation. I'd rather do it all in one "step".

Community
  • 1
TheChymera
  • 754
  • 2
  • 10
  • 24
  • 4
    16 decimal places for the parameter estimates when the data are to 2 d.p.! – Nick Cox Nov 16 '14 at 16:44
  • I don't get your point. The distribution will still look the same if I round the parameters to fewer decimals. – TheChymera Nov 16 '14 at 16:48
  • 2
    That was just a quick initial reaction, but spurious precision often is a bad sign and encouraging better practices in reporting is surely consistent with the aims here. In terms of your question, something appears very wrong with your estimates, which imply a mean of about 0.303, which isn't even within the range of your data. At this point, the question is really whether you found a bug in the software or misapplied the software, and is possibly outside the scope of the forum. I see absolutely no point to fitting a normal first. – Nick Cox Nov 16 '14 at 16:54

1 Answers1

6

It looks like the scipy.stats.beta.fit function is weird and undocumented. According to this question, what it does is fit a generalized beta distribution (i.e. a beta distribuction which has been shifted and stretched to be nonzero on an interval other than $[0,1]$) with max and min included as parameters. If you don't specify the max and min, then it will choose weird values instead of the standard zero and 1 (see the top answer to the linked question.) So it looks like you should try

scipy.stats.beta.fit(data, floc=0, fscale=1)

on your data. Unfortunately, I don't have access to it right now so can't try it, but the correct values should be $\alpha = 168$, $\beta = 144$. You can get them from an online app here.

Community
  • 1
Flounderer
  • 9,575
  • 1
  • 32
  • 43
  • I actually discovered this independently in the mean time, but it's encouraging to also hear it from you, I was afraid I might just have "hacked" it. – TheChymera Nov 16 '14 at 19:59
  • And with such large values of the parameters, the normal approximation is essentially indistinguishable. – Alecos Papadopoulos Nov 16 '14 at 20:00
  • 2
    From a variety of starting points, R's `fitdistr` gives values for $\alpha$ between about 166 and 173 and for $\beta$ between about 142 and 148 (in each case the standard error is a bit under half the estimate). It looks like there's a reasonably broad range of values in the region of the values you give, with similar fits. The two estimates are very highly correlated. – Glen_b Nov 16 '14 at 22:21
  • 1
    Regarding the high correlation - their ratio seems to be pretty well determined (e.g. we can estimate the mean accurately), but their sum is not (we don't have such a good estimate of the variance). – Glen_b Nov 16 '14 at 22:27
  • so you're saying `scipy.stats.beta.fit` is noticeably not as good as R's `fitdistr`? Is that also fitting a beta distribution? – TheChymera Nov 17 '14 at 00:43
  • 1
    Stata program `betafit` from SSC with defaults gives $\alpha = 169, \beta = 145$, rounding to nearest integer, and so matching @Glen_b's results. – Nick Cox Nov 17 '14 at 19:59
  • 1
    TheChymera -- I don't see how you'd get that from what I said. I was confirming Flounderer's results, but indicating that other start values might give what read like different answers (as would different optimizers), but those fits may be almost identical. R's `fitdistr` will fit whatever you tell it to (the beta is one of the built-in densities, but you can supply a density function of your choice). I asked it to fit the usual two-parameter beta. – Glen_b Nov 17 '14 at 21:33