Are we guaranteed that the harmonic series minus infinite random terms always converge?

Consider the known harmonic series $\sum_{n=1}^\infty \frac{1}{n}$ and modify it as follows

$$\sum_{n=1}^\infty a_n\frac{1}{n}$$ where $$a_n \sim \operatorname{Bern} \left({\frac{1}{2}}\right)$$ i.e. each $a_n$ is $0$ or $1$ with probability $\frac{1}{2}$ (basically what one is doing here is removing randomly an infinite number of terms of this series)

Is it possible to know that this series is almost always convergent? Meaning that it converges with probability $1$?

For any such sequence $(a_n)$, the corresponding sequence $(b_n)$ given by $b_n=1-a_n$ is just as likely. So according to your intuition, the series $$\sum_{n=1}^\infty b_n\frac{1}{n}$$ should converge almost surely. But then $$\sum_{n=1}^\infty \frac{1}{n} = \sum_{n=1}^\infty a_n\frac{1}{n} + \sum_{n=1}^\infty b_n\frac{1}{n}$$ would converge too!

Your series is almost surely divergent (i.e., it diverges with probability $1$). Let $X_n = a_n/n$ be the $n$th term. Then $$ E[X_n]=\frac{1}{2n} $$ and $$ {\text{Var}}[X_n]=E[X_n^2]-E[X_n]^2=\frac{1}{4n^2}. $$ The variance of the partial sums is therefore $$ {\text{Var}}\left[\sum_{i=1}^{n}X_i\right]=\frac{1}{4}\sum_{i=1}^{n}\frac{1}{i^2}\rightarrow\frac{\pi^2}{24}; $$ the expectation value of the partial sums, on the other hand, is $$ E\left[\sum_{i=1}^{n}X_i\right]=\frac{1}{2}\sum_{i=1}^{n}\frac{1}{i}=\frac{1}{2}H_n \sim \frac{1}{2}\log n. $$ So the sequence of partial sums diverges logarithmically with probability $1$.

