# How could one estimate $Y = X\beta \times \epsilon$ with $\epsilon \sim \text{Beta}(a,b)$?

Suppose the data generating process is $$Y_i=f(X_i)M_i.$$ where $$M_i$$~$$Beta(a,b)$$ and $$f$$ linear, for example,$$f(X_i)=\alpha_0+\alpha_1X_i$$. This is equivalent to $$M_i=Y_i/f(X_i).$$ I'm trying to think of someway to estimate $$a,b$$ and the $$\alpha'$$s. Do you have any suggestion?

I don't know if I have to go through Bayesian estimation.

• I imagine you could find the maximum likelihood by a computational method, where you can use as a starting point an initial least-squares fit and the moments of the residuals to come up with initial estimates. – Sextus Empiricus Sep 26 '19 at 0:20
• @MartijnWeterings It is Beta distributed. I will remove this comment. You suggested me to estimate $f(X_i)$ with some method and then use it to estimate $a,b$? I just read the title you put. It is exactly what I wanted to ask. – Celine Harumi Sep 26 '19 at 1:47
• I suggest using maximum likelihood to compute all four $a$, $b$, $\alpha_0$ and $\alpha_i$, but use a linear model to obtain initial estimates of $a$, $b$, $\alpha_0$ and $\alpha_i$ (where $a$ and $b$ are to be estimated based on the higher moments of $Y_i/\hat{f}(X_i)$). This linear model can be a regular least squares fit or if you like (not sure it is needed) something where you take heteroscedasticity into account (but I do not know your typical range of $X$, $\alpha_0$, $\alpha_i$ in order to say how you could do this, a GLM with a quasi link function to scale variance might work). – Sextus Empiricus Sep 26 '19 at 1:54
• You seem to be confident on the maximum likelihood. Could you please give me a more detailed answer? I still don't get how exactly I could calculate the likelihood function. Can you describe what you are suggesting as an answer below? – Celine Harumi Sep 26 '19 at 2:20
• I have added a code which works quite well. Note by the way that it is not exactly the same as beta-regression. However it is related to it. Beta-regression is, I believe, solved by using the maximum likelihood (at least that is how the R package betareg does it). It does not allow a simplification like solving it with an iterative reweighted least squares as in GLM (see stats.stackexchange.com/questions/304538/…). – Sextus Empiricus Sep 26 '19 at 11:05

You might try to use find parameters by maximizing the likelihood function computationally (after finding some initial estimate for the parameters).

1. linear model estimate

2. compute parameters based on linear estimates and residuals

