Why can we replace a cavity inside a sphere by a negative density? I have a sphere with radius $R$ and inside this sphere there is a smaller sphere with radius $\frac{R}{3}$. This small cavity has its center at $\frac{R}{2}$, it doesn't matter in which direction.
If I want to find the electric field at the point $(0.0.0)$, why can I not enclose it with a surface and apply Gauss' law, for example with a small little sphere around $(0,0,0)$ but not touching the cavity. Is it because the sphere is not symmetric outside the Gaussian surface?
I find lots of solutions on the internet that say you can replace the cavity with a negative density, why?
http://jkwiens.com/2007/10/24/answer-electric-field-of-a-nonconducting-sphere-with-a-spherical-cavity/
 A: "I find lots of solutions on the internet that say you can replace the cavity with a negative density, why?"
Because they use a trick to calculate the potential easier. They assume that the empty hole is neutral, but composed of a positive charge density equal to that of the sphere plus a negative charge density of the same amount. In this way you can compute the potential of a sphere with the hollow sphere inside as the potential of a uniform large sphere and that of a negatively charged sphere located at the hole.
A: The principle in general is called superposition in physics or linearity in mathematics. It is very useful when you want to study a system for that system to be approximately linear.
What linearity is, in a more general context.
Here "linear" is the property of a function (or an operator, or whatever) to distribute over addition. So for example if you have a 3-dimensional vector space, a linear transform of the points of that space would be a function $f$ such that $f(\vec a + \vec b) = f(\vec a) + f(\vec b)$, usually paired with the property that also if $k$ is a scalar from the scalar field associated with the vector field, then $f(k~\vec v) = k~f(\vec v).$ It turns out that an easy way to think about all of these transformations is matrices. The idea is that we find some basis vectors $\hat e_1,\;\hat e_2\;\hat e_3$ of the space and represent every vector with components $\vec v = v_1~\hat e_1 + v_2~\hat e_2 + v_3~\hat e_3.$ This includes the vectors $f(\hat e_1),\;f(\hat e_2),\;f(\hat e_3),$ so we write these as a matrix $$f(\vec v) = v_1~f(\hat e_1) + v_2~f(\hat e_2) + v_3~f(\hat e_3) = \mathbf{F}~\vec v = \Big[f(\hat e_1)\;f(\hat e_2)\;f(\hat e_3)\Big]\begin{bmatrix}v_1\\v_2\\v_3\end{bmatrix}. $$ If we write the $f(\hat e_i)$ components vertically and number the overall matrix in the obvious way we get $f(\vec e_j) = \sum_i \hat e_i ~ F_{ij}$ which causes the above matrix product to look like $\vec u = \mathbf{F}~\vec v ~\leftrightarrow~ u_i = \sum_j F_{ij}~v_j,$ so it's naturally a sum over the adjacent indices.
How it applies in superposition cases
The fields caused by a charge distribution $\rho$ and current distribution $\vec J$ are the Maxwell equations, which take the form:$$\begin{matrix}\nabla\cdot\vec E&=& \kappa \rho &\;\;& \nabla\times\vec E&=&-(\kappa/\mu c^2) \frac{\partial \vec B}{\partial t}\\
\nabla\cdot\vec B&=& 0 &\;\;& \nabla\times\vec B&=&\mu \vec J + (\mu/\kappa)\frac{\partial \vec E}{\partial t}
\end{matrix}$$for parameters $\mu,\kappa$ which depend on your unit system of choice: for SI units, $\mu = \mu_0$ and $\kappa = 1/\epsilon_0;$ for Gaussian units, $\mu = 4\pi/c$ and $\kappa = 4\pi.$ 
Whatever the units, these equations are linear too -- that means that if I have the solutions $\vec E_1, \vec B_1$ to inputs $\rho_1, \vec J_1$ (with some boundary conditions) and $\vec E_2, \vec B_2$ respectively for $\rho_2, \vec J_2,$ then a valid solution for inputs $\rho_1 + \rho_2, \vec J_1 + \vec J_2$ is the fields $\vec E, \vec B$ where, you guessed it, $\vec E = \vec E_1 + \vec E_2$ and likewise for $\vec B.$ It goes from being a solution to being the solution if those fields combine to have the right total boundary conditions.
In addition if $\vec B$ is not changing, then $E = -\nabla V$ for some "electric potential" (or "voltage field") $V$, and then $V = V_1 + V_2$ is a correct solution for the voltage field.
The many ways this applies to your context
The Coulomb field due to a point charge $q$ in vacuum (with the field vanishing at infinity) is $\vec E = \kappa~q~\hat r/(4\pi~r^2)$ where $r$ is the distance to the particle; this becomes the potential $V = - \kappa q / (4\pi~r).$
So if we try to calculate the potential at a position $(0, 0, z)$ due to a uniform shell of surface charge $\sigma$ and radius $R$, in spherical coordinates the charge is $dq = \sigma~dA = R^2 \sin\phi ~d\phi~d\theta$ where the position $(x, y, z)$ on the sphere is $(R~\sin\phi~\cos\theta, R~\sin\phi~\sin\theta, R~\cos\phi).$ The difference in distances between these two points is $r = \sqrt{R^2\sin^2\phi + (R \cos \phi - z)^2}$ which you can expand as: $$V = \int dV = \int \frac{dq~\kappa}{4\pi~r} = \int_0^{2\pi} d\theta~\int_0^\pi d\phi~\frac{\kappa\sigma~\sin\phi}{4\pi~\sqrt{R^2 + z^2 - 2 R z \cos\phi}}.$$The integral over $\theta$ is just $2\pi;$ the integral over $\phi$ is done by switching to $u = R^2 + z^2 - 2 R z \cos \phi,\; du = 2 R z \sin\phi~d\phi,$ yielding $$V = \frac{\kappa\sigma}2 \int_{(R-z)^2}^{(R+z)^2} \frac{du}{2Rz} ~u^{-1/2} = \frac{\kappa\sigma}{2Rz} \big[|R + z| - |R - z|\big].$$For $z < R$ the signs on both terms are positive and they cancel out to simply $2z$, so the potential inside is $\kappa\sigma/R$, totally independent of $z.$ Such a constant voltage means that there is zero electric field inside the shell. Similarly, for $z > R$ we have $|R + z| = R + z$ while $|R - z| = z - R$ and this cancels out to $\kappa\sigma/z,$ decaying with the classic $1/r$ dependence of a point charge. (In fact if you look closer, or with a Gaussian surface, you see that the shell has the exact external field as if all of its charge is focused at the center of the sphere, but zero internal field.) 
So we just calculated the field inside the spherical shell of charge as zero by superposition (that's what $\int dq~\dots$ is). But also now we can use superposition again to deal with spherically-symmetric charge distributions, and then again to calculate the field inside the cavity.
As bruce smitherson answered in his much shorter response, it's as simple as filling a sphere with charge density $+\rho$ and then superimposing on it a smaller off-center sphere of charge density $-\rho$. When these are superimposed, the net charge density in their overlap is 0 while the net charge density everywhere else in the larger sphere is a positive constant.
