[R] Why two curves and numerical integration look so different?
C W
tmrsg11 at gmail.com
Fri Feb 12 16:36:47 CET 2016
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]]
More information about the R-help
mailing list