Solve and find x for normally distributed functions

Question

I'm trying to solve a equation and find the value for x from it. I'm using the following code:

mu1 = 0
mu2 = 1
sigma1 = 0.5
sigma2 = 0.6
prior1 = 0.3
prior2 = 0.7

boundary = function(x) {
  return(( 1 / sqrt(2 * pi * sigma1)) * exp(-0.5 * ((x - mu1) / sigma1)^2)*prior1) - 
    ((1 / sqrt(2 * pi * sigma2)) * exp(-0.5 * ((x - mu2) / sigma2)^2)*prior2)
}
uniroot(boundary, c(-1e+05, 1e+07))

This does not give me correct answers. I'm pretty new to R and am not sure exactly how uniroot works.

• Is there a better way to solve this equation?
• Are there any packages available to solve this (similar to MATLAB's solve() function)?

Answer 1

You can shorten your code a little bit (and minimize the chance of typos) by using the built-in dnorm() function:

curve(dnorm(x,mean=0,sd=0.5),from=-4,to=4)
curve(dnorm(x,mean=0,sd=0.6),add=TRUE,col=2)

If I try uniroot() over a reasonable range, I get sensible answers:

uniroot(function(x) dnorm(x,0,0.5)-dnorm(x,0,0.6), c(0,5))  ## 0.546
uniroot(function(x) dnorm(x,0,0.5)-dnorm(x,0,0.6), c(-5,0)) ## -0.546

If I try to start from huge values, the calculation is going to underflow (whether you do it by hand or use dnorm(); the Gaussian goes below R's minimum representable valuable (see ?.Machine) between 19 and 20.

dnorm(19,0,0.5)  ## 2e-314
dnorm(20,0,0.5)  ## 0

In some cases you can use log=TRUE in dnorm() to compute the log-probability and use that in computations instead, but here (as in many Bayesian computations) you're a bit stuck because you have to subtract the results.