Return to Answer

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.

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.

library(betareg)
set.seed(1)

# create modelled data
n <- 10^3
x <- runif(n,0,100)
m <- rbeta(n,2,1)
f <- 10+x
y <- f*m

# likelihood to optimize
loglik <- function(par,dat=dat) {
  # linear model pars
  bt0 <- par[1]
  bt1 <- par[2]
  
  # ratio Y/f(X) which should relate to beta distributed var
  f <- bt0 + bt1*dat[,2] 
  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 ~ 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=100)) 
p

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

yf <- dat[,1] / (p$par[1] + p$par[2]*dat[,2])
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[1],
         p$par[2])

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


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

parfin

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[1]
  #bt1 <- par[2]
  
  # 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[1],parfin[2])*n*sig,col=2)

parfin

library(betareg)
set.seed(1)

# create modelled data
n <- 10^3
x <- runif(n,0,100)
m <- rbeta(n,2,1)
f <- 10+x
y <- f*m

# likelihood to optimize
loglik <- function(par,dat=dat) {
  # linear model pars
  bt0 <- par[1]
  bt1 <- par[2]
  
  # ratio Y/f(X) which should relate to beta distributed var
  f <- bt0 + bt1*dat[,2] 
  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
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=100)) 
p

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

yf <- dat[,1] / (p$par[1] + p$par[2]*dat[,2])
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[1],
         p$par[2])

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


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

parfin

library(betareg)
set.seed(1)

# create modelled data
n <- 10^3
x <- runif(n,0,100)
m <- rbeta(n,2,1)
f <- 10+x
y <- f*m

# likelihood to optimize
loglik <- function(par,dat=dat) {
  # linear model pars
  bt0 <- par[1]
  bt1 <- par[2]
  
  # ratio Y/f(X) which should relate to beta distributed var
  f <- bt0 + bt1*dat[,2] 
  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 ~ 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=100)) 
p

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

yf <- dat[,1] / (p$par[1] + p$par[2]*dat[,2])
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[1],
         p$par[2])

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


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

parfin