Derivative of integrated noise Gaussian likelihood

Question

In a Bayesian problem with Gaussian likelihood with mean $\mu$ and a uniform prior on the standard deviation $\sigma$, it is possible to derive the marginal posterior (where $\sigma$ has been integrated out of the joint).

$p(\mu) = \int_a^b \mathcal{N}(x|\mu, \sigma) U(\sigma|a,b) p(\mu) \mathbb{d}\sigma$

which can be done in Mathematica (omitting $p(\mu)$ as it doesn't affect the integration):

Assuming[{a > 0, b > 0, n > 0, sse > 0, b > a}, 
 Integrate[(1/(Sqrt[2 Pi sigma^2]))^n
 Exp[-(x - mu)^2/(2 sigma^2)] PDF[UniformDistribution[{a, b}], sigma], {sigma, a, b}]]

To yield this expression,

(\[Pi]^(-n/2) (1/(mu - x)^2)^(-(1/2) + n/2) (Gamma[1/2 (-1 + n), (mu - x)^2/(2 a^2)] - 
Gamma[1/2 (-1 + n), (mu - x)^2/(2 b^2)]))/(2 Sqrt[2] (a - b))

Letting $g$ equal the log of the above expression, I can determine its derivative wrt x,

FullSimplify@D[g, mu]

which equals this horror,

 fDeriv2[x_, mu_, n_, a_, b_] := (-2^(3/2 - n/2) E^(-((mu - x)^2/(2 a^2))) ((mu - x)^2/a^2)^(1/2 (-1 + n)) +
   2^(3/2 - n/2) E^(-((mu - x)^2/(2 b^2))) ((mu - x)^2/b^2)^(
   1/2 (-1 + n)) - (-1 + n) (Gamma[1/2 (-1 + n), (mu - x)^2/(2 a^2)] - 
   Gamma[1/2 (-1 + n), (mu - x)^2/(2 b^2)]))/((mu - x) (Gamma[1/2 (-1 + n), (mu - x)^2/(2 a^2)] - 
   Gamma[1/2 (-1 + n), (mu - x)^2/(2 b^2)]))

The issue with this expression is that it becomes unstable when $x - mu$ is small and $n$ is large.

For example,

fDeriv2[0.1, 0.2, 100, 2, 4]

yields

ComplexInfinity

I know that the derivative exists but the denominator and the numerator are either really big or small which, with numerical under/over-flow leads the calculation to fail.

Does anyone know how I can derive a more 'friendly' expression for the derivative that doesn't have these pathologies?

Edit: to further illustrate the pathologies of this expression, I plot it as a function of x:

Plot[fDeriv2[x, 10, 100, 2, 4], {x, 10, 100}, PlotRange -> Full]

Answer 1

You can write your function as a Piecewise function or use an If statement that attempts to take care of at least some of the numerical issues.

When x = mu, then you have a few places where you end up dividing by zero. It might make sense to have fDeriv2 use the limit (if it exists) as x -> mu.

Table[Limit[fDeriv2[x, mu, 5, a, b], x -> mu], {n, 1, 10}]
(* {0, 0, 0, 0, 0, 0, 0, 0, 0, 0} *)

So it seems that fDeriv2 should be defined as zero if x = mu. But the bigger issue might be where n is large. To deal with that you might try something like the following:

fDeriv[x_, mu_, n_, a_, b_] := If[x == mu, 0, 
  Block[{$MaxExtraPrecision = 1000}, N[fDeriv2[x, mu, n, a, b], 30]]]

fDeriv[1/10, 2/10, 100, 2, 4]
(* -0.0245049387169514315208885210720 *)

You can make the plot in your example work by using

Plot[fDeriv[Rationalize[x, 0], 10, 100, 2, 4], {x, 1, 100}, PlotRange -> Full]

Answer 2

The answer to this is actually much simpler than it appears since I failed to notice that there is a common term (differences of two incomplete Gammas) in the above. This means that we can avoid the numerical issues above by the following expression (note, not an approximation):

fDeriv[x_, mu_, n_, a_, b_] := ((-2^(3/2 - n/2) E^(-((mu - x)^2/(2 a^2))) ((mu - x)^2/
          a^2)^(1/2 (-1 + n)) + 
      2^(3/2 - n/2) E^(-((mu - x)^2/(2 b^2))) ((mu - x)^2/
          b^2)^(1/2 (-1 + n)))/((mu - x) Gamma[
       1/2 (-1 + n), (mu - x)^2/(2 a^2), (mu - x)^2/(2 b^2)]) - (-1 + 
      n)/(mu - x))

which when plotted produces a graph indistinguishable from that of JimB's answer.