[R] How to double integrate a function in R
I would like to express my gratitude for the great help given by David and
Hans regarding my last query.
Thank you very much for your time, people.
All the best,
Tiago
---
Hello, R users!
I am trying to double integrate the following expression:
# expression
(1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1)))
for y2>y1>0.
I am trying the following approach
# first attempt
library(cubature)
fun <- function(x) {
(1/(2*pi))*exp(-x[2]/2)*sqrt((x[1]/(x[2]-x[1])))} adaptIntegrate(fun,
lower = c(0,0), upper =c(5, 6), tol=1e-8)
However, I don't know how to constrain the integration so that y2>y1>0.
Any ideas?
Tiago
--
Tiago V. Pereira, MSc, PhD
Center for Studies of the Human Genome
Department of Genetics and Evolutionary Biology
University of São Paulo
Rua do Matão, 277
CEP 05508-900
São Paulo - SP, Brazil
__
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] How to double integrate a function in R
Hello, R users!
I am trying to double integrate the following expression:
# expression
(1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1)))
for y2>y1>0.
I am trying the following approach
# first attempt
library(cubature)
fun <- function(x) { (1/(2*pi))*exp(-x[2]/2)*sqrt((x[1]/(x[2]-x[1])))}
adaptIntegrate(fun, lower = c(0,0), upper =c(5, 6), tol=1e-8)
However, I don't know how to constrain the integration so that y2>y1>0.
Any ideas?
Tiago
__
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] How to compute a P-value for a complex mixture of chi-squared distributions in R
Thank you very much, Rui and Peter, for you detailed and helpful tips! It worked like a charm! I would spend more two weeks (or more) to figure out that by myself. Cheers! Tiago > Hello, > > No, nothing wrong. (I feel silly for not having noticed it.) In fact not > only it's much simpler but it's also more accurate than the use of > accurate with the default rel.tol. > It should be better, however, to use lower.tail = FALSE, since the op > wants p-values. > > 0.5 * pchisq(x^2, 1, lower.tail = FALSE) + 0.5 * pchisq(x^2, 2, > lower.tail = FALSE) > > Rui Barradas > > Em 01-06-2013 14:57, peter dalgaard escreveu: >> >> On Jun 1, 2013, at 06:32 , Tiago V. Pereira wrote: >> >>> Hello, R users! >>> >>> I am struggling with the following problem: >>> >>> I need to compute a P-value for a mixture of two chi-squared >>> distributions. My P-value is given by: >>> >>> P = 0.5*prob(sqrt(chi2(1)) <= x) + 0.5*prob(sqrt(chi2(2)) <= x) >>> >>> In words, I need to compute the p-value for 50?50 mixture of the square >>> root of a chi-squared random variable with 1 degree of freedom and the >>> square root of a chi-squared with two degrees of freedom. >>> >>> Although I can quickly simulate data, the P-values I am looking for are >>> at >>> the tail of the distribution, that is, alpha levels below 10^-7. Hence, >>> simulation is not efficient. >>> >>> Are you aware of smart approach? >> >> Er,... >> >> Anything wrong with >> >> 0.5 * pchisq(x^2, 1) + 0.5 * pchisq(x^2, 2) >> >> ??? >> >> -pd >> >> >>> >>> >>> All the best, >>> >>> Tiago >>> >>> __ >>> 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] How to compute a P-value for a complex mixture of chi-squared distributions in R
Hello, R users! I am struggling with the following problem: I need to compute a P-value for a mixture of two chi-squared distributions. My P-value is given by: P = 0.5*prob(sqrt(chi2(1)) <= x) + 0.5*prob(sqrt(chi2(2)) <= x) In words, I need to compute the p-value for 5050 mixture of the square root of a chi-squared random variable with 1 degree of freedom and the square root of a chi-squared with two degrees of freedom. Although I can quickly simulate data, the P-values I am looking for are at the tail of the distribution, that is, alpha levels below 10^-7. Hence, simulation is not efficient. Are you aware of smart approach? All the best, Tiago __ 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] How to compute a P-value for a complex mixture of chi-squared distributions in R
Hello, R users! I am struggling with the following problem: I need to compute a P-value for a mixture of two chi-squared distributions. My P-value is given by: P = 0.5*prob(sqrt(chi2(1)) <= x) + 0.5*prob(sqrt(chi2(2)) <= x) In words, I need to compute the p-value for 5050 mixture of the square root of a chi-squared random variable with 1 degree of freedom and the square root of a chi-squared with two degrees of freedom. Although I can quickly simulate data, the P-values I am looking for are at the tail of the distribution, that is, alpha levels below 10^-7. Hence, simulation is not efficient. Are you aware of smart approach? All the best, Tiago __ 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] How to create a correct matrix in R
Hello Rlisters!
In my codes, I need to import a matrix:
v <- read.table("/home/tiago/matrix.txt", header=FALSE)
v<-as.matrix(v)
v
V1 V2 V3 V4 V5 V6
[1,] 1. -0.89847480 -0.73929292 -0.99055335 -0.04514469 0.04056137
[2,] -0.89847480 1. 0.95986852 0.82978466 0.04056137 -0.04514469
[3,] -0.73929292 0.95986852 1. 0.63996937 0.03337515 -0.04333297
[4,] -0.99055335 0.82978466 0.63996937 1. 0.04471823 -0.03746038
[5,] -0.04514469 0.04056137 0.03337515 0.04471823 1. -0.89847480
[6,] 0.04056137 -0.04514469 -0.04333297 -0.03746038 -0.89847480 1.
[7,] -0.60519045 0.67357531 0.64654374 0.55892246 -0.06244832 0.06950480
V7
[1,] -0.60519045
[2,] 0.67357531
[3,] 0.64654374
[4,] 0.55892246
[5,] -0.06244832
[6,] 0.06950480
[7,] 1.
However, I keep getting the same error after loading that matrix:
`v' is not a covariance matrix
Nonetheless, if I input the matrix directly, there is no error:
x1 = c(1, -0.898474804259413, -0.739292919198965, -0.990553354617789,
-0.0451446949071635, 0.0405613709200646, -0.605190448449146)
x2 = c(-0.89847480425931, 1, 0.959868518981255, 0.829784658203916,
0.0405613709200599, -0.0451446949071635, 0.673575314054563)
x3 = c(-0.739292919198939, 0.959868518981239, 1, 0.639969373426519,
0.0333751532842623, -0.0433329714403989, 0.646543739123876)
x4 = c(-0.990553354617685, 0.82978465820392, 0.639969373426531, 1,
0.0447182289834827, -0.0374603752332609, 0.558922461747364)
x5 = c(-0.0451446949071635, 0.0405613709200646, 0.0333751532842635,
0.0447182289834874, 1, -0.898474804259413, -0.0624483157850655)
x6 =
c(0.0405613709200679, -0.0451446949071635,
-0.0433329714403999,-0.0374603752332612,-0.898474804259486,1,0.0695048046856916)
x7 = c(-0.605190448449077, 0.673575314054563, 0.646543739123887,
0.558922461747361, -0.0624483157850583, 0.0695048046856916, 1)
v <- rbind(x1,x2,x3,x4,x5,x6,x7)
row.names(v)<-NULL
v
[,1][,2][,3][,4][,5][,6]
[1,] 1. -0.89847480 -0.73929292 -0.99055335 -0.04514469 0.04056137
[2,] -0.89847480 1. 0.95986852 0.82978466 0.04056137 -0.04514469
[3,] -0.73929292 0.95986852 1. 0.63996937 0.03337515 -0.04333297
[4,] -0.99055335 0.82978466 0.63996937 1. 0.04471823 -0.03746038
[5,] -0.04514469 0.04056137 0.03337515 0.04471823 1. -0.89847480
[6,] 0.04056137 -0.04514469 -0.04333297 -0.03746038 -0.89847480 1.
[7,] -0.60519045 0.67357531 0.64654374 0.55892246 -0.06244832 0.06950480
[,7]
[1,] -0.60519045
[2,] 0.67357531
[3,] 0.64654374
[4,] 0.55892246
[5,] -0.06244832
[6,] 0.06950480
[7,] 1.
How can one import the data correctly?
I could not figure it out.
Thanks in advance.
Tiago
__
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] How can I access information stored after I run a command in R?
I would like to thank Justin and Matthias for their very helpful on my
query (see it below).
All the best,
Tiago
Dear all,
Supposed I run the following command:
###
#install.packages("Rassoc", dependencies=TRUE)
library("Rassoc")
ca=c(139,249,112)
co=c(136,244,120)
a=rbind(ca,co)
MAX3(a,"asy",1)
##
I get:
The MAX3 test using the asy method
data: a
statistic = 0.5993, p-value = 0.7933
How can one save the result 0.7933 into a file?
say:
foo <- 0.7933
write.table(foo, file ="/home/foo.txt", sep = " ",
row.names=FALSE,col.names=TRUE, quote=FALSE, qmethod = "double")
However, instead of typing the value above, I would like to replace it by
the macro (scalar, local) that has the accurate p-value.
thanks in advance for your help.
Tiago
__
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] How can I access information stored after I run a command in R?
Dear all,
Supposed I run the following command:
###
#install.packages("Rassoc", dependencies=TRUE)
library("Rassoc")
ca=c(139,249,112)
co=c(136,244,120)
a=rbind(ca,co)
MAX3(a,"asy",1)
##
I get:
The MAX3 test using the asy method
data: a
statistic = 0.5993, p-value = 0.7933
How can one save the result 0.7933 into a file?
say:
foo <- 0.7933
write.table(foo, file ="/home/foo.txt", sep = " ",
row.names=FALSE,col.names=TRUE, quote=FALSE, qmethod = "double")
However, instead of typing the value above, I would like to replace it by
the macro (scalar, local) that has the accurate p-value.
thanks in advance for your help.
Tiago
__
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.
