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}.$$
 A: 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:

A: Firstly, without imposing a parameter range for $a$ we are obviously going to get a non-unique MLE, since we can shift it by increments of $2 \pi$ without any change to the distribution (i.e., the parameter is non-identifiable within these equivalence classes).  Given observations $x_1,...,x_n \sim \text{IID GilbertSine}(a,b)$ you can write the resulting likelihood function as:
$$\begin{align}
L_\mathbf{x}(a,b)
&= \prod_{i=1}^n \frac{\pi}{2b} \cdot \sin \bigg( \pi \cdot \frac{x_i-a}{b} \bigg) \cdot \mathbb{I}(a \leqslant x_i \leqslant a+b) \\[6pt]
&= \mathbb{I}(a \leqslant x_{(1)}) \cdot \mathbb{I}(a+b \geqslant x_{(n)}) \cdot \frac{\pi^n}{2^n b^n} \prod_{i=1}^n \sin \bigg( \pi \cdot \frac{x_i-a}{b} \bigg). \\[6pt]
\end{align}$$
Consequently, over the support of this distribution (i.e., when $a < x_{(1)}$ and $a+b > x_{(n)}$) the log-likelihood function is:
$$\ell_\mathbf{x}(a,b)
= \text{const} - n\log(b) + \sum_{i=1}^n \log \sin \bigg( \pi \cdot \frac{x_i-a}{b} \bigg),$$
and the score function is:
$$\begin{align}
s_\mathbf{x}(a,b) \equiv \nabla \ell_\mathbf{x}(a,b)
&= - \frac{\pi}{b} \begin{bmatrix} 
\sum_i \frac{\cos(\pi(x_i-a)/b)}{\sin(\pi(x_i-a)/b)} \\
\frac{n}{\pi} + \sum_i \frac{x_i-a}{b} \frac{\cos(\pi(x_i-a)/b)}{\sin(\pi(x_i-a)/b)} \\
\end{bmatrix}. \\[6pt]
\end{align}$$
Any MLE $(\hat{a},\hat{b})$ will occur at a critical point satisfying the nonlinear equations:
$$\sum_{i=1}^n \frac{\cos(\pi(x_i-a)/b)}{\sin(\pi(x_i-a)/b)} = 0
\quad \quad \quad \quad \quad 
\sum_{i=1}^n \frac{x_i-a}{b} \cdot \frac{\cos(\pi(x_i-a)/b)}{\sin(\pi(x_i-a)/b)} = - \frac{n}{\pi}.$$
You would need to use numeric methods to solve this equation to obtain the MLE.
