項目反応理論やその他の統計モデルで, 「正規分布で期待値をとる」という操作は 時々あるのではないかと思います。IRTの教科書や 井手さんの論文 を読んでいて, これはガウス-エルミート求積法 (Gauss-Hermite quadrature)という方法で簡単かつ高精度に計算できることを知りました。 以下は知っている人には当たり前ですが, 自分用のメモ代わり。

ガウス-エルミート求積法は, 関数空間で直交するエルミート多項式を使って, N個の分点(abscissas)によって関数を(2N-1)次の多項式で近似することで, 積分を効率よく求める 方法のようです。数学的な詳細は置いておくと, 関数 f(x) を exp(-x^2) で期待値を とった値は以下のように書くことができる, という話。

∫^∞_-∞e^{-x^2}f(x)dx ~= Σ_k=1^N w_kf(x_k).

ここで必要な分点 x_kと重み w_kは, MATLABでは http://hips.seas.harvard.edu/content/gauss-hermite-quadrature-matlab の hermquad.m で求めることができます。 Pythonのnumpyでは, numpy.polynomial.hermite.hermgauss で計算することができるようです。
話としてはわかるのですが, 実際にどのくらい精度が出るのか実感がなかったので, 実際に試してみました。

上の式は e^{-x^2}で期待値をとった形になっているため, 標準正規分布の形である e^{-x^2/2}で期待値をとるには, x^2/2 = t^2 すなわち, x=√2tと変数変換する必要があります。 このとき dx/dt = √2から dx=√2dt なので,

I = ∫^∞_-∞ 1/√(2π) e^{-x^2/2} f(x)dx
= ∫^∞_-∞ 1/√(2π) e^{-t^2} f(√2t)√2dt
= ∫^∞_-∞ 1/√π e^{-t^2} f(√2t)dt
= 1/√π Σ_k=1^N w_k f(√2x_k) です。
いま, f(x) = N(1,1) = 1/√(2π) exp[-(x-1)^2/2] としてみます。このとき, I = ∫^∞_-∞ N(x|0,1) f(x)dxは解析的に計算できて, 少し計算すると

I = e^-1/4 / (2√π) = 0.219695644733861

になります。
これを, 上のガウス-エルミート求積法でどのくらい正確に近似できるかMATLABで 計算してみました。ghtest(N) のNが分点の数を表します。

>> format long
>> I = ghtest(10)
I =
   0.219700135209121
>> I = ghtest(20) 
I =
   0.219695644690565
>> I = ghtest(30)
I =
   0.219695644733860
>> I = ghtest(50)
I =
   0.219695644733861
>> I = ghtest(100)
I =
   0.219695644733860

というわけで, すさまじく高精度に近似できて, 分点の数は20個程度でほぼ充分, 30個以上では値がほとんど同じということがわかりました。 (ただし今回は被積分関数が単純なので, もっと複雑な場合はここまでではない 可能性があります。)
ちなみに, 分点の数が20の時に真の値との差は 0.000000000043296 でした。よって, 安心して上のガウス-エルミート求積法で正規分布に関する, ほぼ正しい期待値を計算することができるようです。
なお, 少し調べるとガウス-エルミート求積が積分値の高次のラプラス近似にあたる という話が, 94年のBiometrika にありました。