Method of Lagrange Multipliers with random variables I'll illustrate the issue I'm having with a simple problem.
Let $c_1\in \mathbb{R}$, and $C_2, Z$ be real-valued random variable. Let $u:\mathbb{R} \rightarrow \mathbb{R} $ be a differentiable function, and $f(c_1,Z)$ be a real-valued function that is differentiable with respect to $c_1$.
The problem is:
Maximise: $u(c_1) + \mathbb{E}[u(C_2)]$ 
such that $C_2=f(c_1, Z)$
This problem is easily solved by direct substitution, and the answer is
$u'(c_1) + \mathbb{E}\left[u'(f(c_1,Z)) \frac{\partial f(c_1,Z)}{\partial c_1}\right] = 0 \label{answer}$
The issue is how to write down a Lagrangian whose extrema correspond to the solution of this problem in the normal way.
My first instinct was to write
$\mathcal{L} = u(c_1) + \mathbb{E}[u(C_2)] + \lambda(C_2 - f(c_1,Z))$. 
Now I feel I have to "take the derivative" with respect to the random variable $C_2$, and I'm not sure how to proceed.
Question: How do I write down and solve a Lagrangian whose stationary points correspond to the solution of the above optimisation problem? 
 A: One way to solve is as follows:
Turn the problem into a normal calculus problem, by finding a pdf $g(y)$ folowed by $C_2$ and a real number $c_1$ that is a stationary point of the following Lagrangian
$$ \mathcal{L}(c_1, g) = u(c_1) + \int u(x)g(x) \mathrm{d}x + \int \lambda(y)  \left[ g(y) - \eta(c_1, y)  \right] \mathrm{d}y. $$
Here, $\eta(c_1, y)$ is the pdf followed by the random variable $f(c_1, Z)$. $\lambda(y)$ is a Lagrange multiplier, used to enforce the constraint $g(y)=\eta(c_1,y)$. All integrals are over $\mathbb{R}$. The first order conditions are
$$ u'(c_1) - \frac{\partial}{\partial c_1} \int \lambda(y) \eta(c_1,y) \mathrm{d} y =0 \\
u(y) + \lambda(y) =0 \\
g(y) - \eta(c_1, y) = 0  \\
   $$
Combining the first two gives
$$ u'(c_1) + \frac{\partial}{\partial c_1} \int u(y) \eta(c_1,y) \mathrm{d} y =0 .$$
Now using the 'law of the unconscious statistician' gives
$$ u'(c_1) + \frac{\partial}{\partial c_1} \int u(f(c_1,z)) \phi(z) \mathrm{d} z =0 \\
\Rightarrow u'(c_1) + \int u'(f(c_1,z)) \frac{\partial f(c_1,z)}{\partial c_1} \phi(z) \mathrm{d} z =0 $$
Here, $\phi(z)$ is the pdf followed by the random variable $Z$. This is equivalent to the solution given by OP.
