Deriving the formula for the Möbius function? There is a rich and beautiful theory in combinatorics which deals with the Möbius functions of posets and arrives at the "Classical" formula of the Möbius function. Let's forget all about it, since I'm hoping to find an elementary solution.
I wish to arrive at the formula for the classical Möbius function $\mu$ (which is defined as $1$ on $1$, $0$ on non square-free numbers, and as $(-1)^k$ for square-free numbers with exactly $k$ prime divisors) from the following formula:
$$\sum_{d|n}\mu(d)=\cases{1, & n=1;\\ 0, & n>1.}$$
So let's assume I'm given this formula and nothing else (especially not a pre-hand knowledge of the definition of $\mu$) - is there a simple way to derive the standard definition, without using deeper theory or "insight"?
 A: Even though the question demands forgetting about posets, I'll summarize the general approach of Möbius inversion here, since it is no less elementary than a "direct" approach, does not rely on deep insight, and is in fact quite painless.
Suppose $(P,\preceq))$ are a set $P$ and a partial ordering $\preceq$ on $P$ (for instance $P=\mathbb N_{\geq1}$, and $x\preceq y$ means $x\mid y$), and suppose that for each $y\in P$ there are only finitely many $x\in P$ with $x\preceq y$, so that our summations below will make sense. We want a function $\mu(x,y)$ on $P^2$ with the property
$$
  \sum_{y\preceq z}\mu(x,y)=\delta_{x,z}.
$$
for all $x,z\in P$. This equation defines $\mu(x,z)$ once all other $\mu(x,y)$ that occur are known, so the values $\mu(x,y)$ are uniquely defined by this equation, as can be shown for fixed $x$ by strong induction on $y$ along the partial ordering $\preceq$ (or along any total ordering extending it, if one prefers). This shows in particular that $\mu(x,y)=0$ whenever $x\not\preceq y$, and $\mu(x,x)=1$ for all $x$.
In our example we will get $\mu(x,y)=\mu(y/x)$ with $\mu$ now denoting the classical (arithmetic) Möbius function. Indeed for $x=1$ the equation becomes the usual $\sum_{y\mid z}\mu(1,y)=\delta_{1,z}$, and for $x>1$ one can easily check that $\mu(x,y)=\mu(1,y/x)$ for all multiples $y$ of $x$. It is somewhat deranging that the general setting requires a function of two arguments rather than one, but that setting provides no means to express anything like the quotient operation used to reduce to a single argument; nevertheless in concrete situations there is often a similar simple function of $x$ and $y$ in terms of which $\mu(x,y)$ is defined. Note that the equation is asking $\mu(x,y)$ to be the coefficients of the inverse of the matrix $M$, with rows and columns indexed by elements of $P$, given by $M_{x,y}=1$ if $x\preceq y$ and $M_{x,y}=0$ otherwise. Since a left inverse is also right inverse, this point of view allows us to recast the above equation as
$$
  \sum_{x\preceq y\preceq z}\mu(y,z)=\delta_{x,z}.
