Could you explain how to derive the Pauli matrices?
$$\sigma_1 = \sigma_x = \begin{pmatrix} 0&1 \\ 1&0 \end{pmatrix}\,, \qquad \sigma_2 = \sigma_y = \begin{pmatrix} 0&-i\\ i&0 \end{pmatrix}\,, \qquad \sigma_3 = \sigma_z = \begin{pmatrix} 1&0\\0&-1 \end{pmatrix} $$
Maybe you can also link to an easy to follow tutorial ?
That certainly depends on what exactly you mean. I take your question as "how do you see that the (non-relativistic) electron spin (or more generally, Spin-1/2) is described by the Pauli matrices?"
Well, to start, we know that measuring the electron spin can only result in one of two values. From this we see that we need matrices of at least dimension 2. The simplest choice is then of course exactly dimension 2.
Moreover the spin is an angular momentum, and thus described by three operators obeying the angular momentum algebra: $[L_i, L_j] = \mathrm{i} \hbar\epsilon_{ijk} L_k$. This together with matrix dimension 2 basically restricts the choice to sets of three matrices which are equivalent to $\hbar/2$ times the Pauli matrices (the freedom of choice of those matrices corresponds to the freedom to use three arbitrary orthogonal directions as $x$, $y$ and $z$ direction).
So now why choose from those equivalent choices exactly the Pauli matrices? Well, there's always one measurement direction which is represented by a diagonal matrix; this makes calculations much easier. Of course it makes sense to choose the matrices in a way that this direction is one of the coordinate directions. By convention, the $z$ direction is chosen. This ultimately fixes the matrices.
Spin is an angular momentum, so in the rest frame it is a 3-dimensional vector, or 4-dimensional vector with zero time component:
$\vec{v} = (v_1,v_2,v_3)$
Each 3D vector can be associated with a 2x2 matrix by the following rule:
$V = \begin{vmatrix} v_3 & v_1-iv_2 \\ v_1+iv_2 & -v_3 \end{vmatrix}$
In particular, if you chosse $\vec{v}$ as a basis vector: $\vec{v} = (1,0,0)$, it is associated with matrix
$H_1=\begin{vmatrix} 0 & 1 \\ 1 & 0 \end{vmatrix}$
Similarly, for $\vec{v} = (0,1,0)$ you will get
$H_2=\begin{vmatrix} 0 & -i \\ i & 0 \end{vmatrix}$
and for $\vec{v} = (0,0,1)$ you will get
$H_3=\begin{vmatrix} 1 & 0 \\ 0 & -1 \end{vmatrix}$
These are Pauli matrices. Now arbitrary vector corresponds to linear combination of $H_1, H_2$ and $H_3$. For given $\vec{v} = (v_1,v_2,v_3)$ you will get
$V = v_1 H_1 + v_2 H_2 + v_3 H_3$
Further you can use matrices to act on special objects with 2 complex components. These objects are called "spinors". They are used to construct wave functions of fermions with 1/2 spin. For instance, we can choose spinor of the form
$s_3=\begin{vmatrix} 1 \\ 0 \end{vmatrix}$
and act on it with matrix $H_3$. We will obtain:
$H_3 s_3=\begin{vmatrix} 1 & 0 \\ 0 & -1 \end{vmatrix} \begin{vmatrix} 1 \\ 0 \end{vmatrix} = \begin{vmatrix} 1 \\ 0 \end{vmatrix} = s_3$
As you can see, $s_3$ is an "eigenspinor" of $H_3$ with eigenvalue +1. There is also another "eigenspinor" of $H_3$ with eigenvalue -1:
$\begin{vmatrix} 1 & 0 \\ 0 & -1 \end{vmatrix} \begin{vmatrix} 0 \\ 1 \end{vmatrix} = -\begin{vmatrix} 0 \\ 1 \end{vmatrix}$
Why do we use spinors instead of scalars, vectors and tensors? Because:
Now if, for instance, $s$ is an "eigenspinor" of $H_3$, then direction of the vector constructed from $s$ will be invariant w.r.t. rotations around z axis.
The Hilbert space for spin 1/2 is two-dimensional - there are two possible values spin can take: $\hbar/2$ or $-\hbar/2$ (this is taken from experiment).
Now, in two-dimensional Hilbert space spin operator has to be self-adjoined (this comes from foundations of QM). Furthermore, sum of its eigenvalues has to be 0 - because sum of eigenvalues is just a sum of possible results of measurements - in this case $\hbar/2-\hbar/2 = 0$. Therefore the operator corresponding to the measurement of spin in a given direction has to be 2x2 complex hermitan traceless matrix (no commutation relations used so far!). This family of matrices is a three-parameter one (one constant on diagonal, two off-diagonal for real and imaginary part, rest is determined by hermicity and tracelessness).
Furthermore, every such an operator has to be diagonal in its eigenbasis - which corresponds to the measurement in a given direction. Let us name this direction $z$ direction. We see, that the only possible traceless hermitean diagonal matrix is a multiplicity of Pauli matrix $\sigma_z$. Now, we write the remaining two operators as $L_x$ and $L_y$.
Next, we have to use commutation relations of angular momentum: $[L_y,\sigma_z]=L_x$, $[L_x,\sigma_z]=-L_y$ and $[L_x, L_y] = \sigma_z$. Remembering that $L_y$ as well as $L_x$ are determined by 3 real constants each, we need 6 linear equations to solve for them all at once. The above commutation relations give exactly 6 equations (3 matrix equations, two-dimensional matrices) so we can solve for $L_x$ and $L_y$ which are exactly the remaining Pauli matrices.
Yet another way is from almost purely experimental considerations. Consider a Stern-Gerlach set up; the identity is $$\mathbb{I}=|\uparrow\rangle \langle \uparrow |+|\downarrow\rangle \langle \downarrow | $$ Now experimentally we can only get two different numbers, $\hbar/2$ and $-\hbar/2$, and regardless of how we are oriented, in the natural basis corresponding to our orientation the spin operator must take the form $$S=\frac{\hbar}{2}(|\uparrow\rangle \langle \uparrow |-|\downarrow\rangle \langle \downarrow |$$ So now lets fix a basis, say the z, and hit the identity from one side with a spin vector in a different basis, say the x $$|\uparrow\rangle_{x}=|\uparrow\rangle\langle \uparrow |\uparrow\rangle_{x}+|\downarrow\rangle \langle \downarrow |\uparrow\rangle_{x}$$ Through various arguments, the two inner-products can be evaluated to give $\pm1/\sqrt{2}$ to get $$|\uparrow\rangle_{x}=\frac{1}{\sqrt{2}}|\uparrow\rangle-\frac{1}{\sqrt{2}}|\downarrow\rangle $$ to see the arguments see Ballentine, or Sakurai, or Schwinger. Now to get $S_x$ we use the expression for $S$ from above $$S_x=\frac{\hbar}{2}[(|\uparrow\rangle \langle \uparrow |)_x-(|\downarrow\rangle \langle \downarrow |)_x]$$ where my subscript x means in the x basis, but we know how to express those vectors in the "known" basis, the z basis. Substitution yeilds $$S_x=\frac{\hbar}{2}(|\uparrow\rangle\langle\downarrow |+|\downarrow\rangle\langle\uparrow |)$$ Consider the following equal signs, begotten with entirely using the identity symbol $$S_x=\sum_{n=1}^{2}\sum_{m=1}^{2}|n\rangle\langle n|S_x|m\rangle\langle m|$$ The $n$ and $m$ run from 1 to 2 standing for spin up and down. Now evaluate the guy in the middle. this is a two by two matrix whos elements are $$\langle \uparrow |S_x|\uparrow\rangle\quad\langle \uparrow|S_x|\downarrow\rangle\quad\langle\downarrow|S_x|\uparrow\rangle\quad\langle\downarrow|S_x|\downarrow\rangle$$ to see what these equal just hit the $S_x$ buy above with the two up/down vectors in the z basis and use the orthogonality condition to make them either 1 or zero. The result is one of the Pauli matrices.
Hope this helps.
This link seems to be along the required path of thought.
Please note the "axiomatic" facts: experimental inputs in value of $S^2$, raising and lowering operators, desirability of hermitian operators... that go inside the derivation.
Also, once having chosen them, note that the 3 Pauli matrices along with the 2d identity matrix can be used as a basis to write down any 2x2 matrix as has been mentioned here.