In the code below I do this different and I assume $$a=b$$ and make use of $$Y/f \sim Beta$$. This has a problem that these values may not need to be limited in the range $$(0,1)$$. You could advance this step by making computations of the moments of the residuals and relate those to $$a$$ and $$b$$ (You should do this in order to make the convergence in the next step work better. I didn't do it because it involves annoying equations.)

3. perform computational optimization.

The code below shows a rough start for how this could work.

There are several issues:

• Values $$Y_i/f(X_i)$$ may reach outside the range $$(0,1)$$.
• The algorithm may not always converge well (With adaptations you can make it work well. The code before the edit used scaling of parameters. But more stable is the code after the edit, which uses segregated steps and solves the the problem as two nested optimizing functions).
• The error in the parameters is correlated. You can simultaneously increase the linear part $$f(X)$$ and decrease the mean of $$M_i$$, and this will give a roughly similar result. So while you can approximate the distribution well, the parameters may not be very meaningful and have a large error. (although when you test the model a few times the result does seem to turn out very close each time)

Quantile regression might be another angle to tackle this problem.

Furthermore, you can use various rough estimates instead of this computation. For instance a simple GLM model, not with a beta distribution but with some other distribution, might be sufficient for you to estimate the linear function $$f$$. How far these kind of approximations/simplifications can be made and whether they are sufficient for you, depends on the details of your problem.

example results:

True model
a = 2     b = 1      alpha_0 = 10     alpha_1 = 1

Fit
a = 1.82  b = 0.93   alpha_0 = 10.02  alpha_1 = 1.00 example code:

library(betareg)
set.seed(1)

# create modelled data
n <- 10^3
x <- runif(n,0,100)
m <- rbeta(n,2,1)
X <- cbind(rep(1,n),x,x^2)
f <- X %*% c(10,1,-0.005)
y <- f*m
colnames(X) <- c("1","x","x^2")

# likelihood to optimize
loglik <- function(par,dat=dat) {
# linear model pars
#bt0 <- par
#bt1 <- par

# ratio Y/f(X) which should relate to beta distributed var
f <- dat[,-1] %*% par
yf <- dat[,1] / f

#to prevent values outside (0,1) or values condensing in one point
neg_penalty <- which(abs(yf-0.5) < 0.4999)

# fit beta distribution to Y/f(x)
# note under the hood betareg is a call to optim
modfit <- betareg(yf[neg_penalty] ~ 1, link ="log")
a <-    exp(modfit$$coefficients$$mean)  * modfit$$coefficients$$precision
b <- (1-exp(modfit$$coefficients$$mean)) * modfit$$coefficients$$precision

# return loglikelihood or penalty
penalty_size <- length(dat[-neg_penalty,1])
if (penalty_size == 0) {
result <- -modfit$loglik + sum(log(f)) # the sum(log(f)) term is because we actually do not wish the # probability for yf but instead the probability for y which relates # to a scaled beta distribution and this scaling factor occurs in the density as log(f) } else { result <- penalty_size*10^6 } result } # data dat <- cbind(y,X) # start condition mod <- lm(y ~ 0+X) par <- c(2*mod$coefficients)  # we assume that the starting line is twice the mean

# optimize
p <- optim(par, loglik, dat=dat, control = list(trace=2, maxit=10^3))
p

# outcome
p$$par loglik(par,dat) loglik(p$$par,dat)

yf <- dat[,1] / (X %*% p$$par) modfit <- betareg(yf ~ 1, link ="log") parfin <- c(exp(modfit$$coefficients$$mean)*modfit$$coefficients$$precision, (1-exp(modfit$$coefficients$$mean))*modfit$$coefficients$$precision, p$$par)

# view result
plot(x, y,
pch=21, col=rgb(0,0,0,0.1),bg=rgb(0,0,0,0.1))

sig <- 0.02
brks <- seq(0,1,sig)
hist(yf, breaks = brks,sig,xlim=c(min(hf),max(hf)),xlab = "Y/f(X)",main="")
lines(brks,dbeta(brks, parfin,parfin)*n*sig,col=2)

parfin

• Thanks! I have to say that your estimates are quite good in your example. I need some time to understand precisely what you did. Specially, I need to learn how this function "optim" works. Again, thank you very much. – Celine Harumi Sep 26 '19 at 12:45
• optim is an algorithm that finds the coefficients that optimize (maximizes or minimizes) the function. It can use different iterative methods (iterative : making step-wise improvements of the solution). – Sextus Empiricus Sep 26 '19 at 12:51
• I see that you calculate $a$ and $b$ by taking the exponential of the mean estimate and multiply by the precision estimate. Why is that? As far as I understand, you are modelling $g(\mu)=log(\beta_0)$ and $h(\phi)=\gamma_0$ where $yf$~$Beta(\mu,\phi)$, and then construct a log-likelihood based on this. What are exactly the values that betareg estimates? – Celine Harumi Sep 27 '19 at 2:25
• I take the exponential because the betareg uses a log link function $\mu = exp (\eta)$. The precision parameter corresponds to $\phi =a+b$. Then $\mu = a/(a+b)$ turns into $$exp(\eta) = a/\phi$$ – Sextus Empiricus Sep 27 '19 at 5:41
• My reasoning is that when you scale a variable then $$f_{aX}(X) = \frac{1}{a} f_{X}(X/a)$$ So one can compute the likelihood for $f_X$ and then multiply/divide with this factor (or add/subtract when working with logarithms) – Sextus Empiricus Sep 27 '19 at 14:06