How can I test if my observed PDF follows a binomial distribution? 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.

 A: 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)


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:

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):

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)


