How do I use a beta distribution to model this bimodal data? I am trying to model data that looks like this:

The data reflects frequencies of items at 0, 25, 50, 75, and 100. The majority of the time, values are concentrated at the 0, and 100 points. I believe a beta distribution would be the best model for this. I am having trouble trying to calculate the alpha and beta values to fit this data. Any advice on how to do this would be appreciated.
Edit: For some more context, the data points of 0, 25, 50, 75, 100, reflect continuous data points. To give some more perspective, let's say those values reflect percentages, 0% to 100%. We can assign it to an online course, where many students may sign-up, but never touch the course, resulting in many 0% users, and many end up finishing the course as well (100%).
 A: If the only values the data take on are $\{0,25,50,75,100\}$ then the Beta isn't an appropriate model since it's continuous and the data are highly discrete. If this is the case, and you don't want to make assumptions about the distribution, then just use the empirical distribution (i.e. if the data are $X_1,\ldots,X_n$, then $\hat{p}(k) = \frac{1}{n} \sum_{i=1}^n \mathbb{1}(X_i = k)$) and a bootstrap to assess confidence intervals/sets since your sample size is reasonable but not huge.
If the data take on continuous values from 0 to 100, or integers $0,\ldots,100$, then a Beta may be a good model. Divide all the data by 100 so that they lie in $[0,1]$, and then find the Beta MLE (this R package will do it for you, for example). In the analysis after this, just multiply predictions, means, standard errors, etc. by 100 to get back to the original data scale.
A: I do not disagree with @SheridanGrant's Answer. However, I will directly address your question about using a beta distribution. Even if you are only observing values
"rounded" (or otherwise transformed) to $0, .25, .50, .75, 1,$ there may be a beta distribution determining the probabilities with which these values are observed.
Following up on @BenBolker's Comment about estimating
beta parameters by the method of moments, I digitized (approximately) the counts in the histograms you
provided, to obtain $\mu \approx 0.5$ and $\sigma^2 \approx 0.2.$ 
Because $\mu = \alpha/(\alpha + \beta),$ we observe that
$\mu = 0.4$ imples $\alpha  = \beta.$ Then because
$\sigma^2 =  \frac{\alpha\beta}{(\alpha + \beta)^2(
\alpha+\beta+1)}$ we can set $\alpha = \beta$ to get $\alpha=\beta=1/8).$ 
We can quickly verify this in R by
looking at a million observations from $\mathsf{Beta}(1/8, 1/8),$ to get $\mu \approx 1/2, \sigma^2 \approx 0.2.$ (With a million observations, we can expect 2 or 3-place accuracy.) Even if I didn't digitize your 'data' correctly,
the same method would work with the correct data.
x = rbeta(10^6,1/8,1/8)
mean(x); var(x)
[1] 0.5003867
[1] 0.1999466

Then looking at a smaller sample $(n = 500)$ from this beta distribution we get a histogram that is not much different than the ones you show.
set.seed(2019)  # for reproducibility
x = rbeta(500, 1/8, 1/8)
hist(x, br=5, prob=T, ylim=c(0,3), col="skyblue2", label=T)
  curve(dbeta(x,1/8,1/8), add=T, col="red", n=10001)


Because the width of each bin in this histogram is 0.2, we can multiply the 'density' labels atop the histogram bars by 0.2 to get the probability represented by each bar (for a total of $1).$
Theoretical probabilities (to 4 places) in these five intervals are as follows:
round(diff(pbeta(c(0,.2,.4,.6,.8,1), 1/8, 1/8)), 4)
[1] 0.4269 0.0513 0.0435 0.0513 0.4269

With sample sizes of only $n = 100$ observations, results are more variable as illustrated by the four histograms below.

Note: The density function of $\mathsf{Beta}(1/8, 1/8)$ has vertical asymptotes at $0$ and $1.$ So there is more area under the density curve near $0$ and $1$ than may be apparent in the plots above.
