ML estimation not working for some cases of the Beta distribution

[wleoncio]

See example below:

[wleoncio]

This might be caused by a runaway delta, see example below:

r$> sample.beta <- rtruncbeta(1000, shape1 = 15, shape2 = 4, a = .7, b = .9)

r$> ml_beta <- mlEstimationTruncDist(
      sample.beta, print.iter = TRUE, tol = 1e-7, max.it = 1e3
    )
Estimating parameters for the beta distribution
it:  1 delta:  24.87277  - parm:  43.313 10.702 
it:  2 delta:  61.23992  - parm:  48.136 11.969 
it:  3 delta:  200.9796  - parm:  55.718 13.908 
it:  4 delta:  1332.409  - parm:  69.465 17.373 
it:  5 delta:  376016.3  - parm:  104.87 26.254 
it:  6 delta:  NaN  - parm:  -489.678 -123.842 
Error in while ((delta.L2 > tol) & (it < max.it)) { : 
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In dbeta(y, shape1 = parm[1], shape2 = parm[2]) : NaNs produced
2: In pbeta(a, shape1 = parm[1], shape2 = parm[2]) : NaNs produced
3: In pbeta(b, shape1 = parm[1], shape2 = parm[2]) : NaNs produced

I'll write Rene on this, but I found this function to be strange:

natural2parameters.trunc_beta <- function(eta) {
  # eta: The natural parameters in a beta distribution
  # returns (alpha,beta)
  parms <- c(shape1 = eta[1], shape2 = eta[2])
  class(parms) <- class(eta)
  return(parms)
}

Usually parms does a transformation on eta. Here it just splits it. I suspect this could be causing the issue (not the case, see comment below).

[wleoncio]

There's a nice table on https://en.wikipedia.org/wiki/Exponential_family#Table_of_distributions that suggests two variants of the Beta: one has no transformation and the other has a ±1 transformation between natural and inverted mapping. Neither solution solves this problem, so it might be somewhere else.

[wleoncio]

Hi @rho62,

Is there a reference to the mothodology behind the getYseq() and getGradETinv() functions I could read? Trying to check if the following code chunk is correct:

TruncExpFam/R/beta.R

Lines 77 to 95 in 4ec917e

	getYseq.trunc_beta <- function(y, y.min = 0, y.max = 1, n = 100) {
	# needs chekking
	mean <- mean(y, na.rm = TRUE)
	sd <- var(y, na.rm = TRUE)^0.5
	lo <- max(y.min, mean - 5 * sd)
	hi <- min(y.max, mean + 5 * sd)
	out <- seq(lo, hi, length = n)
	class(out) <- class(y)
	return(out)
	}
	getGradETinv.trunc_beta <- function(eta) {
	# eta: Natural parameter
	# return the inverse of E.T differentiated with respect to eta' : p x p matrix
	term.1 <- sum(1 / (((1:10000) + eta[1]))^2)
	term.2 <- sum(1 / (((1:10000) + eta[2]))^2)
	term.12 <- sum(1 / (((1:10000) + eta[1] + eta[2]))^2)
	return(A = solve(matrix(c(term.1 - term.12, -term.12, -term.12, term.2 - term.12), ncol = 2)))
	}

[wleoncio]

Hi René,

Great news regarding the beta estimation.

I didn't manage to confirm the previous implementation, so I recoded the whole function following this procedure:

1. Calculated the dE(T)/d(eta) matrix by hand:
   

2. Got the moments from here and used the approximation for the digamma described here

3. Calculated the derivatives using Wolfram Alpha.

This seems to solve this issue, the last of those critical ML problems we found. The issues list currently contain only one small annoying bug regarding parameter naming (#74), and the rest are new/improved features. Therefore, I think we should release version 1.0.1 with the current fixes before starting work on the next thing, what do you say?

[rho62]

Smart! Very glad to hear this! 😃 It sounds like a good idea to release v 1.0.1, but I’d like to have a better overview first. Can we touch upon this Friday, too? /René Fra: Waldir Leoncio ***@***.***> Svar til: ocbe-uio/TruncExpFam ***@***.***> Dato: tirsdag 23. august 2022 kl. 10:16 Til: ocbe-uio/TruncExpFam ***@***.***> Kopi: Rene Holst ***@***.***>, Mention ***@***.***> Emne: Re: [ocbe-uio/TruncExpFam] ML estimation not working for some cases of the Beta distribution (Issue #85) Hi René, Great news regarding the beta estimation. I didn't manage to confirm the previous implementation, so I recoded the whole function following this procedure: 1. Calculated the dE(T)/d(eta) matrix by hand: [Bilde er fjernet av sender. image]<https://user-images.githubusercontent.com/8234768/186106651-ed4668d4-f512-4b7d-98d6-380b5d9a94b9.png> 2. Got the moments from here<https://en.wikipedia.org/wiki/Beta_distribution#Moments_of_logarithmically_transformed_random_variables> and used the approximation for the digamma described here<https://en.wikipedia.org/wiki/Digamma_function> 3. Calculated the derivatives using Wolfram Alpha<https://www.wolframalpha.com/>. This seems to solve this issue, the last of those critical ML problems we found. The issues list currently contain only one small annoying bug regarding parameter naming (#74<#74>), and the rest are new/improved features. Therefore, I think we should release version 1.0.1 with the current fixes before starting work on the next thing, what do you say? — Reply to this email directly, view it on GitHub<#85 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AFPRUPVMSEZHV7Q5L5N4KTTV2SCEZANCNFSM5PKCIQMA>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[wleoncio]

Sure thing, we can go over this on Friday.

In the meantime, you can install the latest development version with remotes::install_github("ocbe-uio/TruncExpFam") if you wish to take it for a spin.

W

[rho62]

Thanks! /René Fra: Waldir Leoncio ***@***.***> Svar til: ocbe-uio/TruncExpFam ***@***.***> Dato: tirsdag 23. august 2022 kl. 10:28 Til: ocbe-uio/TruncExpFam ***@***.***> Kopi: Rene Holst ***@***.***>, Mention ***@***.***> Emne: Re: [ocbe-uio/TruncExpFam] ML estimation not working for some cases of the Beta distribution (Issue #85) Sure thing, we can go over this on Friday. In the meantime, you can install the latest development version with remotes::install_github("ocbe-uio/TruncExpFam") if you wish to take it for a spin. W — Reply to this email directly, view it on GitHub<#85 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AFPRUPR3KHJAILHCPYKTIXDV2SDRXANCNFSM5PKCIQMA>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[rho62]

Hi Waldir Yes. We use the first parametrization; the one with (alpa,beta)=(eta1, eta2) . This is the same one as used in R. I checked the calculations and found no errors in them I’ll try to investigate more on the error. But now I have to run for some other meetings – yak Btw, I looked at an older version where the two transformation functions looked like this: #' @export natural2parameters.trunc_beta <- function(eta) { # eta: The natural parameters in a beta distribution # returns (alpha,beta) return(c(shape1 = eta[1], shape2 = eta[2])) } #' @export parameters2natural.trunc_beta <- function(parms) { # parms: The parameters shape and rate in a beta distribution # returns the natural parameters return(c(shape1 = parms[1], shape2 = parms[2])) } Notice that these are just bivariate functions and do not rely on any samples. That’s how it should be. I’ll be in DM tomorrow. Perhaps we can talk then. /René Fra: Waldir Leoncio ***@***.***> Svar til: ocbe-uio/TruncExpFam ***@***.***> Dato: torsdag 7. april 2022 kl. 09:41 Til: ocbe-uio/TruncExpFam ***@***.***> Kopi: Subscribed ***@***.***> Emne: Re: [ocbe-uio/TruncExpFam] ML estimation not working for some cases of the Beta distribution (Issue #85) There's a nice table on https://en.wikipedia.org/wiki/Exponential_family#Table_of_distributions that suggests two variants of the Beta: one has no transformation and the other has a ±1 transformation between natural and inverted mapping. Neither solution solves this problem, so it might be somewhere else. — Reply to this email directly, view it on GitHub<#85 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AFPRUPQEAMMDFIOZLCODSM3VD2GR3ANCNFSM5PKCIQMA>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>