r - Estimating the Standard Deviation of a ratio using Taylor expansion -


i interested build r function can use test limits of taylor series approximation. aware there limits doing, it's limits wish investigate.

i have 2 distributed random variables x , y. x has mean of 7 , standard deviation (sd) of 1. y has mean of 5 , sd of 4.

me.x <- 4; sd.x <- 1 me.y <- 5; sd.y <- 4

i know how estimate mean ratio of y/x, 

# e(y/x) = e(y)/e(x) - cov(y,x)/e(x)^2 + var(x)*e(y)/e(x)^3 me.y/me.x - 0/me.x^2 + sd.x*me.y/me.x^3 [1] 1.328125

i stuck on how estimate standard deviation of ratio? realize have use taylor expansion, not how use it.

doing simple simulation 

 x <- rnorm(10^4, mean = 4, sd = 1);  y <- rnorm(10^4, mean = 5, sd = 4)  sd(y/x)  [1] 2.027593  mean(y/x)[1]  1.362142

there analytical expression pdf of ratio of 2 gaussians, done david hinkley (e.g. see wikipedia). compute momentums, means etc. typed , apparently doesn't have finite second momentum, doesn't have finite standard deviation. note, i've denoted y gaussian x, , x y (formulas assume x/y). i've got mean value of ratio pretty close you've got simulation, last integral infinite, sorry. sample more , more values, sampling std.dev growing well, noted @g.grothendieck

library(ggplot2)  m.x <- 5; s.x <- 4 m.y <- 4; s.y <- 1  <- function(x) {     sqrt( (x/s.x)^2 + (1.0/s.y)^2 ) }  b <- function(x) {     (m.x*x)/s.x^2 + m.y/s.y^2 }  c <- (m.x/s.x)^2 + (m.y/s.y)^2  d <- function(x) {     u <- b(x)^2 - c*a(x)^2     l <- 2.0*a(x)^2     exp( u / l ) }  # pdf ratio of 2 different gaussians pdf <- function(x) {     r <- b(x)/a(x)     q <- pnorm(r) - pnorm(-r)      (r*d(x)/a(x)^2) * (1.0/(sqrt(2.0*pi)*s.x*s.y)) * q + exp(-0.5*c)/(pi*s.x*s.y*a(x)^2) }  # normalization nn <- integrate(pdf, -inf, inf) nn <- nn[["value"]]  # plot pdf p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x)) p <- p + stat_function(fun = function(x) pdf(x)/nn) + xlim(-2.0, 6.0) print(p)  # first momentum m1 <- integrate(function(x) x*pdf(x), -inf, inf) m1 <- m1[["value"]]  # mean print(m1/nn)  # sampling set.seed(32345) n <- 10^7l x <- rnorm(n, mean = m.x, sd = s.x); y <- rnorm(n, mean = m.y, sd = s.y) print(mean(x/y)) print(sd(x/y))  # second momentum - infinite! m2 <- integrate(function(x) x*x*pdf(x), -inf, inf)

thus, impossible test taylor expansion std.dev.


Comments

Post a Comment

Popular posts from this blog

java - pagination of xlsx file to XSSFworkbook using apache POI -

right in code, reading xlsx file, xssfworkbook, , writing database. but, when size of xlsx file increases, causes outofmemory error. can not increase server size, or divide xlsx file pieces. tried loading workbook using file (instead of inputstream), didn't either. i looking way read 10k rows @ time (instead of entire file @ once) , iteratively write workbook , database. is there way apache poi? poi contains called "eventmodel" designed purpose. it's mentioned in faq : the ss eventmodel package api reading excel files without loading whole spreadsheet memory. require more knowledge on part of user, reduces memory consumption more tenfold. based on awt event model in combination sax. if need read-only access, best way it. however, may want double check first if issue somewhere else. check out this item : i think poi using memory! can do? 1 comes quite lot, reason isn't might think. so, first thing check - what's source of prob...
Read more

Unlimited choices in BASH case statement -

consider simple bash snippet: case $option in 1) image=${options[0]%.tar} ;; 2) image=${options[1]%.tar} ;; 3) image=${options[2]%.tar} ;; 4) image=${options[3]%.tar} ;; *) echo "invalid option" exit 1 esac in real script numbers goes 30. makes pretty long. can somehow specify cases variable? something this: case $option in $i) image=${options[$(($i-1))]%.tar} any pointer appreciated. you can match against multiple patterns in single clause: case $option in 1|2|3: echo "$option one, 2 or three" ;; esac if still typing you, can use simple pattern matching: case $option in # following matches # - single-digit numbers 0-9 # - two-digit numbers starting either 1 or 2 # - number 30 [0-9]|[12][0-9]|30) image=${options[$(($option-1))]%.tar} ;; *) echo "invalid option" 1>&2 exit 1 esac
Read more

apache - How do I stop my index.php being run twice for every user -

we've got webapp initial html served using apache , php. i noticed index.php being run twice every request. seems caused / effected rewrite rules in our .htaccess - rewritecond %{http_host} ^(www\.|hotels\.)ourdomain\.com$ [nc,or] rewritecond %{server_name} ^(www\.|hotels\.)ourdomain\.com$ [nc] rewriterule (.*) http://ourdomain.com/$1 [r=301] rewriterule ^hotels/([^/]+)/?\??(.*)$ ?d=$1&$2 [qsa] the final rule moving parameter url path query string. i.e. http://ourdomain.com/hotels/vegas?someparam=1 becomes http://ourdomain.com?d=vegas&someparam=1 if go directly query string version index.php run once. if use url redirected, index.php run twice (i'm checking adding error_log('end of index.php') file). so example going http://ourdomain.com/hotels/paris hits file twice http://ourdomain.com?d=paris hits once. i've seen this question , had @ blog mentioned , , can't find empty string url's (i've tried using yslow process). ...
Read more