Re: [R] R's integrate function
I got the same error.
sessionInfo()
R version 2.6.0 (2007-10-03)
i386-pc-mingw32
locale:
LC_COLLATE=Spanish_Colombia.1252;LC_CTYPE=Spanish_Colombia.1252;LC_MONETARY=Spanish_Colombia.1252;LC_NUMERIC=C;LC_TIME=Spanish_Colombia.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
however, you can modify the maximum number of subdivisions for
integrate (check ?integrate)
z=423
integrate(function(y,z){
+ sapply(y, function(y,z){
+ integrate(function(x,z)
+
1/x*dbeta(0.01,x/(0.005/1.005),(1-x)/(0.005/1.005))*dbeta(y,x/(0.005/1.005),(1-x)/(0.005/1.005))*(1-y)^z,0,1,423)$value
})
+ },0,1,423, subdivisions = 150)
18.97078 with absolute error 0.0013
hth,
Víctor H
2008/9/30 Charles C. Berry [EMAIL PROTECTED]:
What verson of R?
Works for me:
integrate(function(y,z){
+ sapply(y, function(y,z){
+ integrate(function(x,z)
+
1/x*dbeta(0.01,x/(0.005/1.005),(1-x)/(0.005/1.005))*dbeta(y,x/(0.005/1.005),(1-x)/(0.005/1.005))*(1-y)^z,0,1
+ ,423)$value
+ })
+ },0,1,423)
18.9513 with absolute error 0.0011
sessionInfo()
R version 2.7.2 (2008-08-25)
i386-apple-darwin8.11.1
locale:
C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
HTH,
Chuck
On Tue, 30 Sep 2008, Susanne Pfeifer wrote:
Hello,
I am trying to use R's integrate function to calculate the following
integral for z=423:
integrate(function(y,z){
sapply(y, function(y,z){
integrate(function(x,z)
1/x*dbeta(0.01,x/(0.005/1.005),(1-x)/(0.005/1.005))*dbeta(y,x/(0.005/1.005),(1-x)/(0.005/1.005))*(1-y)^z,0,1,423)$value
})
},0,1,423)$value
but I receive an error message saying that the maximum number of
subdivisions is reached. If I choose a smaller z it works fine. Does
anyone know what the problem is?
Thank you very much in advance,
Tiffy
__
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-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] Is there a way to not use an explicit loop?
Hi,
you might try this:
set.seed(100)
m - 10
size.a- 10
prob.a- 0.3
prior.constant = 0
draw1 = rbinom( m , size.a, prob.a )
beta.draws - function(draw, size.a, prior.constant, n) {
rbeta(n, prior.constant + draw, prior.constant + size.a - draw)
}
bdraws - sapply(draw1, beta.draws, size.a = size.a, prior.constant =
prior.constant, n = 1)
beta.post - apply(bdraws, 2, function(x) c(post.mean = mean(x),
post.median = median(x)) )
beta.post
[,1] [,2] [,3] [,4] [,5]
[,6] [,7] [,8]
post.mean 0.2017118 0.1996809 0.2991173 0.10069613 0.3001924
0.2991149 0.4033310 0.2003104
post.median 0.1804893 0.1791630 0.2845427 0.07505278 0.2858155
0.2844503 0.3961419 0.1790511
[,9] [,10]
post.mean 0.3013020 0.1990232
post.median 0.2886199 0.1786447
best
Víctor H Cervantes
2008/9/17 Juancarlos Laguardia [EMAIL PROTECTED]:
I have a problem in where i generate m independent draws from a binomial
distribution,
say
draw1 = rbinom( m , size.a, prob.a )
then I need to use each draw to generate a beta distribution. So, like
using a beta prior, binomial likelihood, and obtain beta posterior, m many
times. I have not found out a way to vectorize draws from a beta
distribution, so I have an explicit for loop within my code
for( i in 1: m ) {
beta.post = rbeta( 1, draw1[i] + prior.constant , prior.constant +
size.a - draw1[i] )
beta.post.mean[i] = mean(beta.post)
beta.post.median[i] = median(beta.post)
etc.. for other info
}
Is there a way to vectorize draws from an beta distribution?
UC Slug
__
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] psychometric functions
Hi, 2008/8/21 Mario Maiworm [EMAIL PROTECTED] Hi, I want to fit some psychophysical data with cumulative gaussians. There is quite a convenient toolbox for matlab called 'psignifit' (formerly known as 'psychofit'). It allows the lower bound of the sigmoid to vary slightly from zero, aswell as the upper bound to vary from one. with these two free parameters, the fitted function is less sensitive to noisy data and outliers. you might also check the function psyfun.2asym from the psyphy package http://cran.r-project.org/web/packages/psyphy/index.html best Víctor H Apart from advertising this toolbox I want to ask for possibilities in R to fit psychometric functions, as I would rather use R than matlab. Is there a comparable package specific for psychophysics in R? otherwise: which function would be a good choice? Can I have functionality equal to 'psignifit' with glm() from MASS. I would be grateful for some suggestions of people who have experience in sigmoid-fitting in R. Best, Mario. __ Mario Maiworm Biological Psychology and Neuropsychology University of Hamburg Von-Melle-Park 11 D-20146 Hamburg Phone: +49 40 42838 8265 Fax: +49 40 42838 6591 http://bpn.uni-hamburg.de/Maiworm_e.html http://cinacs.org __ 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.
