Re: [R] integrate lgamma from 0 to Inf
IMHO, you should consult an advanced text on calculus with integration over infinite intervals, so that you understand what you are trying to do. Joseph F. Lucke Senior Statistician Research Institute on Addictions University at Buffalo State University of New York 1021 Main Street Buffalo, NY 14203-1016 Office: 716-887-6807 Fax: 716-887-2510 http://www.ria.buffalo.edu/profiles/lucke.html Stavros Macrakis Sent by: r-help-boun...@r-project.org 04/27/2009 10:30 AM To Andreas Wittmann cc r-help@r-project.org Subject Re: [R] integrate lgamma from 0 to Inf On Wed, Apr 22, 2009 at 3:28 AM, Andreas Wittmann wrote: > i try to integrate lgamma from 0 to Inf. Both gamma and log are positive and monotonically increasing for large arguments. What can you conclude about the integrability of log(gamma(x))? -s __ R-help@r-project.org mailing list 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@r-project.org mailing list 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.
Re: [R] integrate lgamma from 0 to Inf
On Wed, Apr 22, 2009 at 3:28 AM, Andreas Wittmann wrote: > i try to integrate lgamma from 0 to Inf. Both gamma and log are positive and monotonically increasing for large arguments. What can you conclude about the integrability of log(gamma(x))? -s __ R-help@r-project.org mailing list 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.
Re: [R] integrate lgamma from 0 to Inf
I am late to this discussion so forgive me if I am being redundant. It
appears to me that the integral from 0 to infinity of log-gamma may
diverge to infinity. Equation 6.1.50 of Abramowitz & Stegun shows that a
re-expression of lnGamma(x) and the Stirling approximation involves
(x-1/2)*log(x) -x along with other terms. This term appears to dominate
the integral and itself diverge. It is worth checking out.
> integrate(lgamma, lower = 0, upper = 10)
43.24636 with absolute error < 1.6e-11
> integrate(lgamma, lower = 0, upper = 100)
15438.12 with absolute error < 1.4
> integrate(lgamma, lower = 0, upper = 1000)
2701843 with absolute error < 254
> integrate(lgamma, lower = 0, upper = 1)
385485116 with absolute error < 18464
> gterm <- function(x){(x-.5)*log(x)-x}
> integrate(gterm,lower=0,upper=10)
33.61633 with absolute error < 1.1e-05
> integrate(gterm,lower=0,upper=100)
15345.59 with absolute error < 1.2
> integrate(gterm,lower=0,upper=1000)
2700923 with absolute error < 252
> integrate(gterm,lower=0,upper=1)
385475926 with absolute error < 18462
>
Joseph F. Lucke
Senior Statistician
Research Institute on Addictions
University at Buffalo
State University of New York
1021 Main Street
Buffalo, NY 14203-1016
Office: 716-887-6807
Fax: 716-887-2510
http://www.ria.buffalo.edu/profiles/lucke.html
"Andreas Wittmann"
Sent by: r-help-boun...@r-project.org
04/22/2009 03:28 AM
To
r-help@r-project.org
cc
Subject
[R] integrate lgamma from 0 to Inf
Dear R users,
i try to integrate lgamma from 0 to Inf. But here i get the message
"roundoff error is detected in the extrapolation table", if i use 1.0e120
instead of Inf the computation works, but this is against the suggestion
of integrates help information to use Inf explicitly. Using stirlings
approximation doesnt bring the solution too.
## Stirlings approximation
lgammaApprox <- function(x)
{
0.5*log(2*pi)+(x-(1/2))*log(x)-x
}
integrate(lgamma, lower = 0, upper = 1.0e120)
integrate(lgammaApprox, lower = 0, upper = 1.0e120)
> integrate(lgamma, lower = 0, upper = 1.0e120)
1.374051e+242 with absolute error < 3.2e+235
> integrate(lgammaApprox, lower = 0, upper = 1.0e120)
1.374051e+242 with absolute error < 3.2e+235
integrate(lgamma, lower = 0, upper = Inf)
integrate(lgammaApprox, lower = 0, upper = Inf)
> integrate(lgamma, lower = 0, upper = Inf)
Fehler in integrate(lgamma, lower = 0, upper = Inf) :
roundoff error is detected in the extrapolation table
> integrate(lgammaApprox, lower = 0, upper = Inf)
Fehler in integrate(lgammaApprox, lower = 0, upper = Inf) :
roundoff error is detected in the extrapolation table
Many thanks if you have any advice for me!
best regards
Andreas
--
__
R-help@r-project.org mailing list
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@r-project.org mailing list
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.
[R] integrate lgamma from 0 to Inf
Dear R users,
i try to integrate lgamma from 0 to Inf. But here i get the message "roundoff
error is detected in the extrapolation table", if i use 1.0e120 instead of Inf
the computation works, but this is against the suggestion of integrates help
information to use Inf explicitly. Using stirlings approximation doesnt bring
the solution too.
## Stirlings approximation
lgammaApprox <- function(x)
{
0.5*log(2*pi)+(x-(1/2))*log(x)-x
}
integrate(lgamma, lower = 0, upper = 1.0e120)
integrate(lgammaApprox, lower = 0, upper = 1.0e120)
> integrate(lgamma, lower = 0, upper = 1.0e120)
1.374051e+242 with absolute error < 3.2e+235
> integrate(lgammaApprox, lower = 0, upper = 1.0e120)
1.374051e+242 with absolute error < 3.2e+235
integrate(lgamma, lower = 0, upper = Inf)
integrate(lgammaApprox, lower = 0, upper = Inf)
> integrate(lgamma, lower = 0, upper = Inf)
Fehler in integrate(lgamma, lower = 0, upper = Inf) :
roundoff error is detected in the extrapolation table
> integrate(lgammaApprox, lower = 0, upper = Inf)
Fehler in integrate(lgammaApprox, lower = 0, upper = Inf) :
roundoff error is detected in the extrapolation table
Many thanks if you have any advice for me!
best regards
Andreas
--
__
R-help@r-project.org mailing list
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.
Re: [R] Integrate function
.. but it turned out he wanted; integrate()$value -- David Winsemius On Dec 22, 2008, at 6:26 PM, David Winsemius wrote: If these messages you're hearing are warnings, then the answer might be: ?warnings -- David Winsemius On Dec 22, 2008, at 6:07 PM, glenn roberts wrote: Quick One if any one can help please. On use of integration function ‘integrate’; how do I get the function to return just the value with no messages please Glenn [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
Re: [R] Integrate function
If these messages you're hearing are warnings, then the answer might be: ?warnings -- David Winsemius On Dec 22, 2008, at 6:07 PM, glenn roberts wrote: Quick One if any one can help please. On use of integration function ‘integrate’; how do I get the function to return just the value with no messages please Glenn [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. __ R-help@r-project.org mailing list 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.
[R] Integrate function
Quick One if any one can help please. On use of integration function integrate¹; how do I get the function to return just the value with no messages please Glenn [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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.
[R] integrate problem
i'm in the process of switch from Maple to R and am trying to code the following function: w(d)=\int -Inf_Inf A(x) |int 0_(t(d)-x) Fy dydx can some one point me in the right direction? i don't seem to be able to figure it out on my own. Dr. Wade Winterhalter University of Central Florida Department of Biology ph: 407-260-0134 __ R-help@r-project.org mailing list 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.
[R] integrate: can pass numbers but not variables as additional arguments
I'm trying to integrate one dimension of a bivariate normal. It works
when additional parameters are passed explicitly:
library(mvtnorm)
bivInt <- function(x,y,mx,my,r) { dmvnorm(c(x, y), mean=c(mx, my),
sigma=rbind(c(1, r), c(r, 1))) }
integrate(Vectorize(bivInt), lower=-Inf, upper=2, 1, 2, 2, .5)$value
# [1] 0.1737709
How to make it work when the additional parameters are passed as variables?
I'm getting errors:
y<-1; mx<-2; my<-2; r<-.5
integrate(Vectorize(bivInt), lower=-Inf, upper=2, y, mx, my, r)
# Error in eval(expr, envir, enclos) :
# ..1 used in an incorrect context, no ... to look in
> integrate(Vectorize(bivInt), lower=-Inf, upper=2, y=y, mx=mx, my=my, r=r)
Error in eval(expr, envir, enclos) :
..1 used in an incorrect context, no ... to look in
__
R-help@r-project.org mailing list
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.
Re: [R] Integrate a 1-variable function with 1 parameter (Jose L.Romero)
Hi,
The answer can be obtained in "closed" form using the "pgamma" function,
which is closely related to the incomplete gamma function, as follows:
integrand <- function (t, x) {
exp(-2*t)*(2*t)^x/(10*factorial(x))
}
upper <- 10
x <- 0:44
ans1 <- sapply(x, function(x) integrate(integrand, lower=0, upper=upper,
x=x))
ans2 <- gamma(x+1) * pgamma(q=2*upper, shape=x+1, rate = 1, scale = 1,
lower.tail = TRUE) / (20*factorial(x))# using the "pgamma" function
cbind(x=x, ans1=unlist(ans1[1,]), ans2=ans2) # both answers are identical
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Moshe Olshansky
Sent: Wednesday, August 27, 2008 9:58 PM
To: r-help@r-project.org; [EMAIL PROTECTED]
Subject: Re: [R] Integrate a 1-variable function with 1 parameter (Jose
L.Romero)
This can be done analytically: after changing a variable (2*t -> t) and some
scaling we need to compute
f(x) = integral from 0 to 20 of (t^x*exp(-t))dt/factorial(x)
f(0) = int from 0 to 20 of exp(-t)dt = 1 - exp(-20) and integration by parts
yields (for x=1,2,3,...)
f(x) = -exp(-20)*20^x/factorial(x) + f(x-1) so that
f(x) = 1 - exp(-20)*sum(20^k/factorial(k)) where the sum is for k=0,1,...,x
If I did not a mistake, your original quantity should be f(x)/20.
--- On Thu, 28/8/08, jose romero <[EMAIL PROTECTED]> wrote:
> From: jose romero <[EMAIL PROTECTED]>
> Subject: [R] Integrate a 1-variable function with 1 parameter (Jose L.
> Romero)
> To: r-help@r-project.org
> Received: Thursday, 28 August, 2008, 12:23 AM Hey fellas:
>
> I would like to integrate the following function:
>
> integrand <- function (x,t) {
> exp(-2*t)*(2*t)^x/(10*factorial(x))
> }
>
> with respect to the t variable, from 0 to 10.
> The variable x here works as a parameter: I would like to integrate
> the said function for each value of x in 0,1,..,44.
>
> I have tried Vectorize to no avail.
>
> Thanks in advance,
> jose romero
>
> __
> R-help@r-project.org mailing list
> 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.
__
R-help@r-project.org mailing list
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.
__
R-help@r-project.org mailing list
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.
Re: [R] Integrate a 1-variable function with 1 parameter (Jose L. Romero)
This can be done analytically: after changing a variable (2*t -> t) and some
scaling we need to compute
f(x) = integral from 0 to 20 of (t^x*exp(-t))dt/factorial(x)
f(0) = int from 0 to 20 of exp(-t)dt = 1 - exp(-20)
and integration by parts yields (for x=1,2,3,...)
f(x) = -exp(-20)*20^x/factorial(x) + f(x-1) so that
f(x) = 1 - exp(-20)*sum(20^k/factorial(k)) where the sum is for k=0,1,...,x
If I did not a mistake, your original quantity should be f(x)/20.
--- On Thu, 28/8/08, jose romero <[EMAIL PROTECTED]> wrote:
> From: jose romero <[EMAIL PROTECTED]>
> Subject: [R] Integrate a 1-variable function with 1 parameter (Jose L. Romero)
> To: r-help@r-project.org
> Received: Thursday, 28 August, 2008, 12:23 AM
> Hey fellas:
>
> I would like to integrate the following function:
>
> integrand <- function (x,t) {
> exp(-2*t)*(2*t)^x/(10*factorial(x))
> }
>
> with respect to the t variable, from 0 to 10.
> The variable x here works as a parameter: I would like to
> integrate the said function for each value of x in
> 0,1,..,44.
>
> I have tried Vectorize to no avail.
>
> Thanks in advance,
> jose romero
>
> __
> R-help@r-project.org mailing list
> 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.
__
R-help@r-project.org mailing list
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.
Re: [R] Integrate a 1-variable function with 1 parameter (Jose L. Romero)
To Ravi Varadhan and Peter Dalgaard:
Infinite thanks for your help and suggestions- they were very instructive, as i
have still to learn on the appropriate use of the apply family of functions.
jlrp
--- On Wed, 8/27/08, Ravi Varadhan <[EMAIL PROTECTED]> wrote:
> From: Ravi Varadhan <[EMAIL PROTECTED]>
> Subject: RE: [R] Integrate a 1-variable function with 1 parameter (Jose L.
> Romero)
> To: [EMAIL PROTECTED], r-help@r-project.org
> Date: Wednesday, August 27, 2008, 2:42 PM
> Here is one way:
>
> integrand <- function (t, x) {
> exp(-2*t)*(2*t)^x/(10*factorial(x))
> }
>
> x <- 0:44
> ans <- sapply(x, function(x) integrate(integrand,
> lower=0, upper=10, x=x))
> cbind(x=x, integral=unlist(ans[1,]),
> abs.error=unlist(ans[2,]))
>
>
> Ravi.
>
>
> ---
>
> Ravi Varadhan, Ph.D.
>
> Assistant Professor, The Center on Aging and Health
>
> Division of Geriatric Medicine and Gerontology
>
> Johns Hopkins University
>
> Ph: (410) 502-2619
>
> Fax: (410) 614-9625
>
> Email: [EMAIL PROTECTED]
>
> Webpage:
> http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
>
>
>
>
> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On
> Behalf Of jose romero
> Sent: Wednesday, August 27, 2008 10:24 AM
> To: r-help@r-project.org
> Subject: [R] Integrate a 1-variable function with 1
> parameter (Jose L.
> Romero)
>
> Hey fellas:
>
> I would like to integrate the following function:
>
> integrand <- function (x,t) {
> exp(-2*t)*(2*t)^x/(10*factorial(x))
> }
>
> with respect to the t variable, from 0 to 10.
> The variable x here works as a parameter: I would like to
> integrate the said
> function for each value of x in 0,1,..,44.
>
> I have tried Vectorize to no avail.
>
> Thanks in advance,
> jose romero
>
> __
> R-help@r-project.org mailing list
> 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.
__
R-help@r-project.org mailing list
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.
Re: [R] Integrate a 1-variable function with 1 parameter (Jose L. Romero)
Peter Dalgaard wrote:
> jose romero wrote:
>
>> Hey fellas:
>>
>> I would like to integrate the following function:
>>
>> integrand <- function (x,t) {
>> exp(-2*t)*(2*t)^x/(10*factorial(x))
>> }
>>
>> with respect to the t variable, from 0 to 10.
>> The variable x here works as a parameter: I would like to integrate the said
>> function for each value of x in 0,1,..,44.
>>
>> I have tried Vectorize to no avail.
>>
>>
>>
> Will this do?
>
> sapply(0:44,function(x)integrate(integrand,x=x, lower=0, upper=10)$value
>
>
Oops. Forgot this:
)
--
O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
__
R-help@r-project.org mailing list
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.
Re: [R] Integrate a 1-variable function with 1 parameter (Jose L. Romero)
Here is one way:
integrand <- function (t, x) {
exp(-2*t)*(2*t)^x/(10*factorial(x))
}
x <- 0:44
ans <- sapply(x, function(x) integrate(integrand, lower=0, upper=10, x=x))
cbind(x=x, integral=unlist(ans[1,]), abs.error=unlist(ans[2,]))
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of jose romero
Sent: Wednesday, August 27, 2008 10:24 AM
To: r-help@r-project.org
Subject: [R] Integrate a 1-variable function with 1 parameter (Jose L.
Romero)
Hey fellas:
I would like to integrate the following function:
integrand <- function (x,t) {
exp(-2*t)*(2*t)^x/(10*factorial(x))
}
with respect to the t variable, from 0 to 10.
The variable x here works as a parameter: I would like to integrate the said
function for each value of x in 0,1,..,44.
I have tried Vectorize to no avail.
Thanks in advance,
jose romero
__
R-help@r-project.org mailing list
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.
__
R-help@r-project.org mailing list
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.
Re: [R] Integrate a 1-variable function with 1 parameter (Jose L. Romero)
jose romero wrote:
> Hey fellas:
>
> I would like to integrate the following function:
>
> integrand <- function (x,t) {
> exp(-2*t)*(2*t)^x/(10*factorial(x))
> }
>
> with respect to the t variable, from 0 to 10.
> The variable x here works as a parameter: I would like to integrate the said
> function for each value of x in 0,1,..,44.
>
> I have tried Vectorize to no avail.
>
>
Will this do?
sapply(0:44,function(x)integrate(integrand,x=x, lower=0, upper=10)$value
--
O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
__
R-help@r-project.org mailing list
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.
[R] Integrate a 1-variable function with 1 parameter (Jose L. Romero)
Hey fellas:
I would like to integrate the following function:
integrand <- function (x,t) {
exp(-2*t)*(2*t)^x/(10*factorial(x))
}
with respect to the t variable, from 0 to 10.
The variable x here works as a parameter: I would like to integrate the said
function for each value of x in 0,1,..,44.
I have tried Vectorize to no avail.
Thanks in advance,
jose romero
__
R-help@r-project.org mailing list
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.
Re: [R] Integrate() error message, I am at a loss
On Wed, 12 Sep 2007, Sergey Goriatchev wrote:
> Hello!
>
> I have a problem with integrate() in my function nctspa(). Integrate
> produces an error message "evaluation of function gave a result of
> wrong length". I don't know what that means. Could anyone suggest me
> what is wrong with my function?
Sure, but you can do this yourself.
For one thing, setting
options( error = recover )
before you run your function will help you to see what is happening.
(viz. after the error message select 2, then figure out what f() and ff()
are, then try ff(1))
>
> These are the examples of function calls that work OK:
> nctspa(a=1:10,n=5)
> nctspa(a=1:10, n=5, mu=2, theta=3, renorm=0)
>
> This does not work:
> nctspa(a=1:10, n=5, mu=2, theta=3, renorm=1)
>
> Many thanks in advance for your help!
> please, send reply also to [EMAIL PROTECTED]
>
> /Sergey
>
> Here is the function:
>
> #Computes the s.p.a. to the doubly noncentral t at value x.
> #degrees of freedom n, noncentrality parameters mu and theta.
> #==#
> nctspa <- function(a,n,mu=0,theta=0,renorm=0,rec=0){
> #Pass renorm=1 to renormalize the SPA to the pdf.
> #There is a last argument called rec. DO NOT PASS it!
>
> alpha <- mu/sqrt((1+theta/n))
> normconst <- 1
> if(renorm==1 & rec==0){
> term1 <- integrate(nctspa, -Inf, alpha, n=n, mu=mu, theta=theta)$value
> term2 <- integrate(nctspa, alpha, Inf, n=n, mu=mu, theta=theta)$value
Whoa! Let's read the help page for integrate:
Arguments:
f: an R function taking a numeric first argument and returning a
numeric vector of the same length.
..^.
which is not what nctspa() does. It returns a list of length 2 not matter
what the first arguement.
Did you mean something like
integrate(function(x,...) nctspa(x,...)$PDF, ??
HTH,
Chuck
> normconst <- 1/(term1+term2)
> }
> cdf <- numeric()
> pdf <- cdf
> c3 <- n^2+2*n*a^2+a^4
> c2 <- (-2*mu*(a^3+n*a))/c3
> c1 <- (-n^2-n*a^2-n*theta+a^2*mu^2)/c3
> c0 <- (n*a*mu)/c3
> q <- c1/3-(c2^2)/9
> r <- 1/6*(c1*c2-3*c0)-1/27*c2^3
> b0 <- sqrt(-4*q)*cos(acos(r/sqrt(-q^3))/3)-c2/3
> t1 <- -mu+a*b0
> t2 <- -a*t1/b0/n/2
> nu <- 1/(1-2*t2)
> w <- sqrt(-mu*t1-n*log(nu)-2*theta*nu*t2)*sign(a-alpha)
> u <- sqrt((a^2+2*n*t2)*(2*n*nu^2+4*theta*nu^3)+4*n^2*b0^2)/(2*n*b0^2)
> pdf <- normconst*dnorm(w)/u
>
> nz <- (abs(t1*b0)>=1e-10)
> iz <- (abs(t1*b0)<=1e-10)
> if(any(nz)){
> d <- numeric()
> d[nz] <- 1/(t1[nz]*b0[nz])
> cdf[nz] <- pnorm(w[nz])+dnorm(w[nz])*(1/w[nz]-d[nz]/u[nz])
> }
> if(any(iz)){
> n <- sum(iz==1)
> rez <- nctspa(c(a[iz]-1e-4, a[iz]+1e-4),n,mu,theta,0,rec+1)
> if(rec>5){
> cdf[iz] <- 0.5
> warning('Too many recursions')
> } else {
> cdf[iz] <- 0.5*(rez$CDF[1:n]+rez$CDF[(n+1):length(rez$CDF)])
> }
> }
> list(PDF=pdf, CDF=cdf)
> }
> #==
>
> __
> R-help@r-project.org mailing list
> 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.
>
Charles C. Berry(858) 534-2098
Dept of Family/Preventive Medicine
E mailto:[EMAIL PROTECTED] UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901
__
R-help@r-project.org mailing list
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.
[R] Integrate() error message, I am at a loss
Hello!
I have a problem with integrate() in my function nctspa(). Integrate
produces an error message "evaluation of function gave a result of
wrong length". I don't know what that means. Could anyone suggest me
what is wrong with my function?
These are the examples of function calls that work OK:
nctspa(a=1:10,n=5)
nctspa(a=1:10, n=5, mu=2, theta=3, renorm=0)
This does not work:
nctspa(a=1:10, n=5, mu=2, theta=3, renorm=1)
Many thanks in advance for your help!
please, send reply also to [EMAIL PROTECTED]
/Sergey
Here is the function:
#Computes the s.p.a. to the doubly noncentral t at value x.
#degrees of freedom n, noncentrality parameters mu and theta.
#==#
nctspa <- function(a,n,mu=0,theta=0,renorm=0,rec=0){
#Pass renorm=1 to renormalize the SPA to the pdf.
#There is a last argument called rec. DO NOT PASS it!
alpha <- mu/sqrt((1+theta/n))
normconst <- 1
if(renorm==1 & rec==0){
term1 <- integrate(nctspa, -Inf, alpha, n=n, mu=mu, theta=theta)$value
term2 <- integrate(nctspa, alpha, Inf, n=n, mu=mu, theta=theta)$value
normconst <- 1/(term1+term2)
}
cdf <- numeric()
pdf <- cdf
c3 <- n^2+2*n*a^2+a^4
c2 <- (-2*mu*(a^3+n*a))/c3
c1 <- (-n^2-n*a^2-n*theta+a^2*mu^2)/c3
c0 <- (n*a*mu)/c3
q <- c1/3-(c2^2)/9
r <- 1/6*(c1*c2-3*c0)-1/27*c2^3
b0 <- sqrt(-4*q)*cos(acos(r/sqrt(-q^3))/3)-c2/3
t1 <- -mu+a*b0
t2 <- -a*t1/b0/n/2
nu <- 1/(1-2*t2)
w <- sqrt(-mu*t1-n*log(nu)-2*theta*nu*t2)*sign(a-alpha)
u <- sqrt((a^2+2*n*t2)*(2*n*nu^2+4*theta*nu^3)+4*n^2*b0^2)/(2*n*b0^2)
pdf <- normconst*dnorm(w)/u
nz <- (abs(t1*b0)>=1e-10)
iz <- (abs(t1*b0)<=1e-10)
if(any(nz)){
d <- numeric()
d[nz] <- 1/(t1[nz]*b0[nz])
cdf[nz] <- pnorm(w[nz])+dnorm(w[nz])*(1/w[nz]-d[nz]/u[nz])
}
if(any(iz)){
n <- sum(iz==1)
rez <- nctspa(c(a[iz]-1e-4, a[iz]+1e-4),n,mu,theta,0,rec+1)
if(rec>5){
cdf[iz] <- 0.5
warning('Too many recursions')
} else {
cdf[iz] <- 0.5*(rez$CDF[1:n]+rez$CDF[(n+1):length(rez$CDF)])
}
}
list(PDF=pdf, CDF=cdf)
}
#==
__
R-help@r-project.org mailing list
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.
