Distribution of ratio of sample means from two independent normal variables?

Question

The Question

We have a sample of size $N$ with mean $\bar{x}$ and SD $\bar{\sigma_x}$ from a random variable $X \sim \mathcal{N} (\mu, \sigma^2)$

We have a sample of size $M$ with mean $\bar{y}$ and SD $\bar{\sigma_y}$ from a random variable $Y \sim \mathcal{N} (c\mu, c^2\sigma^2)$

We wish to find estimates of $\mu$ and $c$, along with the distributions of those estimates

My Progress So Far

• $\bar{x}$ is one obvious estimate of $\mu$, and we know it has a t-distribution, but using only this statistic ignores the information about $\mu$ contained in $\bar{y}$
• $\frac{\bar{y}}{\bar{x}}$ would give us an estimate of $c$, but what distribution would it have? $\frac{Y}{X}$ has a Cauchy distribution, but what is the analogous distribution when using the ratio of sample means? Said another way:

Normal : t-distribution :: Cauchy : ???

• Once we have an estimate of $c$, we could divide $\bar{y}$ by that estimate to get another estimate of $\mu$, and thereby extract the additional info about $\mu$ contained in $\bar{y}$. But what is the distribution of that estimate, and how to combine it with our $\bar{x}$ estimate? Things seem to be getting confusing since now we'd have a t-distribution random variable divided by our sample-equivalent-of-a-Cauchy random variable.... Is there a more straightforward way to do this?

Answer 1

This framework is a particular case of Cox's model

http://www.jstor.org/stable/2530661

studied here

http://onlinelibrary.wiley.com/doi/10.1002/bimj.200310009/abstract

Answer 2

If you could only divide Y by c, all of your data would come from $N(\mu, \sigma^2)$. This suggests to me an iterative approach. Estimate c, then use the pooled data to estimate $\mu$ and $\sigma^2$; then use these improved estimates to get a better estimate of c, and repeat until it converges. This ducks the question of the theoretical best estimator but might still be a useful approach.

You could use simulation based on your model (if you are confident in it) to work out the approximate distribution of any estimator, or mix of estimators, you choose.

Then I'd use a bootstrap to estimate the variances of your estimates of $\mu$ and $c$. This has the advantage of not being as dependent on the distributional assumptions of your model.

It's easier for me to illustrate this general approach than to try to explain:

###
# Create a function that does the iterative thing
RatioEst <- function(x,y, verbose=FALSE){
    mu_latest <- mean(x)
    sigma2_latest <- var(x)
    for (i in 1:5){
        c_latest <- mean(c(
            mean(y / mu_latest),
            sqrt(var(y)/sigma2_latest)))
        mu_latest <- mean(c(x, y/c_latest))
        sigma2_latest <- var(c(x, y/c_latest))
        if(verbose){print(c(mu_latest, c_latest, sigma2_latest))}
    }
    return(c(mu_latest, c_latest))
}

#### Simulation to get an idea of the distribution of estimates.
# Simulate data many times and see the results of our estimation technique. 
# True values of mu and c are 30 and 2

reps <- 10000
results <- matrix(0, nrow=reps, ncol=2)

for (i in 1:reps){
    x <- rnorm(20,30,5)
    y <- rnorm(30,60,10)
    results[i,] <- RatioEst(x,y, verbose=FALSE)
}

summary(results)

par(mfrow=c(1,2))
plot(density(results[,1]), bty="l", main="Simulated estimates of mu",
    xlab="True value=30")
plot(density(results[,2]), bty="l", main="Simulated estimates of c",
    xlab="True value=2")

This gives the results below which suggest that the estimators I've chosen are biased (for mu upwards; for c downwards) although the median of repeated estimates is very good.

       mu              c        
 Min.   :24.43   Min.   :0.5937  
 1st Qu.:28.85   1st Qu.:1.8256  
 Median :30.01   Median :2.0072  
 Mean   :31.21   Mean   :1.9340  
 3rd Qu.:31.87   3rd Qu.:2.1284  
 Max.   :73.57   Max.   :2.6688 

So that was a simulation to show the properties of the estimators I'd chosen (which you'll see included a funny sort of estimate of c that is an average of two estimates). Now below is how you'd go about the actual estimation, if you used this approach:

#### Actual estimation
set.seed(123)
x <- rnorm(20,30,5)
y <- rnorm(30,60,10)

# point estimates
RatioEst(x, y, verbose=TRUE)

which gives these results (including showing how the iteration works):

[1] 31.12087  1.89926 22.66501
[1] 31.050508  1.906381 22.529121
[1] 31.001155  1.911407 22.438041
[1] 30.967360  1.914864 22.377693
[1] 30.944615  1.917198 22.337999
[1] 30.944615  1.917198

To get a confidence interval here is the bootstrap:

# bootstrap
# Simulate data *once* and then resample from it many times.
# Has the advantage that will work even if original specification
# of distribution is incorrect
reps <- 699
boot.results <- matrix(0, nrow=reps, ncol=2)
for (i in 1:reps){
    boot.results[i,] <- RatioEst(
        x=sample(x, replace=TRUE), 
        y=sample(y, replace=TRUE))
}
summary(boot.results)
apply(boot.results, 2, quantile, probs=c(0.025, 0.975))

which gives these results for a (non symmetrical) 95% confidence interval:

           mu      c
2.5%  28.02008 1.109987
97.5% 44.38868 2.236229