n <- 30
s <- 1.8
r <- 1/2.2
G <- 0.5
E <- s*(1/r)
D <- G^2 + s*(1/r^2)
x <- seq(-20, 20, 0.1)
t <- c(rnorm(n, 0, 0.5))+c(rgamma(n, shape=s, rate=r))
dfun3 <- function(x) {
  dfun <- function(g,y) {
  dnorm(g, 0, 0.5) * dgamma(y-g, shape=s, rate=r)}
  b <- x;
  for (i in 1:length(x))
    b[i] <- integrate(dfun, lower = -Inf, upper = Inf, y = x[i])$value
  return(b)
}
pfun3 <- function(x){
for (i in 1:length(x)) {
  b[i] <- integrate(dfun3, -Inf, x[i])$value
  return(b)
  }
}
plot(x, dfun3(x), type='l', main='Normal+Gamma pdf', xlab='', ylab='')
lines(rep(E,2), c(0, dfun3(E)), col='red', lty=2)
lines(rep((E-sqrt(D)),2), c(0, dfun3((E-sqrt(D)))), col='blue', lty=2)
lines(rep((E+sqrt(D)),2), c(0, dfun3((E+sqrt(D)))), col='blue', lty=2)
legend('topleft', legend=c(sprintf( 'E = %s', round(E, 2) ),sprintf( 'D = %s', round(D, 2) )), inset=0.03, box.lty=0, col=c('red', 'blue'), lty=2)
plot(x, pfun3, type='l', main='Normal+Gamma cdf', xlab='', ylab='')
lines(rep(E,2), c(0, pfun3(E)), col='red', lty=2)
lines(rep((E-sqrt(D)),2), c(0, pfun3(E-sqrt(D))), col='blue', lty=2)
lines(rep((E+sqrt(D)),2), c(0, pfun3(E+sqrt(D))), col='blue', lty=2)
legend('topleft', legend=c(sprintf( 'E = %s', round(E, 2) ),sprintf( 'D = %s', round(D, 2) )), inset=0.03, box.lty=0, col=c('red', 'blue'), lty=2)
hist(t, freq=FALSE, main='30 random values', xlab='', ylab='')
lines(x, dfun3(x), col='red')
plot(ecdf(t), main = 'Empirical cdf', xlab='', ylab='')
lines(x, pfun3, col = 'red')