Bert Gunter
2016-Feb-12 15:53 UTC
[R] Why two curves and numerical integration look so different?
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.
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]]
peter dalgaard
2016-Feb-12 16:09 UTC
[R] Why two curves and numerical integration look so different?
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