5
$\begingroup$

I have got an observed PDF and I would like to test if it fits well a binomial distribution, i.e. if I can affirm that my observed data follow a binomial distribution.

I am using R, is there any package that can help me doing so?

The picture shows my observed PDF in bold black and a theoretical binomial distribution.

enter image description here

$\endgroup$
6
  • 2
    $\begingroup$ Do a chi-squared test. $\endgroup$ May 8, 2017 at 18:59
  • 1
    $\begingroup$ I keep an anonymous site with notes, and I wonder if example 7 in the link fits your question. Let me know... $\endgroup$ May 8, 2017 at 19:26
  • 1
    $\begingroup$ You can also use the R package fitdistrplus as in here. $\endgroup$ May 8, 2017 at 19:40
  • $\begingroup$ @David A chi-squared test will have low power against smooth deviations compared to other alternatives. $\endgroup$
    – Glen_b
    May 9, 2017 at 0:06
  • 1
    $\begingroup$ Possible duplicate of Get p-value from fidistrplus $\endgroup$ May 9, 2017 at 18:09

1 Answer 1

8
$\begingroup$

Please note that this question is duplicated, and there is a better answer by Glen_b over here.


Here is an example using the the R package fitdistrplus. The simulated data is strictly binomial, and hence we should fail to reject the null hypothesis.

require("fitdistrplus")
set.seed(10)
n = 25
size = 27
prob = .4
data = rbinom(n, size = size, prob = prob)
fit = fitdist(data = data, dist="binom", 
                   fix.arg=list(size = size), 
                   start=list(prob = 0.1))

summary(fit)
Fitting of the distribution ' binom ' by maximum likelihood 
Parameters : 
      estimate Std. Error
prob 0.3822225 0.01870338
Fixed parameters:
     value
size    27
Loglikelihood:  -52.24948   AIC:  106.499   BIC:  107.7178 

plot(fit)

enter image description here

This estimation of the probability of success ($\Pr(S) = 0.38$) is very close to $\Pr(S) =0.4$ set up when the data was generated. However, you are not just estimating the parameter, but also the distribution, as Glen_b explains on the parallel post:

... not only are you estimating parameters (which may be able to be dealt with - some tests adapt to that more or less easily) but frequently you're also selecting between several choices of distributional model.

You may want, therefore, to contrast the result with other distributions. For example here is the output selecting a Poisson distribution:

fit2 = fitdist(data, dist = "pois")
summary(fit2)
plot(fit2)
Fitting of the distribution ' pois ' by maximum likelihood 
Parameters : 
       estimate Std. Error
lambda    10.32  0.6424951
Loglikelihood:  -55.95267   AIC:  113.9053   BIC:  115.1242 

Predictably, the AIC increases: we have set up the data as binomial, so it would be expected that the better fitting distribution (lower AIC) is binomial, and not Poisson. Here are the corresponding plots:

enter image description here

Notice that we have to supply the probability of success to then estimate the goodness of fit. In the simulation above, we asked fitdistrplus() to start estimating MLE from a value of $\small \Pr(\text{success})=0.1$.

If we were to run a chi square goodness of fit test we could try several values. Here is a simulation with the same dataset as above:

par(mfrow = c(1,3))
obs.counts  = table(factor(data, levels = 0:27)) # Tabulated results of dataset
plot(obs.counts,"h",  col = "turquoise", main = "Obs'd counts", cex.main = .9)
values = seq(0.3, 0.7, 0.1) # P(S) chosen (.3,.4,.5,.6,.7)
targ.freq = matrix(0, length(values), size + 1)  # Starting an empty matrix
# PMF for different p (Prob) values:
for(i in 1:length(values)) targ.freq[i,] = dbinom(0:size, size = size, prob = values[i]) 
plot(targ.freq[1,], type = "l", col = "darkblue", 
main = "Exp'd freq", cex.main = .9, ylab = 'Freq', xlab="")
for(i in 1:nrow(targ.freq)) lines(targ.freq[i,], col = "darkblue")
for(i in 1:nrow(targ.freq)) text(which.max(targ.freq[i,]), 
                                 max(targ.freq[i,] + .005), values[i],
                                 col= "darkblue")
exp.counts = matrix(0, length(values), size + 1) # Starting empty matrix
# Calculating expected counts...
for(i in 1:nrow(targ.freq)) exp.counts[i,]  = round(targ.freq[i,] * length(data)) 
colnames(exp.counts) = 0:size
vec = c(1,3,5)
matplot(t(exp.counts)[ ,vec],type="h",col=c(2:nrow(exp.counts)+1),lty=1,
        ylab ="Exp. counts", main = "Exp counts 3 Pr(S)")
lab = values[c(1,3,5)]
x.positions = c(7,15,22)
for(i in 1:nrow(targ.freq)) text(x.positions[i], 
                            max(exp.counts[i,] + .1), lab[i], col= "darkblue")

The expected frequencies and counts (and consequentially the fit to the observed data) will change with the different probabilities of success (on the right plot only three of them were plotted):

enter image description here

Likewise, the p values of the GOF chi square will show a peak around the probability we set up in the generation of the data:

p.values = NULL
for(i in 1:nrow(targ.freq)) p.values[i] = chisq.test(obs.counts , p = targ.freq[i,])$p.value
par(mfrow = c(1,1))
plot(values, p.values, type = "h", main = "Chi square p values", cex.main = .8)

enter image description here

$\endgroup$
3
  • $\begingroup$ thank you very much. regarding the fitdist example, why you set start=list(prob = 0.1) ? $\endgroup$
    – aaaaa
    May 9, 2017 at 7:36
  • $\begingroup$ in other words, I don't understand the start argument and eventually the relative probability that needs to be specified. $\endgroup$
    – aaaaa
    May 9, 2017 at 8:12
  • $\begingroup$ and also, how can I get the p-value from fitdist? thanks $\endgroup$
    – aaaaa
    May 9, 2017 at 8:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.