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]



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]

enter image description here


2 Answers 2


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]

Improved figure from OP

  • $\begingroup$ thanks for your answer. I'm trying to convert this Mathematica code into another language (Python). Any way to output something from Mathematica that is easily usable in other languages? For example, the Rationalize function doesn't seem to be available in Python. Sorry, I know this is may be off topic for a Mathematica stackexchange question but I frequently find myself trying to use Mathematica results elsewhere and have trouble translating... $\endgroup$
    – ben18785
    Feb 19, 2019 at 18:11
  • 1
    $\begingroup$ You'll certainly need to use a multiple precision arithmetic package if you translate it to Python or R. I don't know Python but a google search finds "mpmath" used with "gmpy". $\endgroup$
    – JimB
    Feb 19, 2019 at 18:28

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.

  • 2
    $\begingroup$ As you've now noticed, the key is to use the three-argument form of Gamma[]. Here's a further simplification of your formula: fderiv[x_, μ_, n_, a_, b_] := (n - 1)/(x - μ) + ((Abs[x - μ]/a)^(n - 1) Exp[-((x - μ)^2/(2 a^2))] - (Abs[x - μ]/b)^(n - 1) Exp[-((x - μ)^2/(2 b^2))])/(2^((n - 3)/2) (x - μ) Gamma[(n - 1)/2, (x - μ)^2/(2 a^2), (x - μ)^2/(2 b^2)]) $\endgroup$ Mar 12, 2019 at 2:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.