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

peter dalgaard pdalgd at gmail.com
Fri Feb 12 17:09:34 CET 2016 

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

More information about the R-help mailing list