Solve for and Plot the Relationship Between Mean and Standard Deviation of a Normal Distribution Conditional on Satisfaction of A System of Equations

Question

I am trying to use Mathematica, R, or Matlab to solve for (since it cannot seem to be solved analytically) and plot the relationship between mean and standard deviation of a normal distribution conditional on the satisfaction of below system of equations (the motivation is that I am trying to investigate the relationship between mean and standard deviation conditional on the normal CDF is tangent with the y=x line in the [0, 1] domain):

$$ 1=\frac{1}{\sigma \sqrt{2\pi}}e^{-\frac{-(x-\mu)^2}{2\sigma^2}}\\ x=\frac{1}{2}\left[ 1 + erf(\frac{x-\mu}{\sqrt{2}\sigma}) \right] $$

The first equation specifies that the normal CDF is tangent to the 45 degree line y=x (which is equivalent to setting the normal PDF, which is the derivative of normal CDF, equal to 1). The second equation specifies that the normal CDF intersect with the 45 degree line y=x. The objective is to ascertain the relationship between $\sigma$ and $\mu$, when the condition as specified by the system of equations is satisfied.

So this is a system of two equation with three unknowns: $x$, $\sigma$, and $\mu$, and the objective is to eliminate x and reveal the relationship between $x$ and $\sigma$. The problem is with the error function erf(), which is an integral. Thus, when you solve for x as a function of $\mu$ and $\sigma$, and then substitute x into the second equation in the hope of revealing the relationship between $\mu$ and $\sigma$ under the specified condition, the equation cannot be solved analytically.

I wonder if anyone could offer help as to how to solved this using mathematica, R, or Matlab and to plot the relationship between $\mu$ and $\sigma$ with $\mu$ on the horizontal axis and within the [0, 1] domain only.

Graphically, it should look like a "little hill", with $\mu$ on the horizontal axis within domain [0, 1], and $\sigma$ on the vertical axis within domain (0, $\frac{1}{\sqrt{2\pi}}$]. The top of the little hill is at $\mu = 0.5, \sigma=\frac{1}{\sqrt{2\pi}}$.

Thank you for your help. I would appreciate any input in solving this problem.

Answer 1

Let $z=\frac{x-\mu}{\sigma}$, so your first equation gives $$z = \pm\sqrt{-\log(2\pi\sigma^2)} \tag{1} $$ and your second equation gives $$\mu=\Phi(z)-\sigma\,z\tag{2} $$ where $\Phi$ is the standard normal CDF. Substituting (1) into (2), we obtain (two branches of) $\mu$ as a function of $\sigma$:

$$\mu_a (\sigma) = \Phi\left(\sqrt{-\log(2\pi\sigma^2)}\right)-\sigma\,\sqrt{-\log(2\pi\sigma^2)} $$ $$\mu_b (\sigma) = \Phi\left(-\sqrt{-\log(2\pi\sigma^2)}\right)+\sigma\,\sqrt{-\log(2\pi\sigma^2)} = 1-\mu_a(\sigma) $$ where the latter equality follows from the fact that $\Phi(-t) = 1- \Phi(t)$ for all $t$.

This allows a plot of $\sigma$ vs. $\mu$ to be generated by forming appropriate lists of coordinate pairs $(\mu(\sigma),\,\sigma)$ with $\sigma$ varying in the interval $(0,\frac{1}{\sqrt{2\,\pi}}]$:

Sage code:

def Phi(z): return (1 + erf(z/sqrt(2)))/2
def mu_a(sigma): z = sqrt(-log(2*pi*sigma^2)); return (Phi(z) - sigma*z)
def mu_b(sigma): return 1 - mu_a(sigma)

maxsig = 1/sqrt(2*pi)
npoints = 500
La = [( mu_a((k/npoints)*maxsig), (k/npoints)*maxsig ) for k in [1..npoints]]
Lb = [( mu_b((k/npoints)*maxsig), (k/npoints)*maxsig ) for k in [1..npoints]]
g = Graphics()
g += line(La, color="blue")
g += line(Lb, color="red")
g.show(gridlines=["major","major"],  title='', axes_labels=['mu','sigma'], figsize=3)

---

R code producing the same result:

mu_a <- function(sigma) { 
   z <- sqrt(-log(2*pi*sigma^2))
   m <- pnorm(z) - sigma*z
   return(m)}
mu_b <- function(sigma) {1-mu_a(sigma)} 

maxsig <- 1/sqrt(2*pi)
vsig <- (c(1:500)/500)*maxsig
vmu_a <- lapply(vsig,mu_a)
vmu_b <- lapply(vsig,mu_b)

plot(vmu_a,vsig,type="l",xlim=c(0, 1), ylim=c(0, 0.4), 
     xlab="mu", ylab="sigma", col="blue")
lines(vmu_b,vsig,type="l", col="red")