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

peter dalgaard pdalgd at gmail.com
Fri Feb 12 21:54:35 CET 2016 

> On 12 Feb 2016, at 17:44 , C W <tmrsg11 at gmail.com> wrote:
> 
> On a side note, is it ok to do?
> 
> > which(max(p_x)) 
> and use that instead of numerical integration to get E[X]?

Now THAT makes absolutely no sense. max() is a number, so which(max()) usually returns 1.

If you mean whether the mode is equal to the mean: Only if the distribution is symmetric and unimodal.

-pd

> 
> I will try both and report back!  Thank you expeRts
> 
> On Fri, Feb 12, 2016 at 11:29 AM, C W <tmrsg11 at gmail.com> wrote:
> Hi Peter,
> 
> Great, let me try that and get back to you on my findings in a few hours!  :)
> 
> On Fri, Feb 12, 2016 at 11:09 AM, peter dalgaard <pdalgd at gmail.com> wrote:
> I don't see here FAQ 7.31 comes in either (for once!...)
> 
> However, either the density is unnormalized and the integral is not 1, or the integral is 1 and it is normalized. The one in the picture clearly does not integrate to one. You can fit a rectangle of size 0.1 by 1e191 under the curve so the integral should be > 1e190 .
> 
> As depicted, I don't see why a plain integral from .5 to 1.5 shouldn't work.
> 
> -pd
> 
> On 12 Feb 2016, at 16:57 , C W <tmrsg11 at gmail.com> wrote:
> 
> > 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]]
> >
> > ______________________________________________
> > 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.
> 
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Office: A 4.23
> Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com
> 
> 
> 

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com

More information about the R-help mailing list