I am trying to find this integral in R. But I am getting 0 because the constant part is close to 0 (but not zero). How can I solve this problem without getting the final answer of 0? Please see the code I tried. None of them is working.
set.seed(1)
options(digits = 16)
M=1000 #Replicate
n=100 #sample size
s=1.5 # sample standard deviation
xbar=18 # sample mean
alpha=0.05 # significance level
observed_values<- function(xbar,s,n,M,alpha) {
r=list()
for (i in 1:M){
Z=rnorm(1,0,1)
U=sqrt(rchisq(1,df=n-1))
r[[i]]=xbar+(Z*sqrt(n-1)/U)*(s/sqrt(n))+.5*((n-1)*s^2)/(U^2)
}
r2=sort(unlist(r))
quantile(r2,na.rm = T,probs = c(alpha/2,1-alpha/2))
}
values1=observed_values(xbar,s,n,M,alpha)
values2=log(observed_values(xbar,s,n,M,alpha))
> values2
# 2.5% 97.5%
# 2.930742746832908 2.976810026452883
#### Integration
fun0 = function(w,t,xbar,s,n,ll,ul){
(sqrt(n)/(sqrt(2*pi)*gamma((n-1)/2)*2^((n-1)/2)*s*sqrt(n-1)))*(w^((n/2)-1))*
exp(-((((t-xbar-(((n-1)*s^2)/(2*w)))^2)*(n/((n-1)*s^2))+1)*w/2))
}
fun01 <- function(xbar,s,n,ll,ul) {
integrate(function(w) {
sapply(w, function(w) {
integrate(function(t) fun0(w,t,xbar,s,n,ll,ul),0,Inf)$value
})
}, ll, ul)$value
}
fun01(xbar,s,n,values2[1],values2[2])
# [1] 8.352234203850945e-190
fun0 = function(w){
(sqrt(n)/(sqrt(2*pi)*gamma((n-1)/2)*2^((n-1)/2)*s*sqrt(n-1)))*(w[2]^((n/2)-1))*
exp(-((((w[1]-xbar-(((n-1)*s^2)/(2*w[2])))^2)*(n/((n-1)*s^2))+1)*w[2]/2))
}
pcubature(fun0, c( 0, ll), c(Inf, ul))
$integral
[1] 0
$error
[1] 0
$functionEvaluations
[1] 9
$returnCode
[1] 0
> hcubature(fun0, c( 0, ll), c(Inf, ul))
$integral
[1] 0
$error
[1] 0
$functionEvaluations
[1] 85
$returnCode
[1] 0
Do I need to use double precision? Or can anyone suggest any other methods?
