Convergence of Bernoulli sampling procedure 
*

*Let $X$ be a variable which has density $f_{X}$ and mean $E[X]=\mu$

*Define $B$, a Bernoulli variable, which has a probability $P(b_{i}=1)=p_{i}=\frac{1}{N^\gamma}$


Consider the following expression:
\begin{equation}
\frac{\sum_{i=1}^{N}b_{i}x_{i}}{\sum_{i=1}^{N}b_{i}}
\end{equation}
This can be interpreted as a sampling procedure where the probability of being included in the sample is a Bernoulli draw. The above expression is therefore simply an average.
Is the following derivation correct:
\begin{equation}
\lim_{N \rightarrow \infty}\frac{\frac{1}{N}\sum_{i=1}^{N}b_{i}x_{i}}{\frac{1}{N}\sum_{i=1}^{N}b_{i}}=\frac{E[BX]}{E[B]}=\frac{E[E[BX|X]]}{E[B]}=\frac{E[\frac{1}{N^{\gamma}}X]}{E[\frac{1}{N^{\gamma}}]}=E[X]
\end{equation}
I am unsure whether this is correct because of the way the Bernouilli probability is specified: $p_{i}$ tends to zero as $N \rightarrow \infty$. But I do not see where this would invalidate the above derivation - although I have a feeling it might.
Any help is much appreciated.
 A: A not irrelevant digression, to validate also my comment about $\gamma$: We are considering an expression where in the numerator we have the sum of  Bernoulli random variables, i.e. dichotomous discrete random variables taking the values $\{0,1\}$. So we have to consider whether the event "sum of Bernoullis is zero" has strictly positive probability.  
The probability that all draws will be zero is equal to 
$$\Pr \left(\text{\{All b_i's are zero\}}\right) = \left(1-\frac {1}{N^{\gamma}}\right)^N$$
Then for example, if $\gamma =1$ we have
$$\lim_{n\rightarrow \infty}\Pr \left(\text{\{All b_i's are zero\}}\right) = \lim_{n\rightarrow \infty}\left(1-\frac {1}{N}\right)^N = 1/e \approx 0.37$$
and a very high probability indeed. Informally you can check that if $\gamma >1$ the sum will have limiting probability of being zero equal to unity, while only if $\gamma<1$ the limiting probability of the sum being zero, will equal zero.
Another way to explore this is to remember that the sum of i.i.d Bernoullis is a Binomial random variable,
$$\sum_{i=1}^{N}b_{i} \sim \text{Bin}\big(\mu_n = N^{1-\gamma}, \sigma^2_n = N^{1-\gamma}-N^{1-2\gamma}\big)$$
If $\gamma <1$, as a comment indicated that it is the case of interest here,  we avoid the case where the sum in the denominator takes the value zero, at the limit (because for finite $N$ this can still very well be the case).
Assuming from here on that $\sum_{i=1}^{N}b_{i} \neq 0$, we can define
$$w_i = \frac {b_i}{\sum_{i=1}^{N}b_{i}},\;\; \sum_i^N w_i = 1,\;\; \forall N$$ 
and the expression of interest becomes
$$\sum_{i=1}^N w_ix_i$$
which is a weighted average, and more importantly, a convex combination of the $x_i$'s, but with the weights being random variables. This is a rather advanced technical issue, that nevertheless has been studied. A bit informally, and exploiting the fact that the weights relate to Bernoullis, note that the values each weight can take are two,
$$w_i \in \left\{0, \frac 1{\sum_{i=1}^Nb_i}\right\}$$
We wonder whether
$$\sum_{i=1}^N w_ix_i - \frac 1N\sum_{i=1}^N x_i \xrightarrow{p}0 \;\;???$$
Denote $M$ the number of Bernoullis that take the value $1$. This means that $\sum_{i=1}^N b_i =M$.  
Then we can write
$$\sum_{i=1}^N w_ix_i - \frac 1N\sum_{i=1}^N x_i = \sum_{w_i \neq 0}w_ix_i - \frac 1N\sum_{i=1}^N x_i $$
$$=\frac {1}{\sum_{i=1}^N b_i}\sum_{w_i \neq 0}x_i - \frac 1N\sum_{i=1}^N x_i$$
$$=\frac {1}{M}\sum_{w_i \neq 0}x_i - \frac 1N\sum_{i=1}^N x_i$$
Now $M$ is not a number but the Binomial random variable we described before. Still we have that $M\leq N$, and if $M\rightarrow \infty$ the two sums above will have the same probability limit, $\mu$, and the assertion of the question is verified. Intuitively, if, as $N$ becomes "very large", $M$ will tend to take very large values also, the first sum will behave like a proper sample average of the $x$'s.
So the condition we require is that
$$M=\sum_{i=1}^Nb_i \rightarrow \infty$$ 
By looking at its mean and variance, we see that as long as $\gamma <1$ it will diverge. So it appears that 
$$\gamma <1 \implies \frac{\sum_{i=1}^{N}b_{i}x_{i}}{\sum_{i=1}^{N}b_{i}} \xrightarrow{p} \mu$$
A quick simulation supports this: I generated $50,000$ samples from a chi-square with $3$ degrees of freedom (so its expected value is $3$), and the same number of samples from a Bernoulli with $p = (50,000)^{-1/2}$ (so $\gamma = 1/2$ here). I then calculated the expression of interest and I obtained
$$\frac{\sum_{i=1}^{N}b_{i}x_{i}}{\sum_{i=1}^{N}b_{i}} = 3.016$$
