Let's clarify some of the issues.

###Convexity

Consider a real-valued function function $f$ defined on an interval $I$ of real numbers. (Everything that follows extends naturally and obviously to functions defined on sets $I\subset \mathbb{R}^d$ provided $I$ is a convex set, but I won't belabor that point. The key property enjoyed by $I$ is that "mixtures" of elements of $I$ are in $I;$ this will be elaborated below.)

$f$ has a graph (the set of points $(x,f(x))$ where $x\in I$). The convexity of $f$ (in the sense given in the question) is equivalent to the graph always "curving upwards." This concept can be made rigorous in various ways; one nice one is that all points on any line segment connecting any two points on the graph of $f$ lie on or above the graph. This is just a geometric restatement of convexity.

One consequence of convexity is this: given any finite number $n$ of points $x_i,$ $i=1,2,\ldots, n$ in $I$ and any set of non-negative numbers ("weights") $w_i$ summing to unity, $$\sum_{i=1}^n w_i\,f(x_i) \ge f\left(\sum_{i=1}^n w_i\,x_i\right).\tag{1}$$

(The right hand side is defined because the conditions on the $w_i$ imply the argument of $f$ lies inside $I,$ allowing us to apply $f$ to it. For brevity, let's call linear combinations where the coefficients are non-negative and sum to unity "mixtures" of the values.)

Proof: Use induction on $n,$ beginning with the "base case" $n=2,$ which is just the defining property of $f.$ Assume hypothetically that $(1)$ holds for a given $n\ge 2$ and all applicable $x_i$ and $w_i.$ Let $x_i$ and $w_i,$ $i=1,2,\ldots, n+1,$ be $n+1$ points in $I$ with $n+1$ weights, respectively. If the last weight $w_{n+1}$ is zero then (1) applies immediately to the other $n$ points. It remains to treat the case where $w_{n+1}\ne 0.$ This enables us to focus on the first $n$ weights and renormalize them to sum to unity upon division by $1-w_{n+1}:$

$$\eqalign{ \sum_{i=1}^{n+1} w_i\,f(x_i) &= \sum_{i=1}^{n} w_i\,f(x_i) + w_{n+1}\,f(x_{n+1})\\ &= (1-w_{n+1})\sum_{i=1}^{n} \frac{w_i}{1-w_{n+1}}\,f(x_i) + w_{n+1}\,f(x_{n+1})\\ &\ge (1-w_{n+1})f\left(\sum_{i=1}^{n} \frac{w_i}{1-w_{n+1}}x_i\right) + w_{n+1}\,f(x_{n+1}) }$$

by the inductive hypothesis. But this last expression is a mixture of just two values of $f.$ The convexity of $f$ immediately implies

$$\eqalign{ &(1-w_{n+1})f\left(\sum_{i=1}^{n} \frac{w_i}{1-w_{n+1}}x_i\right) + w_{n+1}\,f(x_{n+1})\\ &\ge f\left(\sum_{i=1}^{n} (1-w_{n+1})\frac{w_i}{1-w_{n+1}}x_i + w_{n+1}\,x_{n+1}\right)\\ &= f\left(\sum_{i=1}^{n+1} w_i\,x_i\right). }$$

From the assumption that $(1)$ holds for a particular $n\ge 2$ we have deduced $(1)$ also holds for $n+1.$ Thus, $(1)$ holds for all $n\ge 2,$ QED.

###Indexing random variables

Recall that a random variable $X$ is function defined on a set of outcomes $\Omega.$ Let $f$ be the function in the question and let $I$ be its domain. For expressions like "$f(X)$" to make sense we must assume the image of $X$ is a subset of $I,$ so that we may define the random variable $f(X)$ to be the function on $\Omega$ given by

$$f(X)(\omega) = f(X(\omega))$$

for all $\omega\in\Omega.$ (This is the composition of the function $f$ with the function $X.$)

Consider $n$ random variables defined on a common set $\Omega.$ Two hundred years ago, expositors would have described them by giving with an indefinite sequence of letters, like this: "Let $A,B,C,\ldots, Z$ be a finite collection of random variables." Since then, mathematical terminology has evolved toward using indexing to give such descriptions. An "index" is merely an element of a specified set $S,$ often taken to be $S=\{1,2,\ldots, n\}$ or $S=\{0,1,\ldots, n-1\}.$ Our "$n$ random variables" are then described by supposing we are given a function

$$X: S \to \mathcal{R}(\Omega)\tag{2}$$

where $\mathcal{R}(\Omega)$ is the set of random variables defined on $\Omega.$ When we are indexing and $i\in S$ we don't write $X(i)$ for the random variable associated with $i:$ instead, we write $X_i.$ This is merely a matter of notation; it is not conceptual. Thus,

For each $i\in S,$ $X_i:\Omega\to I$ is a random variable (and these variables are not necessarily distinct).

Often the entire indexed collection of objects is denoted as $(X_i)_{i\in S}.$ This is just another way of writing $(2).$

###Convexity applied to random variables

Let $f$ be the function in the question, defined on an interval $I$ and assumed to be convex (and measurable), and let $(X_i)_{i\in S}$ be a finite collection of random variables. Composing $f$ with each $X_i$ produces a random variable $f(X_i):\Omega\to \mathbb{R}.$

Consider any particular outcome $\omega\in\Omega$ and write $x_i = X_i(\omega)$ for each $i\in S.$ We have seen that $f(x_i)=f(X_i(\omega))$ can be considered the value of $f(X_i)$ at $\omega.$ Because $f$ is convex and $x_i\in I,$ for any mixture with coefficients $w_i$ the inequality $(1)$ says

$$\sum_{i\in S} w_i\, f(X_i)(\omega) = \sum_{i\in S} w_i\, f(x_i) \ge f\left(\sum_{i\in S} w_i\,x_i\right) = f\left(\sum_{i\in S} w_i\,X_i(\omega)\right).$$

Take expectations of both sides with respect to the probability measure defined on $\Omega.$ Because expectation respects inequalities (which comes down to the basic and obvious fact that the integral of a non-negative function is a non-negative number), we obtain

$$\mathbb{E}\left[\sum_{i\in S} w_i\, f(X_i)\right] \ge \mathbb{E}\left[f\left(\sum_{i\in S} w_i\,X_i\right)\right].$$

Linearity of expectation lets us re-express the left hand side as

$$\sum_{i\in S} w_i\, \mathbb{E}\left[f(X_i)\right] = \mathbb{E}\left[\sum_{i\in S} w_i\, f(X_i)\right]\ge \mathbb{E}\left[f\left(\sum_{i\in S} w_i\,X_i\right)\right],$$

giving the assertion in the question.