0

I am generating two distributions (say, mean1, sd1 and mean2, sd2), using convolution on them and then calculating the integral for the newly generated distribution from a given threshold (a given maximum value) to Inf. Here, I am using functions from library distr.

Interestingly, when integrating, if my original distributions present values fairly larger than the threshold the output of the integral = 0 (should be = 100, indicating that the whole distribution falls after the threshold)

MWE

library(distr)

# Two distributions
mean1 <- 560 
sd1   <- 4.1

mean2 <- 1570  # values to be tested: # 1570 # 1571 # 1970
sd2   <- 3.219

# Getting Normal Distribution for both groups
Np1  <- Norm(mean = mean1, sd = sd1)  
Np2  <- Norm(mean = mean2, sd = sd2) 

# Convolution of distributions
conv <- convpow(Np1 + Np2, 1) 

# Sanity check
plot(conv)

# Distribution function
f.Z  <- d(conv) 

# This works well with any values. 
min_threshold <- 1925.3
integrate(f.Z, 0, min_threshold, abs.tol = 0L)$value * 100   

# The statement below provides abnormal output depending on the tested value presented above for mean2:
max_threshold <- 1968.9
integrate(f.Z, max_threshold, Inf, abs.tol = 0L)$value * 100

Here, I would expect that any value for mean2 > 1570 would have 100 as output, indicating a 100% chance of the distribution being after the max_threshold. Instead, output sharply decreases after 1570 until it reaches 0. Would anyone understand this behavior? Am I missing something here?

My machine has Windows 10, 64, R version 4.0.3, and library distr 2.8.0.

This code is based on the code provided here: integrate() in R gives terribly wrong answer

Many thanks!

MRC
  • 1
  • 1
  • The second `integrate` call is giving me `[1] 100`. R4.0.3 on Ubuntu 20.04, distr 2.8.0. – Rui Barradas Oct 30 '20 at 16:54
  • Would you please confirm that when you run with mean2 = 1970, the second integrate: max_threshold <- 1968.9 integrate(f.Z, max_threshold, Inf, abs.tol = 0L)$value is giving output = 1? For me, in the referred machine I have output = 0. Thanks for your feedback!! – MRC Oct 30 '20 at 18:33
  • You are right, it fails. But this one returns 1: `cubature::hcubature(f.Z, max_threshold, Inf, tol = .Machine$double.eps^0.5)$integral`. See [this SO post](https://stackoverflow.com/questions/62122505/integrate-gives-totally-wrong-number?noredirect=1&lq=1) for more examples on how to solve similar problems. – Rui Barradas Oct 30 '20 at 18:44

0 Answers0