Why isn't this estimator unbiased? Suppose we have a IID sample $X_1, X_2, \cdots, X_n$ with each $X_i$ distributed as $\mathcal{N}(\mu, \sigma^2)$. Now suppose we construct (a rather peculiar) estimator for the mean $\mu$: we only choose values from the sample that are greater than a pre-decided value, say $1$, and then take the sample average of only those values:
$$\hat{\mu}=\frac{1}{n_1}\sum\limits_{X_i > 1} X_i$$
Here $n_1$ is the number of values that are greater than $1$. Now I was expecting this estimator to be highly biased. However, we have:
$$\mathbb{E}(\hat{\mu}) = \frac{1}{n_1}\sum\limits_{X_i>1}\mathbb{E}(X_i)=\frac{1}{n_1}n_1\mu=\mu$$
This simply would mean that the estimator is unbiased! But obviously if I do this by generating many numbers and choosing only those which are greater than $1$, I will never get the estimator value less than $1$ and so the estimator has to be biased. What am I missing?
 A: You're hiding random variables in a few places, then applying linearity of expectation selectively. To avoid these mistakes, rewrite your estimator as
$$
\hat{\mu}=\frac{1}{n_i}\sum_{X_i>1}X_i=\frac{\sum_{i=1}^n X_i \cdot 1[X_i>1]}{\sum_{i=1}^n 1[X_i>1]},
$$
where the random variable $1[X_i>1]$ is 1 when $X_i>1$ and 0 otherwise.
Both the numerator and the denominator are random variables. Furthermore, in the numerator, both $X_i$ and $1[X_i>1]$ are random variables, so the numerator's expectation isn't $N \cdot \mu$.
A: An intuitive way to compute this $E(\hat{\mu})$ is to consider the following alternative, but equivalent, sampling process.
Step 1: determine $n_1$ independently based on a binomial distribution (with probability $p$ equal to $pr(X>1)$).
Step 2: draw $n_1$ variables $Y_1,\dots,Y_{n_1}$ from a normal distribution truncated on the left at 1.
The variable $\bar{Y} = \frac{1}{n_1}\sum_{i=1}^{n_1} Y_i$ has the same distribution as your $\hat\mu$.

In this example we see that $E(\bar{Y}\vert n_1) = E(Y_\text{truncated})$ independent* from the sample size $n_1$ and so we also have $E(\bar{Y})= E(Y_\text{truncated})$
This means that the expectation value of your $\hat\mu$ does not equal the expectation of the normal distribution but instead the expectation value of a truncated normal distribution.

Thus the estimator is not unbiased. Your equation
$$\mathbb{E}(\hat{\mu}) = \frac{1}{n_1}\sum\limits_{X_i>1}\mathbb{E}(X_i)=\frac{1}{n_1}n_1\mu=\mu$$
Should have been
$$\mathbb{E}(\hat{\mu}) = \frac{1}{n_1}\sum\limits_{X_i>1}\mathbb{E}(X_i|X_i>1)=\frac{1}{n_1}n_1E(Y_\text{truncated}) = E(Y_\text{truncated})$$
The parameter $n_1$ is a variable, as Michael M, mentions in the comments, so it is problematic to eliminate it. But because you get the term $n_1/n_1$ you do not have this problem. The expectation of the average of $n_1$ i.i.d. variables is independent from the sample size $n_1$.

