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