How to calculate the MLE for Gilbert's Sine Distribution? $$ f(x) = \frac{\pi}{2b}\sin\left(\pi\frac{x-a}{b}\right),~~x\in[a,a+b]$$
So it's a pretty simple distribution, which is just a scaled sine wave. As an experiment, I was trying to find the MLE estimates of $a$ and $b$ but I can't find it for $b$ and the result for $a$ doesn't seem to work. Intuitively, you want to move $a$ and $b$ closer to the point $x. $
$$\frac{\partial\textrm{ MLE}}{\partial a}=-\frac{\pi^2\cos\left(\pi\frac{x-a}{b}\right)}{2b^2};$$
root/zero found at $a = \frac{2x+b}2.$
$$\frac{\partial\textrm{ MLE}}{\partial b}=-\frac{\pi\cdot \left(\sin\left(\pi\frac{x-a}{b}\right)b+ (\pi x-\pi a) \cos\left(\pi\frac{x-a}{b}\right)\right)}{2b^3}.$$

The log of the likelihood for a random sample from that distribution is
$$\text{logL}=\sum _{i=1}^n \log \left(\sin \left(\frac{\pi  (x_i-a)}{b}\right)\right)-n \log (b)-n \log (2)+n \log (\pi )$$
So we choose values of $a$ and $b$ to maximize the likelihood.  We'll need to do this in an iterative manner and need good starting values.
The mean and variance of Gilbert's Sine distribution are $a+b/2$ and $\frac{\left(\pi ^2-8\right) b^2}{4 \pi ^2}$, respectively.  The method of moments estimators are then found by equating the theoretical moments and the sample moments so that
$$a_0=\bar{x}-\frac{\pi  s}{\sqrt{\pi ^2-8}}$$
$$b_0=\frac{2 \pi  s}{\sqrt{\pi ^2-8}}$$
where $\bar{x}$ is the sample mean and $s$ is the sample standard deviation.
These two values can be used for starting values in the iterative procedure that finds the maximum likelihood estimates.  (Probably a safer set of starting values would be $\min\left(\min(x),a_0\right)$ and $\max\left(\max(x)-\min(x),b_0\right)$, respectively, but I'll ignore that for now.)
Some R code to generate a random sample and then find the maximum likelihood estimates of $a$ and $b$.
# Set values for a and b
  a <- 3
  b <- 2
  
# Sample size
  n <- 100
  
# Generate a random sample
  set.seed(12345)
  u <- runif(n)
  x <- (a*pi + 2*b*asin(sqrt(u)))/pi
  
# Define a function for the log of the likelihood
  logL <- function(parms, x) {
    a <- parms[1]
    b <- parms[2]
    n <- length(x)
    # Check for legitimate parameter values 
    if (a < min(x) & (a+b) > max(x)) {
       -n*log(2) - n*log(b) + n*log(pi) + sum(log(sin(pi*(-a + x)/b)))
    } else {
       -10^200
    }
  }
  
# Find maximum of log likelihood
  # Starting values
  a0 <- mean(x) - sd(x)*pi/sqrt(pi^2-8)
  b0 <- 2*sd(x)*pi/sqrt(pi^2-8)
  
  results <- optim(c(a0, b0), logL, x=x, method="L-BFGS-B",
    lower=c(min(x)-2*(max(x)-min(x)), max(x)-min(x)), 
    upper=c(min(x), 2*(max(x)-min(x))),
    control=list(fnscale=-1), hessian=TRUE)

# Covariance matrix
  covmat <- -solve(results$hessian)

# Show results
  cat("Parameter\tEstimate\tStd. Err.\n",
      "a\t\t", results$par[1], " \t", covmat[1,1]^0.5, "\n",
      "b\t\t", results$par[2], " \t", covmat[2,2]^0.5, "\n", sep="")
#Parameter       Estimate        Std. Err.
#a               2.950283        0.0577424
#b               2.120144        0.1106378

Whatever optimization routine you use, good starting values are your friends.  The method-of-moments estimator should work most of the time.  If you use optim in R and the "L-BFGS-B" method, this restricts the possible parameter values to a rectangular region but the values that are legitimate are in a "triangular" region so one could possibly end up with estimates that aren't supported by the data and probability model.  Here's how that triangular region is constructed:

