[R] Why two curves and numerical integration look so different?

Fri Feb 12 16:57:32 CET 2016

Hi Bert,

Yay fantasyland!

In all seriousness, You are referring to this?
https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f

In particular, you mean this: .Machine$double.eps ^ 0.5?

Thanks!

On Fri, Feb 12, 2016 at 10:53 AM, Bert Gunter <bgunter.4567 at gmail.com>
wrote:

> You are working in fantasyland. Your density is nonsense.
>
> Please see FAQ 7.31 for links to computer precision of numeric
> calculations.
>
>
> Cheers,
> Bert
> Bert Gunter
>
> "The trouble with having an open mind is that people keep coming along
> and sticking things into it."
> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>
>
> On Fri, Feb 12, 2016 at 7:36 AM, C W <tmrsg11 at gmail.com> wrote:
> > Hi David,
> >
> > This is the Gaussian looking distribution I am trying to integrate.
> >
> https://drive.google.com/file/d/0B2xN0-A6iTB4NThIZ2tYdGxHc00/view?usp=sharing
> >
> > Notice the unnormalized density goes up as high as 2.5*101^191.
> >
> > I tried to create small intervals like
> >> seq(0.5, 1.3, by = 10^(-8))
> >
> > but that doesn't seem to be good enough, as we know, it should integrate
> to
> > 1.
> >
> > On Thu, Feb 11, 2016 at 3:32 PM, David Winsemius <dwinsemius at comcast.net
> >
> > wrote:
> >
> >>
> >> > On Feb 11, 2016, at 11:30 AM, C W <tmrsg11 at gmail.com> wrote:
> >> >
> >> > Hi David,
> >> >
> >> > My real function is actually a multivariate normal, the simple toy 1-d
> >> normal won't work.
> >> >
> >> > But, you gave me an idea about restricting the bounds, and focus
> >> integrating on that.  I will get back to you if I need any further
> >> assistance.
> >>
> >> You'll probably need to restrict your bounds even more severely than I
> did
> >> in the 1-d case (using 10 SD's on either side of the mean) . You might
> need
> >> adequate representation of points near the center of your
> hyper-rectangles.
> >> At least that's my armchair notion since I expect the densities tail off
> >> rapidly in the corners. You can shoehorn multivariate integration around
> >> the `integrate` function but it's messy and inefficient. There are other
> >> packages that would be better choices. There's an entire section on
> >> numerical differentiation and integration in CRAN Task View: Numerical
> >> Mathematics.
> >>
> >> --
> >> David.
> >>
> >>
> >> >
> >> > Thank you so much!
> >> >
> >> > On Thu, Feb 11, 2016 at 2:06 PM, David Winsemius <
> dwinsemius at comcast.net>
> >> wrote:
> >> >
> >> > > On Feb 11, 2016, at 9:20 AM, C W <tmrsg11 at gmail.com> wrote:
> >> > >
> >> > > I want to do numerical integration w.r.t. mu: P(mu) × N(mu, 0.00001)
> >> > >
> >> > > Because the variance is small, it results in density like:
> 7.978846e+94
> >> > >
> >> > > Is there any good suggestion for this?
> >> >
> >> > So what's the difficulty? It's rather like the Dirac function.
> >> >
> >> > > integrate( function(x) dnorm(x, sd=0.00001), -.0001,0.0001)
> >> > 1 with absolute error < 7.4e-05
> >> >
> >> >
> >> > --
> >> > David.
> >> >
> >> > >
> >> > > Thanks so much!
> >> > >
> >> > >
> >> > > On Thu, Feb 11, 2016 at 9:14 AM, C W <tmrsg11 at gmail.com> wrote:
> >> > >
> >> > >> Wow, thank you, that was very clear.  Let me give it some more runs
> >> and
> >> > >> investigate this.
> >> > >>
> >> > >> On Thu, Feb 11, 2016 at 12:31 AM, William Dunlap <
> wdunlap at tibco.com>
> >> > >> wrote:
> >> > >>
> >> > >>> Most of the mass of that distribution is within 3e-100 of 2.
> >> > >>> You have to be pretty lucky to have a point in sequence
> >> > >>> land there.  (You will get at most one point there because
> >> > >>> the difference between 2 and its nearest neightbors is on
> >> > >>> the order of 1e-16.)
> >> > >>>
> >> > >>> seq(-2,4,len=101), as used by default in curve, does include 2
> >> > >>> but seq(-3,4,len=101) and seq(-2,4,len=100) do not so
> >> > >>> curve(..., -3, 4, 101) and curve(..., -2, 4, 100) will not show
> the
> >> bump.
> >> > >>> The same principal holds for numerical integration.
> >> > >>>
> >> > >>>
> >> > >>> Bill Dunlap
> >> > >>> TIBCO Software
> >> > >>> wdunlap tibco.com
> >> > >>>
> >> > >>> On Wed, Feb 10, 2016 at 6:37 PM, C W <tmrsg11 at gmail.com> wrote:
> >> > >>>
> >> > >>>> Dear R,
> >> > >>>>
> >> > >>>> I am graphing the following normal density curve.  Why does it
> look
> >> so
> >> > >>>> different?
> >> > >>>>
> >> > >>>> # the curves
> >> > >>>> x <- seq(-2, 4, by=0.00001)
> >> > >>>> curve(dnorm(x, 2, 10^(-100)), -4, 4)  #right answer
> >> > >>>> curve(dnorm(x, 2, 10^(-100)), -3, 4)  #changed -4 to -3, I get
> wrong
> >> > >>>> answer
> >> > >>>>
> >> > >>>> Why the second curve is flat?  I just changed it from -4 to -3.
> >> There is
> >> > >>>> no density in that region.
> >> > >>>>
> >> > >>>>
> >> > >>>> Also, I am doing numerical integration.  Why are they so
> different?
> >> > >>>>
> >> > >>>>> x <- seq(-2, 4, by=0.00001)
> >> > >>>>> sum(x*dnorm(x, 2, 10^(-100)))*0.00001
> >> > >>>> [1] 7.978846e+94
> >> > >>>>> x <- seq(-1, 4, by=0.00001) #changed -2 to -1
> >> > >>>>> sum(x*dnorm(x, 2, 10^(-100)))*0.00001
> >> > >>>> [1] 0
> >> > >>>>
> >> > >>>> What is going here?  What a I doing wrong?
> >> > >>>>
> >> > >>>> Thanks so much!
> >> > >>>>
> >> > >>>>        [[alternative HTML version deleted]]
> >> > >>>>
> >> > >>>> ______________________________________________
> >> > >>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more,
> see
> >> > >>>> https://stat.ethz.ch/mailman/listinfo/r-help
> >> > >>>> PLEASE do read the posting guide
> >> > >>>> http://www.R-project.org/posting-guide.html
> >> > >>>> and provide commented, minimal, self-contained, reproducible
> code.
> >> > >>>>
> >> > >>>
> >> > >>>
> >> > >>
> >> > >
> >> > >       [[alternative HTML version deleted]]
> >> > >
> >> > > ______________________________________________
> >> > > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >> > > https://stat.ethz.ch/mailman/listinfo/r-help
> >> > > PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> > > and provide commented, minimal, self-contained, reproducible code.
> >> >
> >> > David Winsemius
> >> > Alameda, CA, USA
> >> >
> >> >
> >>
> >> David Winsemius
> >> Alameda, CA, USA
> >>
> >>
> >
> >         [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>

	[[alternative HTML version deleted]]