[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