0

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)?
Prav
  • 15
  • 1
  • 4
  • So this function plots as a normal curve, but it will never actually reach 0 will it? It will just get smaller and smaller fading from x=0 to more extreme negative or positive values. I'm not sure I understand what you are trying to "solve". – thelatemail Feb 06 '14 at 02:59
  • yes,those are the equations of two normal curves. But there is only one unknown x in the equation. which means i should be able to solve it and find its value. – Prav Feb 06 '14 at 03:03
  • 1
    It isn't clear what you are trying to solve. What is x doing? I suspect there will be an easier way to do whatever it is you are trying to do, once we know what that is. Also, note that your function is named `boundary`, whereas you are calling `decisionBoundary` w/i `unitroot`. Lastly, the `prior`s don't seem to be playing any role. – gung - Reinstate Monica Feb 06 '14 at 03:10
  • I'm trying to do bayesian analysis on normal distributions. The equation represents my discriminant function. solving for x will give me the decision boundary which separate's the two normal distributions. I have updated my question with the priors and corrected the uniroot function as well. – Prav Feb 06 '14 at 03:14
  • Re Q1:....analytically? – mathematical.coffee Feb 06 '14 at 03:15

1 Answers1

2

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.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453