4
$\begingroup$

This code was kindly recommended to me in my original question. It returned the same parameter estimates as the software called CRAFT by Aon Benfield. I have also managed to replicate it for the Weibull distribution.
I was wondering if anyone could help me replicate if for Pareto? Based on CRAFT, I would expect a shape of ~0.59 and scale of ~55.15 (b~93 and q~1.7). Finally, I would also like to do the same but using the SSQ method instead of MLE?

My data = c(28.744, 385.714, 20.595, 99.350,31.864, 77.713, 264.408, 21.204, 31.937, 0.900, 18.762, 173.276, 23.707) 


constrained_mle <- function(x, p, q) {
  lnL <- function(shape) {
    # Solving q = qgamma(p, shape)*scale for the necssary scale
    scale <- q/qgamma(p, shape)
    sum(dgamma(x, shape, scale=scale, log=TRUE))
  }
  res <- optimise(lnL, lower=0, upper=1e+3, maximum=TRUE, tol=1e-8)
  scale <- q/qgamma(p, res$maximum)
  c(scale=scale, shape=res$maximum)
}
par <- constrained_mle(x, .95, 500.912)
par
#>       scale       shape 
#> 231.0561574   0.6038756
# checking that the solution is correct
qgamma(.95, scale=par[1], shape=par[2])
#> [1] 500.912 ```
Improve this question
$\endgroup$
8
  • $\begingroup$ can you be more specific? $\endgroup$
    – utobi
    Commented May 5, 2023 at 10:54
  • $\begingroup$ So I have a data set: loss <- c(28.744,385.714,20.595,99.350,31.864,77.713,264.408,21.204,31.937,0.900,18.762,173.276,23.707). I know I can fit a gamma curve to this by using fitdist(loss,distr="gamma",method="mle"). However, I would also like the fitted curve to be forced through the value of 500.612 at a percentile of 95.646%. I would then like to obtain the parameters/view the plot of this adjusted gamma fit. $\endgroup$
    – Tom
    Commented May 5, 2023 at 11:21
  • 1
    $\begingroup$ Hi Tom, please add any relevant information by editing your post. Comments in the future may get removed or may be overlooked. $\endgroup$
    – utobi
    Commented May 5, 2023 at 11:54
  • 1
    $\begingroup$ Welcome to CV, Tom. That might not be a great way to shape the tail. But if you're sure that's what you want to do, then you likely are adept enough with whatever your fitting method might be to frame it as an optimization problem and solve it subject to this constraint. Evidently there are plenty of details lurking here that you need to spell out before we can provide any specific help: in particular, what method you are using and what distribution family specifically you want to fit. Please, then, edit your post to supply adequate information about your problem. $\endgroup$
    – whuber ♦
    Commented May 5, 2023 at 11:58
  • $\begingroup$ I have added some more relevant information to the post $\endgroup$
    – Tom
    Commented May 5, 2023 at 12:46

2 Answers 2

Reset to default
5
$\begingroup$

You can fit a gamma distribution in R with optim, e.g. by minimizing the sum of squared errors for the mean and the 95th percentile, which actually sets both errors to zero:

obs = c(28.744,385.714,20.595,99.350,31.864,77.713,       
264.408,21.204,31.937,0.900,18.762,173.276,23.707)
tail = 500.912

error = function(x)(qgamma(0.95, shape=x[1], scale=x[2]) - tail)^2 + (x[1]*x[2]-mean(obs))^2

params = optim(c(1,1), error)$par
params

The code uses the formula for the gamma’s mean and gets the desired mean and percentile exactly with a shape parameter of 0.147 and a scale parameter of 616.6.

You can then check the quantile with

qgamma(0.95, shape=params[1], scale=params[2])

You can also check the log-likelihood of the result with

sum(log(dgamma(shape = params[1], scale = params[2], obs)))

and compare it with the log-likelihood from other two-parameter distributions.