*A tricky thing is what your expression means when $n_1 = 0$. It is undefined in that case because you have a division $0/0$.
A: As indicated in my comment (and then in later answers), the error in the reasoning leading to the apparent paradox is to treat the selected or surviving $X_i$'s that we should denote differently, e.g., as $X_j^+$'s,  as distributed from the $\mathcal N(\mu,\sigma^2)$ distribution. This accept-reject selection changes the distribution of these $X^+_j$'s to a truncated $\mathcal N^+(\mu,\sigma^2,1)$ distribution, making their expectation equal to $\mathbb E[X_i|X_i>1]$. The number $n_1$ of such $X_j^+$'s is also a Bernoulli random variable (and brings information about $\mu$) but this does not impact the overall expectation (when $n_1\ne 0$).
The estimator $\hat\mu$ need be defined separately when $N_1=0$, for instance setting $\hat \mu(0)=1$ or $\hat \mu(0)=\bar X_n$. Conditional on $N_1=n_1>0$, the selected $X^+_j$'s are thus truncated Normal variates on $(1,\infty)$ with expectation
$$\mu+\dfrac{\phi(\sigma^{-1}(1-\mu))}{\Phi(\sigma^{-1}(\mu-1))}\sigma$$
Since $\mathbb P(N_1=0)=\Phi(\sigma^{-1}(1-\mu))^n$,
$$\mathbb E[\hat\mu] = \Phi(\sigma^{-1}(1-\mu))^n\hat\mu(0)+[1-\Phi(\sigma^{-1}(1-\mu))^n]\left\{\mu+\dfrac{\phi(\sigma^{-1}(1-\mu))}{\Phi(\sigma^{-1}(\mu-1))}\sigma\right\}$$
and the expectation of $1/N_1$ is not needed.
Note that, for a specific choice of $\hat\mu(0)$ like $\hat\mu(0)=1$, the expectation of $\hat\mu$ may be equal to $\mu$ for some exceptional values of $\mu$, $n$, and $\sigma$, but is biased for almost every parameter value, and definitely so when $\mu<1$. However, $\hat\mu(0)$ can be chosen for the purpose of making $\hat\mu$ unbiased, for instance by using the $X_i$'s (that are all less than $1$ when $n_1=0$).
A: Given that $n_1$ is a random variable (as pointed out already in the comments), the expected value can be computed as $E(\hat\mu)=E_{n_1}[E_{\hat \mu}(\hat\mu|n_1)]$. For the inner expectation, note that one can't just write
$E_{\hat \mu}(\hat\mu|n_1)=\frac{1}{n_1}\sum_{X_i>1}E(X_i),$ because the expected value cannot depend on specific values of certain $X_i$, as would be required for the sum. So $$E_{\hat \mu}(\hat\mu|n_1)=\frac{1}{n_1}E\left[\sum_{X_i>1} X_i|n_1\right].$$ For given $n_1$, we can write, with appropriate renumbering of indexes, $\sum_{X_i>1} X_i=\sum_{j=1}^{n_1} X_j^*$, where $X_j^*$ are random variables distributed according to a truncated normal distribution between $a=1$ and $b=\infty$. Let $E_{\mu,\sigma^2,a,b}X$ denote the expectation of such a truncated normal. For $a=1, b=\infty,$ $$E_{\mu,\sigma^2,1,\infty}X=\mu+\frac{\varphi\left(\frac{1-\mu}{\sigma}\right)}{1-\Phi\left(\frac{1-\mu}{\sigma}\right)}\sigma=t>\mu,$$ see https://en.wikipedia.org/wiki/Truncated_normal_distribution . Conditioning on $n_1$, we have $$E_{\hat \mu}(\hat\mu|n_1)=\frac{1}{n_1}\sum_{j=1}^{n_1} E(X_j^*)=\frac{1}{n_1}n_1 E_{\mu,\sigma^2,1,\infty}(X)=t>\mu.$$ This does not depend on $n_1$ (unless $n_1=0$, in which case the sum is empty and $E_{\hat \mu}(\hat\mu|n_1=0)=0$), so ultimately $$E(\hat \mu)=P\{n_1>0\}t.$$
This is $>\mu$ (bias!) if $\mu\le 0$, and also if $P\{n_1=0\}$ is small enough that $P\{n_1>0\}t>\mu$, which should hold unless $n$ is very small (potentially resulting in a large $P\{n_1=0\}$, the value of which is given in Xi'an's solution).
PS: I corrected this seeing Xi'an's solution, who got a thing right that I had forgotten about. That solution is perfectly right as far as I can see, however my different way of getting there may also help.
PPS: I take $\hat \mu=0$ in case $n_1=0$, which isn't entirely clear in the question.