$$
But enough general theory.
If $\preceq$ is actually a total ordering on $P$, which means (given the finiteness condition) that $P$ is isomorphic to an initial segment of $\mathbb N$ (possibily all of $\mathbb N$), then it is easy to see that one has $\mu(x,y)=\epsilon(y-x)$ where $\epsilon:\mathbb Z\to\{-1,0,1\}$ satisfies $\epsilon(i)=(-1)^i$ if $i\in\{0,1\}$ and $\epsilon(i)=0$ otherwise. This almost trivial computation will allow deducing many important Möbius functions, including the classical Möbius function, using the following
Theorem. If $(P_i,\preceq_i)$ are appropriate posets for $i\in I$, and $(P,\preceq)$ is the product poset (restricted if $I$ is infinite), then the Möbius function for $P$ is the product of the Möbius functions for the posets $P_i$: one has $\mu((x_i)_{i\in I},(y_i)_{i\in I})=\prod_{i\in I}\mu_i(x_i,y_i)$.
Here a product poset is the Cartesian product $\prod_{i\in I}P_i$ with a component-wise partial ordering: $(x_i)_{i\in I}\leq(y_i)_{i\in I}$ if and only if $x_i\leq y_i$ in $P_i$ for all $i\in I$. Also "restricted" means restriction to the subset of the Cartesian product of those $(x_i)_{i\in I}$ such that $x_i$ is a minimal element of $P_i$ for all but finitiely many indices $i$; this restriction is necessary to ensure that there are only finitely many elements below any given element of $P$, so that its Möbius function is well defined. The example that will interest us is with $I$ the set $\mathrm{Pr}$ of prime numbers, and every $(P_i,\preceq_i)$ a copy of $(\mathbb N,\leq)$; then elements of $P$ are collections $(m_p)_{p\in\mathrm{Pr}}$ with $m_p\in\mathbb N$ and $m_p=0$ for all but finitiely many $p$. Such an element can be associated to the product $\prod_{p\in\mathrm{Pr}}p^{m_p}$, which is a bijection of $P$ to the set $\mathbb N_{\geq1}$, under which bijection the partial order $\preceq$ on $P$ corresponds to the divisibility relation on $\mathbb N_{\geq1}$. Thus we see that the theorem describes in particular the classical Möbius function. Indeed it gives $\mu(x,y)=\prod_{p\in\mathrm{Pr}}\epsilon(m_p(y)-m_p(x))$ where $m_p(x)$ denotes the multiplicity of the prime $p$ in the factorisation of $x$; the product equals $\prod_{p\in\mathrm{Pr}}\epsilon(m_p(y/x))$ whenever $x\mid y$, which is the usual definition of $\mu(y/x)$.
Proof. This is simply a formal verification that the proposed Möbius function for $P$ satifies the required equation. Let $x=(x_i)_{i\in I}$ and $z=(z_i)_{i\in I}$ be elements of $P$, and let $J=\{i_1,\ldots,i_n\}$ be a finite set containing all indices for which $z_i$ is not minimal. Then
$$
\begin{align}
  \sum_{y\preceq z}\mu(x,y)
   &=\sum_{y_{i_1}\preceq z_{i_1}}\cdots\sum_{y_{i_n}\preceq z_{i_n}}
    \mu_{i_1}(x_{i_1},y_{i_1})\cdots\mu_{i_n}(x_{i_n},y_{i_n}) \\
&
   =\prod_{j\in J}\,\sum_{y_j\preceq z_j}\mu_j(x_j,y_j)
   =\prod_{j\in J}\delta_{x_j,z_j}
   =\delta_{x,z}.\quad\mbox{QED}\\
