6
$\begingroup$

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?
$\endgroup$
4
  • $\begingroup$ An additional opportunity/complication you don't mention is that you could also use the sample variances to estimate c (and doubtless would get a different estimate to that which comes from the sample means). $\endgroup$ Jan 5, 2013 at 5:00
  • $\begingroup$ Yes that is a good point. The real question is what is the best way to extract every drop of information possible about $\mu$ and $c$ from our samples, while knowing their distributions so that we can create confidence intervals for our estimates $\endgroup$
    – Jonah
    Jan 5, 2013 at 5:17
  • $\begingroup$ Are the two samples independent of each other, too? $\endgroup$
    – cardinal
    Jan 5, 2013 at 19:19
  • $\begingroup$ Yes, they are independent $\endgroup$
    – Jonah
    Jan 5, 2013 at 19:59

2 Answers 2

4
$\begingroup$

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

$\endgroup$
3
  • 1
    $\begingroup$ Can you briefly summarize in your answer the content contained in the papers in case the link goes dead? $\endgroup$
    – jonsca
    Jan 6, 2013 at 0:37
  • $\begingroup$ Hmmm... That is not so easy to do, actually. If you'd like to take a shot at it, I will happily send you the article though. $\endgroup$
    – Jonah
    Jan 6, 2013 at 17:46
  • $\begingroup$ For posterity, the reference to the article in the first link is "Biometrics, Vol. 41, No. 1 (Mar., 1985), pp. 261-265, Interval Estimates for the Ratio of the Means of Two Normal Populations with Variances Related to the Means" $\endgroup$
    – Jonah
    Jan 6, 2013 at 17:53
3
$\begingroup$

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 

enter image description here

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
$\endgroup$
2
  • $\begingroup$ Hi Peter, can you elaborate on your idea of bootstrapping to estimate the variances of $\mu$ and $c$? Also that still leaves open the question of what distribution our estimates would be following... what were you assuming for that? $\endgroup$
    – Jonah
    Jan 5, 2013 at 5:13
  • $\begingroup$ I have greatly expanded this with a demo of the bootstrapping and also of simulation to test the distribution of the estimates, if your model is correct. Hopefully this complements the approach suggested by @Leotardo. $\endgroup$ Jan 7, 2013 at 6:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.