David Winsemius
2016-Feb-11 19:06 UTC
[R] Why two curves and numerical integration look so different?
> 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
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. 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 > >[[alternative HTML version deleted]]
David Winsemius
2016-Feb-11 20:32 UTC
[R] Why two curves and numerical integration look so different?
> 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