progress on CI of ratio with Fieller method

R/analysis/ratio.CI.R

 ##
 
-qt(p= 0.025, df = 100000)
-qnorm(0.025)
 
 library(MASS)
 
@@ -29,7 +27,7 @@ mvrnorm <-
 
 ## define param
 mean.v <- c(0.15, 0.25)
-Sigma <- matrix(c(0.002,0.001,0.001,0.005),2,2)
+Sigma <- matrix(c(0.002,0.0001,0.0001,0.005),2,2)
 # MC
 res.mc <- mvrnorm(n = 10000, mean.v, Sigma)
 ratio.mc <- res.mc[,1]/res.mc[,2]
@@ -38,21 +36,23 @@ hist(res.mc[,1], breaks = 100)
 abline(v = quantile(res.mc[, 1], probs = c(0.025, 0.975)))
 hist(res.mc[,2], breaks = 100)
 abline(v = quantile(res.mc[, 2], probs = c(0.025, 0.975)))
-hist(ratio.mc, breaks = 1000, xlim = c(0, 1.4))
+hist(ratio.mc, breaks = 1000, xlim = c(0, 1.5))
 abline(v = quantile(ratio.mc, probs = c(0.025, 0.975)))
 abline(v = mean(res.mc[, 1])/mean(res.mc[, 2]), col = 'red')
 
-## Fieller's theorem
+# compute Ci of ratio based on Fieller's methods
+ratio.ci <- function(mean.v, Sigma, alpha = 0.05){
+ta <-  -qnorm(alpha/2)
+theta <- mean.v[1]/mean.v[2]
+k <-  ta^2*Sigma[2,2]/(theta)^2
+mm <- theta+(k/(1-k))*(theta+Sigma[1,2]/Sigma[2,2])
+ci <- ta/(mean.v[2]*(1-k))*sqrt(Sigma[1,1]+2*theta^2*Sigma[2,2] -
+                                k*(Sigma[1,1]- Sigma[1,2]^2/Sigma[2,2]))
 
+return(c(q.l = mm-ci, q.h =  mm+ci))
+}
 
-##because lot of data approximate teh student by a normal
-ta <-  qnorm(0.025)
-g <- mean.v[1]*mean.v[2] - ta^2*Sigma[1,2]
-d <- mean.v[2]^2 - (ta*Sigma[2,2])^2
-b <- mean.v[1]^2 - (ta*Sigma[1,1])^2
-
-ci.b <- sqrt(g^2-d*b)
+abline(v = ratio.ci(mean.v, Sigma), col = 'green')
 
-q.l <- (g-ci.b)/d
-q.h <- (g+ci.b)/d
-abline(v = c(q.l, q.h), col = 'green')
+quantile(ratio.mc, probs = c(0.025, 0.975))
+ratio.ci(mean.v, Sigma)