On the Wikipedia page for Grover's algorithm, it is mentioned that:
"Grover's algorithm can also be used for estimating the mean and median of a set of numbers"
So far I only knew how it can be used to search a database. But not sure how to implement that technique to estimate the mean and median of a set of numbers. Moreover, there's no citation (as far as I noticed) on that page which explains the technique.
The idea for estimating the mean is roughly as follows:
For any $f(x)$ that gives outputs in the reals, define a rescaled $F(x)$ that gives outputs in the range 0 to 1. We aim to estimate the mean of $F(x)$.
Define a unitary $U_a$ whose operation is $$U_a:|0\rangle|0\rangle\mapsto\frac{1}{2^{n/2}}\sum_x|x\rangle(\sqrt{1-F(x)}|0\rangle+\sqrt{F(x)}|1\rangle).$$ It is important to note that this unitary is easily implemented. You start with a Hadamard transform on the first register, perform a computation of $f(x)$ on an ancilla register, use this to implement a controlled-rotation of the second register, and then uncompute the ancilla register.
Define the unitary $G=U_a (\mathbb{I}-2|0\rangle\langle 0|\otimes |0\rangle\langle 0|)U_a^\dagger \mathbb{I}\otimes Z$.
Starting from a state $U_a|0\rangle|0\rangle$, use $G$ much like you would use the Grover iterator to estimate the number of solutions to a search problem.
The main bulk of this algorithm is amplitude amplification, as described here. The main idea is that you can define two states $$ |\psi\rangle=\frac{1}{\sqrt{\sum_x F(x)}}\sum_x\sqrt{F(x)}|x\rangle|1\rangle \qquad |\psi^\perp\rangle=\frac{1}{\sqrt{\sum_x 1-F(x)}}\sum_x\sqrt{1-F(x)}|x\rangle|0\rangle, $$ and this defines a subspace for the evolution. The initial state is $U_a|0\rangle|0\rangle=(\sqrt{\sum_x F(x)}|\psi\rangle+\sqrt{\sum_x 1-F(x)}|\psi^\perp\rangle)2^{-n/2}$. The amplitude of the $|\psi\rangle$ term clearly contains the information about the mean of $F(x)$, if we could just estimate it. You could just repeatedly prepare this state and measure the probability of getting a $|1\rangle$ on the second register, but Grover's search gives you a quadratic improvement. If you compare to the way Grover's is usually set up, the amplitude of this $|\psi\rangle$ which you can 'mark' (in this case by applying $\mathbb{I}\otimes Z$) would be $\sqrt{\frac{m}{2^n}}$ where $m$ is the number of solutions.
Incidentally, this is interesting to compare to the "power of one clean qubit", also known as DQC1. There, if you apply $U_a$ to $\frac{\mathbb{I}}{2^n}\otimes|0\rangle\langle 0|$, the probability of getting the 1 answer is just the same as the non-accelerated version, and gives you an estimate of the mean.
For the median, it can apparently be defined as the value $z$ that minimises $$ \sum_x|f(x)-f(z)|. $$ There are two steps here. The first is to realise that the function we're trying to minimise over is basically just a mean. Then the second step is to use a minimisation algorithm which can also be accelerated by a Grover search. The idea here is to use a Grover's search, and mark all items for which the function evaluation gives a value less than some threshold $T$. You can estimate the number of inputs $x$ that give $f(x)\leq T$, then repeat for a different $T$ until you localise the minimum value sufficiently.
Of course, I am skipping over some details of precise running times, error estimates etc.