\end{align}
$$
A: Setting $n=1$ we deduce $\mu(1)=1$ we could continue with 
n=p prime to get $\mu(1)+\mu(p)=0 \implies \mu(p)=-1$, try the product of two primes to get $\mu(p\cdot q)=1$ and so on...
We could too define the arithmetical function u : u(k)=1 for all k and rewrite your definition as $\mu * u = I$ (of course this is just a shortcut for $\sum_{d|n}\mu(d) u(\frac{n}{d})=I(n)$ with $I$ the function of n defined at the right of your equation).
This and $\mu(1)=1$ is enough to prove the Möbius inversion formula (details in Apostol's excellent introduction to A.N.T. page 30 to 32 : the proofs there only require your definition). You could prove that u and I are multiplicative ($f(m n)=f(m)f(n)$ if $\gcd(m,n)=1$) so that $\mu$ will also be multiplicative (using/proving Theorem 2.16 page 36). Of course this is only the draft of an action plan and I'll let you or others concretize it! :-) 
A: Here's another solution using Dirichlet series. Consider the ring of sequences $a_n, n \in \mathbb{N}$ under Dirichlet convolution
$$a_n * b_n = \sum_{d | n} a_d b_{n/d}.$$
This ring has identity the sequence $a_1 = 1, a_n = 0, n \ge 2$. We are asked to find the inverse of the sequence $a_n = 1$ in this ring. A basic observation is that the map
$$a_n \mapsto \sum_{n \in \mathbb{N}} \frac{a_n}{n^s}$$
is a ring homomorphism (where the product on power series is the standard product). The sequence $a_n = 1$ is the Riemann zeta function
$$\zeta(s) = \sum_{n \in \mathbb{N}} \frac{1}{n^s}.$$
Another basic observation, the translation to Dirichlet series of the product poset observation in Marc van Leeuwen's answer, is that the zeta function has an Euler product
$$\zeta(s) = \prod_{p \text{ prime}} \left( \frac{1}{1 - p^{-s}} \right)$$
and hence its inverse has Euler product
$$\zeta(s)^{-1} = \sum_{n \in \mathbb{N}} \frac{\mu(n)}{n^s} = \prod_{p \text{ prime}} \left( 1 - \frac{1}{p^s} \right).$$
The conclusion follows by expanding the RHS. 
A: Yes, there is! It is proved in the same way you prove the Möbius inversion formula:
Going directly from the formula you give, it is clear that
$$ \mu(1) = 1, \mu(p) = -1 \text{ for $p$ a prime}. $$
In general, for square-free numbers, we can show the desired formula by principle of inclusion-exclusion. For instance, for $n = pq, p, q$ primes, we have
$$ \mu(pq) = \sum_{d|pq} \mu(d) - \sum_{d|p} \mu(d) - \sum_{d|q} \mu(d) + \sum_{d|1} \mu(d) = 1 $$
and for three primes,
$$ \mu(pqr) = \left (\sum_{d|pqr} - \sum_{d|pq} - \sum_{d|qr} - \sum_{d|pr} + \sum_{d|p} + \sum_{d|q} + \sum_{d|r} - \sum_{d|1} \right )\mu(d) = -1. $$
Now PIE works for the general case as well. If $n = p_1^{e_1} \cdots p_s^{e_s}$, let $m = p_1^{e_1 - 1 } \cdots p_s^{e_s -1}$, i.e. $n$ divided by its square-free part. Then we get that
$$\small \mu(n) = \sum_{d|n} \mu(d) - \sum_{p_1|n} \left( \sum_{d|n/p_1} \mu(d) \right ) + \sum_{p_1, p_2 |n, p_1 \neq p_2} \left ( \sum_{d|n/(p_1 p_2)} \mu(d) \right ) -\cdots + (-1)^s \sum_{d|m} \mu(d) = 0 $$
if $m \neq 1$, as desired.
A: Indeed one can derive it.
So, from the first case of the given formula, we see that
$\mu(1)=1$.
This proves the first condition of the Möbius function.
Let p be a prime.
$\sum_{d|p}\mu(d)=0$ from 2nd case of given formula.
Thus, $\mu(1)+\mu(p)=0$.
i.e. $ 1+ \mu(p)=0$ hence $\mu(p)=-1$ for any prime.
Further, $\sum_{d|p^2}\mu(d)=0$ Thus, $\mu(1)+\mu(p)+\mu(p^2)=0$ hence $1-1+\mu(p^2)=0$ thus $\mu(p^2)=0$.
For $k>2$, a simple induction with $k=3$ as base case will suffice.
Consider $n={p^k}{q^l}$ with $k,l>2$.
$\sum_{d|{p^k}{q^l}}\mu(d)=\sum_{d|{p^k}}\mu(d)+\sum_{d|{q^l}}\mu(d)+\mu({p^k}{q^l})=0,$
Hence $\mu({p^k}{q^l})=0$. A simple induction will suffice to prove the thing for product of $k$ prime powers. Thus, the Möbius function is zero for any number divisible by a square. This proves the second condition of the Möbius function.
Now, consider two primes $p$ and $q$.
$\sum_{d|pq}\mu(d)=0$, thus $\mu(1)+\mu(p)+\mu(q)+\mu(pq)=0$ i.e. $1-1-1+\mu(pq)=0$; thus,
$\mu(pq)=1$.
Thus, for a product of two primes $p$ and $q$, we have $\mu(pq)=(-1)^2$.
Furthermore, for the product of $k$ primes, we have 
$\sum_{d|{\prod_{i=1}^k p_i}}\mu(d)=\sum_{i=1}^k \binom{k}{i}\mu(p_1p_2\cdots p_i)=0=(1-1)^k=\sum_{i=1}^k \binom{k}{i}{(-1)^i}$. 
Hence, we get $\mu(p_1p_2\cdots p_k)=(-1)^k$. This proves the third condition for the Möbius function, and we are done.