Improve this answer
$\endgroup$
10
  • $\begingroup$ Thanks for this Matt. If you were then to check if this has worked by labelling the optim and using qgamma, it does not return 500 for the 95th percentile? $\endgroup$
    – Tom
    Commented May 5, 2023 at 14:29
  • 1
    $\begingroup$ I get 500.8 from qgamma(0.95, shape=0.147, scale=616.6) $\endgroup$
    – user225256
    Commented May 5, 2023 at 15:48
  • $\begingroup$ @Tom Increasing the penalty/cost for the difference in the tail can get the 95th percentile closer. Alternatively you could place a fixed boundary although that might make the optimization algorithm run badly $\endgroup$
    – Sextus Empiricus
    Commented May 5, 2023 at 15:55
  • $\begingroup$ Thanks Matt, I can now see that. I want to do the same with a weibull distribution. How could I then compare the two with regards to their fit to the original data? $\endgroup$
    – Tom
    Commented May 5, 2023 at 16:17
  • 1
    $\begingroup$ I see -- so why not describe your solution exactly that way, since the use of the loss function is merely a computational tool? Indeed, since the mean of any Gamma function is a simple function of its shape and scale parameters, why even do the computation with two parameters? Just use one parameter and adjust it to match the percentile. That's a simple and much more efficient root-finding problem. $\endgroup$
    – whuber ♦
    Commented May 5, 2023 at 19:55
5
$\begingroup$

For distribution with a scale parameter such as the gamma and the weibull, the likelihood can be maximised numerically with respect to the remaining parameters (the shape in case of the gamma and the weibull) since for any given value of those parameters (the shape), there will be only a single value of the scale parameter for which the 95% quantile of the distribution matches the target value. The following R code implements this method:

x <- c(28.744,385.714,20.595,99.350,31.864,77.713,       
  264.408,21.204,31.937,0.900,18.762,173.276,23.707)
constrained_mle <- function(x, p, q) {
  lnL <- function(shape) {
    # Solving q = qgamma(p, shape)*scale for the necssary scale
    scale <- q/qgamma(p, shape)
    sum(dgamma(x, shape, scale=scale, log=TRUE))
  }
  res <- optimise(lnL, lower=0, upper=1e+3, maximum=TRUE, tol=1e-8)
  scale <- q/qgamma(p, res$maximum)
  c(scale=scale, shape=res$maximum)
}
par <- constrained_mle(x, .95, 500.912)
par
#>       scale       shape 
#> 231.0561574   0.6038756
# checking that the solution is correct
qgamma(.95, scale=par[1], shape=par[2])
#> [1] 500.912

The constained MLEs are evidently quite different from the estimates produced by matching the first moment and the target quantile (shape=0.1470374 and scale=616.5818235) in the answer by Matt F.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Hi Jarle, this is brilliant thank you. So I assume the same function could apply to Weibull, and what about Pareto2? $\endgroup$
    – Tom
    Commented May 9, 2023 at 15:11
  • $\begingroup$ @Tom Yes, sure. But note that for the Pareto, the constained likelihood may be zero for some values of the shape parameter since the support is of the pdf is $[x_m,\infty)$, that is, it depends on the scale parameter $x_m$. $\endgroup$
    – Jarle Tufto
    Commented May 9, 2023 at 15:51
  • $\begingroup$ When trying the function for Weibull, I receive this error: Warning messages: 1: In optimise(lnL, lower = 0, upper = 1000, maximum = TRUE, ... : NA/Inf replaced by maximum positive value. When trying it with pareto, i receive this error: Error in qpareto(p, shape) : argument "shape" is missing, with no default. Sorry for the basic questions, I am a beginner in r! $\endgroup$
    – Tom
    Commented May 9, 2023 at 16:06
  • 1
    $\begingroup$ I have managed to replicate it for Weibull but I am still struggling with pareto? $\endgroup$
    – Tom
    Commented May 12, 2023 at 9:19
  • $\begingroup$ Also, is there an alternative function that I could use for the SSQ method? $\endgroup$
    – Tom
    Commented May 12, 2023 at 9:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.