Yes, if your alternative is a shift alternative, the corresponding estimate is the Hodges-Lehmann.
The estimate is easy enough.
The Hodges Lehmann estimate is the sample median of the cross-sample pairwise differences.
For really large samples, there are efficient methods (for example, ones that sort the two samples and then use the information in that to speed up the search for the median).
I'm going to assume your sample is not so large that its worth the effort of doing that, in which case, you simply want to find all pairwise differences and find the median of the result.
For example, if your x-sample is:
16.388 6.775 18.270 17.034 18.825 16.197 16.709 18.188 9.999
and your y-sample is:
9.305 9.799 11.103
then the matrix of differences, x-y is
7.083 6.589 5.285
-2.530 -3.024 -4.328
8.965 8.471 7.167
7.729 7.235 5.931
9.520 9.026 7.722
6.892 6.398 5.094
7.404 6.910 5.606
8.883 8.389 7.085
0.694 0.200 -1.104
(here columns correspond to particular y's and rows to particular y's). You then
take the median of the whole set of differences.
(In R this is just a matter of:
> median(outer(x,y,"-"))
[1] 6.91
It's doubtless just about as easy in Matlab; would something like [X,Y] = meshgrid(x,y); Z = X - Y; work to generate the matrix of differences?)
If you have samples so large that this isn't feasible, you should indicate this.
I don't know why it didn't click before, but the basic concept
of finding an interval for something like the Hodges-Lehmann shift
parameter in the Wilcoxon-Mann-Whitney (WMW) test is pretty obvious;
the same idea is widely used for other permutation tests.
The underlying idea is this: you shift the second sample left and right
until you hit the $\alpha/2$ and $1-\alpha/2$ quantiles of the WMW statistic
(being discrete, these will cut off less than $\alpha/2$ in the tails, but
you can legitimately have one a bit larger than $\alpha/2$ if the sum of
tail areas is less than $\alpha$; however, most people will prefer symmetry
even at the cost of a more-conservative-than-necessary interval).
This approach can be easily implemented in software via use of a canned
root-finding algorithm that doesn't require derivatives. Bounds are obvious
enough (if you haven't hit it by the time you get complete separation
of samples, you never will, so check those 'maximum shift' values first
and then use them as bounds; similarly, the HL estimate itself is the bound
for the other end of each search).
There are other approaches to doing this that will probably be more efficient,
but this should work well enough with a good root-finder (perhaps involving
only a dozen or so evaluations of statistics at different shift values).
In response to the comments about the issue being a change of scale, not location:
One possibility - you could consider looking at estimating the location shift in the logs, which should correspond to the log of the scale change (ratio of scales) in the untransformed values. That is if $X∼F(x)$ and $Y∼F(y/a)$ then if $X^∗=\ln X∼G(x), Y^∗=\ln Y∼G(y−\ln a)$, where $G(x)= F(\exp(x))$.
Intervals, at least should carry over nicely.
Here's an example of that in use:
The data are as follows. The x sample is 60 observations from a distribution with scale parameter $\mu=18$, while the y sample is 50 observations with scale parameter $\mu=25$. They have the same distribution, aside from a shift in scale.
The increase in scale is a factor of $25/18 \approx 1.3889$.
Here's a boxplot of the data on the original scale and on the log-scale:
Here we do a Wilcoxon-Mann-Whitney on the log-data:
> lx=log(x)
> ly=log(y)
> wilcox.test(ly,lx,conf.int=TRUE)
Wilcoxon rank sum test with continuity correction
data: ly and lx
W = 1743, p-value = 0.1455
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-0.09044406 0.68813695
sample estimates:
difference in location
0.2917924
If we exponentiate the interval, it should contain the true ratio of scales 95% of the time. It does contain the population ratio on this occasion:
> exp(wilcox.test(ly,lx,conf.int=TRUE)$conf.int)
[1] 0.9135254 1.9900046
attr(,"conf.level")
[1] 0.95
If we exponentiate the estimated shift in the logs, it will be a (somewhat) downward-biased estimate of the scale-ratio, but it's fairly good:
> exp(wilcox.test(ly,lx,conf.int=TRUE)$estimate)
difference in location
1.338825
[The data used in the example were generated as follows:
> set.seed(94394381)
> x=rexp(60,1/18)
> y=rexp(50,1/25)
... just in case anyone wants to try it]
There are other ways to estimate the scale change; one could, for example, divide the y-sample by a constant, $c$, and try different values of $c$ until the WMW statistic was equal to its expected value under the null hypothesis. One could also try different values until the statistic was equal to the $\alpha/2$ and $1-\alpha/2$ quantiles of the null distribution of the WMW. As above, this will be the point estimate of and interval for the ratio of scales.
(That would have the advantage of not introducing a small bias into the central estimate.)
Best Answer
One can construct a confidence interval for the median difference using Walsh averages. It does assume that the distribution of the differences is symmetric, but it seems to be OK here. See this document for an explanation of the procedure. It appears that somebody even wrote a Matlab function for this calculation.
EDIT: implementation in Matlab The
psignrankfunction in R calculates the cumulative distribution of the Sign-rank test statistic. If Matlab does not have that function, then you probably can't easily calculate an exact confidence interval. However you can calculate an approximate interval as described here. I think that is the calculation attempted at the end of the Matlab code.