As I'm getting some down votes, I wanted to add a few steps to show that this argument is sound. I've got random variables $S_n$ such that $\mu_n = E[S_n]\rightarrow\infty$ and $\sigma_n^2={\text{Var}}[S_n]\rightarrow \sigma^2 < \infty$. I claim that $S_n$ is unbounded with probability $1$. If $S_n$ is bounded, then for some $x\in\mathbb{R}$ and some $N\in\mathbb{N}$, $S_n \le x$ for all $n\ge N$. In particular, for any sequence of indices $n_i\rightarrow\infty$ and sequence of bounds $A_i\rightarrow\infty$, $$ S_n {\text { bounded}} \implies S_{n_1} < A_1 \vee S_{n_2} < A_2 \vee \ldots, $$ and hence $$ {\text{Pr}}[S_n{\text { bounded}}] \le {\text{Pr}}[S_{n_1}<A_1] + {\text{Pr}}[S_{n_2}<A_2] + \ldots. $$ Let $M>0$ be a fixed, large number. Because $\mu_n\rightarrow\infty$ and $\sigma^2_n\rightarrow\sigma^2$, we can choose $n_i$ large enough that $\mu_{n_i} > 2^{i/2+1} M \sigma_{n_i}$ and $n_i\rightarrow\infty$; and choose $A_i=\mu_{n_i} - 2^{i/2} M\sigma_{n_{i-1}}$. But $$ {\text{Pr}}[S_{n_i} < \mu_{n_i} - 2^{i/2} M\sigma_{n_{i-1}}] \le {\text{Pr}}\left[\frac{(S_{n_i}-\mu_{n_i})^2}{\sigma^2_{n_{i-1}}} > 2^{i} M^2\right] \le \frac{1}{2^{i}M^2}, $$ so $$ {\text{Pr}}[S_n{\text { bounded}}] \le \sum_{i=1}^{\infty}\frac{1}{2^{i}M^2} = \frac{1}{M^2}. $$ Since $M$ was arbitrary, we can choose it to be as large as we like; we conclude that ${\text{Pr}}[S_n{\text { bounded}}]=0$, and hence that $S_n$ converges with probability $0$ as well. In our case, we need only the additional fact that $S_n$ is monotonically increasing to conclude that $S_n$ a.s. diverges to infinity (since that's the only other option).

Why would it ever converge ?

Taking the expectation,

$$E\left(\sum_{n=1}^\infty a_n\frac{1}{n}\right)=\sum_{n=1}^\infty E(a_n)\frac{1}{n}=\frac12\sum_{n=1}^\infty \frac{1}{n}.$$

Even if the argument isn't rigorous because none of these series converge, it should be enough to shed doubt.

In the harmonic series, removing all but every millionth term is not enough to restore convergence, because $\sum\frac1{1000000n}=\frac1{1000000}\sum\frac1n$.

This answer combines and elaborates on two comments to the OP's post. Let

$$a_n \sim \operatorname{Bern} \left({f(n)}\right),\;\;\; f(n) \in (0,1)$$

set $X_n \equiv n^{-1}a_n$ and consider

$$\sum_{i=1}^\infty X_i $$

with the random variables assumed independent.

Checking the conditions of Kolmogorov's Three-series Theorem, set $A > 1$. Then

$$\Pr(|X_n|\geq A)=0 \implies \sum_{i=1}^{\infty} \Pr(|X_i|\geq A) =0 \tag{1}$$

converges (this is just a special case that holds in general for any random variable with support bounded from above).

Also $I_{\{|X_n|\leq A\}} =1\;\; \forall\, n$, so

$$\sum_{i=1}^{\infty}E[X_i\cdot I_{\{|X_i|\leq A\}}] = \sum_{i=1}^{\infty}E[X_i] =\sum_{i=1}^{\infty}\frac {f(i)}{i} \tag{2}$$

A choice for $f(n)$ for this to converge is $f(n) = n^{-\alpha},\, \alpha >0$, since this will give a hyperharmoninc series, which converges.

Finally,

$$\sum_{i=1}^{\infty}\text{Var}[X_i\cdot I_{\{|X_i|\leq A\}}] = \sum_{i=1}^{\infty}\text{Var}[X_i] =\sum_{i=1}^{\infty}\frac {f(i)[1-f(i)]}{i^2} \tag {3}$$

which converges for all and any $f(i) \in (0,1)$, since it is strictly smaller than $\sum_{i=1}^{\infty}i^{-2} = \pi^2/6$.

So the critical condition, is the convergence of the sum of the expected values (condition $2$), since it is this condition that imposes restrictions on the form of $f(n)$ in order to obtain convergence.

This post directly proves that

$$ P\left( \sum_{n=1}^{\infty} \frac{a_n}{n} = \infty \right) = 1$$

rather than study about different quantities so as to make a heuristic argument that helps convince that the result is plausible.

The method is

• Decompose the sum into a sequence of independent, disjoint partial sums with rapidly decreasing variance
• Find a lower bound on the probability that each partial sum is greater than $0.01$
• Find a lower bound on the probability that infinitely many of the partial sums are all greater than $0.01$ (and thus the sum diverges)
• Deduce from the form of the final result that the probability of divergence is $1$

Borrowing the methods of the other answer,

Let $X_n = a_n / n$ be the $n$-th term.

Define a new sequence of independent random variables

$$ S_n = \sum_{i=2^n}^{2^{n+1} - 1} X_n $$

so that the total sum is

$$ \sum_{i=1}^{\infty} \frac{a_i}{i} = \sum_{n=0}^{\infty} S_n $$

We have

$$ E[S_n] = \sum_{i=2^n}^{2^{n+1} - 1} E[X_n ] = \frac{1}{2} (H_{2^{n+1} - 1} - H_{2^n - 1} ) \approx \frac{1}{2} \ln(2)$$

where $H_n$ is the $n$-th harmonic number. We also have the variance

$$ Var[S_n] = \sum_{i=2^n}^{2^{n+1} - 1} \frac{1}{4i^2} < 2^{-n} $$

(we can prove a slightly stronger bound, but I'm lazy!)

Chebyshev's inequality then says that, for any positive $k$,

$$ P\left(\left|S_n - \frac{1}{2} \ln 2\right| < k 2^{-n/2} \right) \geq 1 - \frac{1}{k^2} $$

For sufficiently large $k$, I claim that

$$ \exp(-2/k^2) = 1 - \frac{2}{k^2} + \frac{2}{k^4} - \ldots \leq 1 - \frac{1}{k^2} $$

and therefore,

$$ P\left(\left|S_n - \frac{1}{2} \ln 2\right| < \frac{1}{4} \right) \geq 1 - \frac{1}{2^{n-4}} \geq \exp\left(-\frac{1}{2^{n-5}} \right) $$

These events are jointly independent, so we have

$$ P\left(\forall n \geq m: \left|S_n - \frac{1}{2} \ln 2\right| < \frac{1}{4} \right) \geq \exp\left(-\sum_{n = m}^{\infty} \frac{1}{2^{n-5}} \right) = \exp(-2^{6-m}) $$

In particular, this means

$$ P(\forall n \geq m: S_n > 0.01) \geq \exp(-2^{6-m}) $$

from which it follows

$$ P\left( \sum_{n=0}^{\infty} S_n = \infty \right) \geq \exp(-2^{6-m}) $$

by taking the limit as $m \to \infty$, we get

$$ P\left( \sum_{n=0}^{\infty} S_n = \infty \right) \geq 1 $$