Commit e1cb0154 authored by kunstler's avatar kunstler
Browse files

ratio CI

parent 8593643c
......@@ -29,20 +29,30 @@ mvrnorm <-
## define param
mean.v <- c(0.15, 0.25)
Sigma <- matrix(c(0.02,0.01,0.01,0.05),2,2)
Sigma <- matrix(c(0.002,0.001,0.001,0.005),2,2)
# MC
res.mc <- mvrnorm(n = 10000, mean.v, Sigma)
ratio.mc <- res.mc[,1]/res.mc[,2]
hist(ratio.mc, breaks = 100)
lines(v = quantile(ratio.mc, probs = c(0.025, 0.975)))
par(mfrow = c(2,2))
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))
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
##because lot of data approximate teh student by a normal
ta <- qnorm(0.025)
g <- mean.v[1]*mean.v[2] - ta*Sigma[1,2]
d <- mean.v[2]^2 - ta*Sigma[2,2]
b <- mean.v[1]^2 - ta*Sigma[1,1]
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)
q.l <- (g-ci.b)/d
q.h <- (g+ci.b)/d
abline(v = c(q.l, q.h), col = 'green')
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment