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

[David Winsemius]

> 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