0

I would like to create a function in order to integrate the product of 2 distributions (a PDF and CDF) with respect to x. This is in order to perform a probabilistic risk assessment and determine an 'expected total risk'.

My code I have so far, which does not work, is...

F <- function(x) {dgamma(x, shape=0.4259325, rate=8.1490741)*plnorm(x, meanlog=3.4906738, sdlog=0.2938556)}
c <- integrate( F, lower = 0, upper = inf)

Where and what have I done wrong here?

Additionally, Ive fitted my distributions using 'fitdist' function and for now am hard coding the distribution parameters to try and get the integration above working. However, how would I take the shape/rate/meanlog etc directly from the fitted distribution - the 'Fitted' object in the example code below?

x_vect <- seq(0.1,10, length.out = 1e3)
example <- rgamma(x_vect, shape=0.4259325, rate=8.1490741)
plot(x_vect, example, type = "l", log="x")
Fitted <- fitdist(example, "gamma") 
plot(Fitted)
Gareth Le Page
  • 37
  • 7
  • Try changing the upper limit of the integration to Inf (lower case inf is not infinity in R) – Calvin Aug 03 '18 at 08:48
  • 1
    If you follow that advice, you'll see that precision of the integral is not very good. It might be sufficient to calculate a finite interval. E.g., check out `curve(F(x), from = 0, to = 10)`. It might be even better to talk to a mathematician or statistician. – Roland Aug 03 '18 at 08:51
  • You should also scale your function. Values of 1e-27 magnitude result in numerical issues. – Roland Aug 03 '18 at 08:55
  • Damn - should have spotted the I in Inf. thanks. @Roland by scale do you mean multiply the whole function within the code? and if so, do you know at what point do the numerical issues disappear? – Gareth Le Page Aug 03 '18 at 09:11
  • Do you really need product of two PDFs? Or you have product of two random values? – Severin Pappadeux Aug 03 '18 at 15:07

0 Answers0