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

[Bert Gunter]

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.