Re: [R] two dimensional integration
On 16-02-2013, at 18:01, julia cafnik julia.caf...@gmail.com wrote:
Dear R-users,
I'm wondering how to calculate this double integral in R:
int_a^b int_c^y g(x, y) dx dy
where g(x,y) = exp(- alpha (y - x)) * b
A very similar question was asked about nine years ago:
http://tolstoy.newcastle.edu.au/R/help/04/10/5831.html
The final message in the thread gave a workable answer.
In your case define function g (leave out the multiply by b since you can
always do that outside the integral).
g - function(x,y,alpha=1) exp(-alpha*(y-x))
Assume a=0, b=1 and c=0.
Then following the final message in the thread mentioned above there are two
ways of getting an answer:
if function g is fully vectorized:
integrate( function(y) {
sapply(y, function(y) {
integrate(function(x) g(x,y),0,y)$value
})
},0,1)
and if g is not vectorized and only takes scalars as input
integrate(function(y) {
sapply(y, function(y) {
integrate(function(x) {
sapply(x, function(x) g(x,y)) },0,y)$value
})
},0,1)
Both of the approaches result in
0.3678794 with absolute error 4.1e-15
which corresponds to the exact analytical answer 1/exp(1) (if I didn't make a
mistake)
You can also do this with package cubature with Gabor's approach (from the
mentioned thread) like this
library(cubature)
h - function(z) g(z[1],z[2])*(z[1]z[2])
adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4)
but it's also very slow with higher tolerances.
adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4)
$integral
[1] 0.3678723
$error
[1] 3.677682e-05
$functionEvaluations
[1] 268413
$returnCode
[1] 0
Berend
__
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] two dimensional integration
thank for your help. already solved it.
Cheers,
J.
On Sun, Feb 17, 2013 at 9:41 AM, Berend Hasselman b...@xs4all.nl wrote:
On 16-02-2013, at 18:01, julia cafnik julia.caf...@gmail.com wrote:
Dear R-users,
I'm wondering how to calculate this double integral in R:
int_a^b int_c^y g(x, y) dx dy
where g(x,y) = exp(- alpha (y - x)) * b
A very similar question was asked about nine years ago:
http://tolstoy.newcastle.edu.au/R/help/04/10/5831.html
The final message in the thread gave a workable answer.
In your case define function g (leave out the multiply by b since you can
always do that outside the integral).
g - function(x,y,alpha=1) exp(-alpha*(y-x))
Assume a=0, b=1 and c=0.
Then following the final message in the thread mentioned above there are
two ways of getting an answer:
if function g is fully vectorized:
integrate( function(y) {
sapply(y, function(y) {
integrate(function(x) g(x,y),0,y)$value
})
},0,1)
and if g is not vectorized and only takes scalars as input
integrate(function(y) {
sapply(y, function(y) {
integrate(function(x) {
sapply(x, function(x) g(x,y)) },0,y)$value
})
},0,1)
Both of the approaches result in
0.3678794 with absolute error 4.1e-15
which corresponds to the exact analytical answer 1/exp(1) (if I didn't
make a mistake)
You can also do this with package cubature with Gabor's approach (from the
mentioned thread) like this
library(cubature)
h - function(z) g(z[1],z[2])*(z[1]z[2])
adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4)
but it's also very slow with higher tolerances.
adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4)
$integral
[1] 0.3678723
$error
[1] 3.677682e-05
$functionEvaluations
[1] 268413
$returnCode
[1] 0
Berend
[[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] two dimensional integration
On 17-02-2013, at 10:01, julia cafnik julia.caf...@gmail.com wrote: thank for your help. already solved it. Show us how. So that others looking for answers to similar problems in future can find an answer. Berend __ 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] two dimensional integration
This function is separable.
If calculating by hand with bounds a -0 and b - 1 i got the result:
theta / kappa * ( 1 + exp( - kappa) / kappa - 1 / kappa)
by putting
theta - 0.1
kappa -0.3
in the above result I got 0.04535
I implemented it in R this way:
integrate(function(y) {
sapply(y, function(y) {
integrate(function(x) exp(- kappa * (y - x)) * theta, a, y)$value
})
}, a, b)$value
)
Result: [1] 0.04535358
Thanks for helping!
On Sun, Feb 17, 2013 at 10:04 AM, Berend Hasselman b...@xs4all.nl wrote:
On 17-02-2013, at 10:01, julia cafnik julia.caf...@gmail.com wrote:
thank for your help. already solved it.
Show us how.
So that others looking for answers to similar problems in future can find
an answer.
Berend
[[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] Multidimensional correlation matrix question
Hello, I previously sent this message which was stripped off due to HTML mail. Sorry! Thanks. Hello All, I am new to the list. I have been learning to use R recently and its been great to see so much help available. However I must admit, I have stumbled upon a problem with correlation matrices and was hoping if someone could help. I have an excel workbook with a sheet per drug. Each sheet has 40 patients' drug response measured over different time points (several days). I have already been able to create a cor() matrix of all the drugs for all the patients on each day. However I have been asked to prepare a matrix of all days of all people for all drugs. So since there are 26 drugs, I am expected to form a 26x26 matrix for all this data. I am not even sure if this is possible. Can anyone point me in the right direction? PS: If I make individual matrix and combine them using cbind(), it is not a 26x26 sqaure matrix. Any help would be greatly appreciated. Thanks __ 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] Getting WinBUGS Leuk example to work from R using R2winBUGS
I am trying to learn to use winBUGS from R, I have experience with R.
I have managed to successfully run a simple example from R with no
problems. I have been trying to run the Leuk: Survival from winBUGS
examples Volume 1. I have managed to run this from winBUGS GUI with no
problems. My problem is try as I might ( and I have been trying and
searching for days) I cannot get it to run using R2winBUGS.I am sure
it is something simple.
The error message I get if I try and set inits in the script is
Error in bugs(data = L, inits = inits,
parameters.to.save = params, model.file model.txt, :
Number of initialized chains (length(inits)) != n.chains
I know this means I have not initialised some of the chains, but I am
pasting the inits code from winbugs examples manual and all other
settings seem to me to be the same as when run on the winBUGS GUI.
If I try inits=NULL I get another error message
display(log)
check(C:/BUGS/model.txt)
model is syntactically correct
data(C:/BUGS/data.txt)
data loaded
compile(1)
model compiled
gen.inits()
shape parameter (r) of gamma dL0[1] too small -- cannot sample
thin.updater(1)
update(500)
command #Bugs:update cannot be executed (is greyed out)
set(beta)
Which indicates to me I will still have problems after solving the
first one!! I am about to give up on using winBUGS, please can someone
save me? I know I am probably going to look stupid, but everyone has
to learn:-)
I have also tried changing nc-2 (On advice, which doesnt work and
gives an uninitialised chain error)
I am using winBUGS 1.4.3 on Windows XP 2002 SP3
My R code is below, many thanks for at least reading this far.
rm(list = ls())
L-list(N = 42, T = 17, eps = 1.0E-10,
obs.t = c(1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12,
12, 15, 17, 22, 23, 6,
6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 32,
32, 34, 35),
fail = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0,
0),
Z = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
-0.5),
t = c(1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 22, 23, 35))
### 5.4. Analysis using WinBUGS
library(R2WinBUGS) # Load the R2WinBUGS library CHOOSE to use WinBUGS
#library(R2OpenBUGS)# Load the R2OpenBUGS library CHOOSE
to use OpenBUGS
setwd(C://BUGS)
# Save BUGS description of the model to working directory
sink(model.txt)
cat(
model
{
# Set up data
for(i in 1:N) {
for(j in 1:T) {
# risk set = 1 if obs.t = t
Y[i,j] - step(obs.t[i] - t[j] + eps)
# counting process jump = 1 if obs.t in [ t[j], t[j+1] )
# i.e. if t[j] = obs.t t[j+1]
dN[i, j] - Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * fail[i]
}
}
# Model
for(j in 1:T) {
for(i in 1:N) {
dN[i, j] ~ dpois(Idt[i, j]) # Likelihood
Idt[i, j] - Y[i, j] * exp(beta * Z[i]) * dL0[j] # Intensity
}
dL0[j] ~ dgamma(mu[j], c)
mu[j] - dL0.star[j] * c # prior mean hazard
# Survivor function = exp(-Integral{l0(u)du})^exp(beta*z)
S.treat[j] - pow(exp(-sum(dL0[1 : j])), exp(beta * -0.5));
S.placebo[j] - pow(exp(-sum(dL0[1 : j])), exp(beta * 0.5));
}
c - 0.001
r - 0.1
for (j in 1 : T) {
dL0.star[j] - r * (t[j + 1] - t[j])
}
beta ~ dnorm(0.0,0.01)
}
,fill=TRUE)
sink()
params- c(beta,S.placebo,S.treat)
inits-list( beta = 0.0,
dL0 = c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0))
# MCMC settings
nc -1 # Number of chains
ni - 1000 # Number of draws from posterior (for each chain)
ns-1000 #Number of sims (n.sims)
nb - floor(ni/2) # Number of draws to discard as burn-in
nt - max(1, floor(nc * (ni - nb) / ns))# Thinning rate
Lout-list()
# Start Gibbs sampler: Run model in WinBUGS and save results in object
called out
out - bugs(data = L, inits =inits, parameters.to.save = params,
model.file = model.txt,
n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni, debug = T, DIC
= TRUE,digits=5,
codaPkg=FALSE, working.directory = getwd())
__
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] Extracting Numeric Columns from Data Fram
Dear Barry, I saw that you received several nice answers on how to 'pull out' numeric columns. You also wrote that an alternative could be to 'just operate only on' numerics. Here is one possibility: library(plyr) # some dummy data df - data.frame(nonnum1 = letters[1:5], num1 = 1:5, nonnum2 = letters[6:10], num2 = rnorm(5)) # operate only on numeric variables numcolwise(.fun = mean)(df) # or only on 'discrete' variables (following the terminology on ?colwise) catcolwise(.fun = toupper)(df) Best regards, Henrik Hello, I've got a data frame with a mix of numeric, integer and factor columns. I'd like to pull out (or just operate only on) the numeric/integer columns. Every thing I've found in searches is about how to subset by rows, or how to operate assuming you have the column names.? I'd like to pull by type. Thanks! Barry -- Henrik Pärn Centre for Conservation Biology Department of Biology Norwegian University of Science and Technology NO-7491 Trondheim NORWAY Office: +47 735 96084 Mobile: +47 909 89 255 Fax: +47 735 96100 __ 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] reading data
HI Vera,
No problem. I am cc:ing to r-help.
A.K.
From: Vera Costa veracosta...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Sunday, February 17, 2013 5:44 AM
Subject: Re: reading data
Hi. Thank you. It works now:-)
And yes, I use windows.
Thank you very much.
No dia 17 de Fev de 2013 00:44, arun smartpink...@yahoo.com escreveu:
Hi Vera,
Have you tried the suggestion?
Are you using Windows?
Thanks,
Arun
From: Vera Costa veracosta...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Saturday, February 16, 2013 7:10 PM
Subject: Re: reading data
Thank you.
In mine, I have an error 'what' must be a character string or a function.
I need to do equivalent in my system.
Thank you and sorry one more time.
No dia 16 de Fev de 2013 23:53, arun smartpink...@yahoo.com escreveu:
Hi,
You didn't mention what the error message or whether you are reading file
names which are not m11kk.txt.
It is workiing on my system as I run it again.
?c() combine values into a vector or list.
sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_CA.UTF-8 LC_COLLATE=en_CA.UTF-8
[5] LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8
[7] LC_PAPER=C LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] stringr_0.6.2 reshape2_1.2.2
loaded via a namespace (and not attached):
[1] plyr_1.8
#code
res-do.call(c,lapply(list.files(recursive=T)[grep(m11kk,list.files(recursive=T))],function(x)
{names(x)-gsub(^(.*)\\/.*,\\1,x); lapply(x,function(y)
read.table(y,header=TRUE,stringsAsFactors=FALSE,fill=TRUE))})) #it seems
like one of the rows of your file doesn't have 6 elements, so added fill=TRUE
names(res)-paste(group_,gsub(\\d+,,names(res)),sep=)
res2-split(res,names(res))
res3- lapply(res2,function(x)
{names(x)-paste(gsub(.*_,,names(x)),1:length(x),sep=);x})
#result
res3
#$group_a
#$group_a$a1
Id M mm x b u k j y p v
1 aAA 1 2 739 0.1257000 2 2 AA 2 8867 8926
2 a 1 2 2263 0.0004000 2 2 AR 4 7640 8926
3 aA 2 1 1 0.0845435 2 AA 2 6790 734,1092 NA
4 aAA 1 2 1965 0.0007000 4 3 AR 2 11616 8926
5 aAAA 1 3 3660 0.0008600 18 3 AA 2 20392 496
6 AA na 2 1972 0.0007000 11 3 AR 25 509 734
$group_a$a2
Id M mm x b u k j y p v
1 aAA 1 2 739 0.1257000 2 2 AA 2 8867 8926
2 a 1 2 2263 0.0004000 2 2 AR 4 7640 8926
3 aA 2 1 1 0.0845435 2 AA 2 6790 734,1092 NA
4 aAA 1 2 1965 0.0007000 4 3 AR 2 11616 8926
5 aAAA 1 3 3660 0.0008600 18 3 AA 2 20392 496
6 AA na 2 1972 0.0007000 11 3 AR 25 509 734
$group_a$a3
Id M mm x b u k j y p v
1 aAA 1 2 739 0.1257000 2 2 AA 2 8867 8926
2 a 1 2 2263 0.0004000 2 2 AR 4 7640 8926
3 aA 2 1 1 0.0845435 2 AA 2 6790 734,1092 NA
4 aAA 1 2 1965 0.0007000 4 3 AR 2 11616 8926
5 aAAA 1 3 3660 0.0008600 18 3 AA 2 20392 496
6 AA na 2 1972 0.0007000 11 3 AR 25 509 734
$group_b
$group_b$b1
Id M mm x b u k j y p v
1 aAA 1 2 739 0.1257000 2 2 AA 2 8867 8926
2 a 1 2 2263 0.0004000 2 2 AR 4 7640 8926
3 aA 2 1 1 0.0845435 2 AA 2 6790 734,1092 NA
4 aAA 1 2 1965 0.0007000 4 3 AR 2 11616 8926
5 aAAA 1 3 3660 0.0008600 18 3 AA 2 20392 496
6 AA na 2 1972 0.0007000 11 3 AR 25 509 734
$group_b$b2
Id M mm x b u k j y p v
1 aAA 1 2 739 0.1257000 2 2 AA 2 8867 8926
2 a 1 2 2263 0.0004000 2 2 AR 4 7640 8926
3 aA 2 1 1 0.0845435 2 AA 2 6790 734,1092 NA
4 aAA 1 2 1965 0.0007000 4 3 AR 2 11616 8926
5 aAAA 1 3 3660 0.0008600 18 3 AA 2 20392 496
6 AA na 2 1972 0.0007000 11 3 AR 25 509 734
$group_c
$group_c$c1
Id M mm x b u k j y p v
1 aAA 1 2 739 0.1257000 2 2 AA 2 8867 8926
2 a 1 2 2263 0.0004000 2 2 AR 4 7640 8926
3 aA 2 1 1 0.0845435 2 AA 2 6790 734,1092 NA
4 aAA 1 2 1965 0.0007000 4 3 AR 2 11616 8926
5 aAAA 1 3 3660 0.0008600 18 3 AA 2 20392 496
6 AA na 2 1972 0.0007000 11 3 AR 25 509 734
A.K.
From: Vera Costa veracosta...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Saturday, February 16, 2013 6:32 PM
Subject: Re: reading data
Sorry again... In:
res-do.call(c,lapply(list.files(recursive=T)[grep(...
What
Re: [R] list of matrices -- array
abind() (from package 'abind') can take a list of arrays as its first argument, so in general, no need for do.call() with abind(). As another poster pointed out, simplify2array() can also be used; while abind() gives more options regarding which dimension is created and how dimension names are constructed. x - list(A=cbind(X=c(a=1,b=2,c=3,d=4),Y=5:8,Z=9:12), B=cbind(X=c(a=13,b=14,c=15,d=16),Y=17:20,Z=21:24)) $A X Y Z a 1 5 9 b 2 6 10 c 3 7 11 d 4 8 12 $B X Y Z a 13 17 21 b 14 18 22 c 15 19 23 d 16 20 24 dim(abind(x, along=3)) [1] 4 3 2 dim(abind(x, along=1.5)) [1] 4 2 3 dim(abind(x, along=0.5)) [1] 2 4 3 dim(abind(x, along=1, hier.names=T)) # construct rownames in a hierarchical manner A.a, A.b, etc [1] 8 3 dim(abind(x, along=2, hier.names=T)) # construct colnames in a hierarchical manner [1] 4 6 abind(x, along=2, hier.names=T) A.X A.Y A.Z B.X B.Y B.Z a 1 5 9 13 17 21 b 2 6 10 14 18 22 c 3 7 11 15 19 23 d 4 8 12 16 20 24 On 2/14/2013 3:53 AM, Rolf Turner wrote: require(abind) do.call(abind,c(my_list,list(along=0))) # Gives 2 x 4 x 5 do.call(abind,c(my_list,list(along=3))) # Gives 4 x 5 x 2 The latter seems more natural to me. cheers, Rolf Turner On 02/14/2013 07:03 PM, Murat Tasan wrote: i'm somehow embarrassed to even ask this, but is there any built-in method for doing this: my_list - list() my_list[[1]] - matrix(1:20, ncol = 5) my_list[[2]] - matrix(20:1, ncol = 5) now, knowing that these matrices are identical in dimension, i'd like to unfold the list to a 2x4x5 (or some other permutation of the dim sizes) array. i know i can initialize the array, then loop through my_list to fill the array, but somehow this seems inelegant. i also know i can vectorize the matrices and unlist the list, then build the array from that single vector, but this also seems inelegant (and an easy place to introduce errors/bugs). i can't seem to find any built-in that handles this already... but maybe i just haven't looked hard enough :-/ __ 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] positioning a light source within a rgl-plot
Am 27.01.2013 16:48, schrieb Duncan Murdoch: On 13-01-27 10:37 AM, Alexander Senger wrote: Hello useRs, I would like to draw a 3D-surface using rgl with a point-like light-source within the scene, that is with finite distance of the light-source to the surface to be lit. The rgl package doesn't support that. From the help to the 'light3d' command I read: They [the light-sources] are positioned either in world space or relative to the camera using polar coordinates. which *could* be understood as if such a thing would be possible. But probably this is wishful thinking as my naive approach: light3d(x = 0, y = 0, z = 1) There are no x, y, z arguments to light3d. Only directional sources at infinite distance are supported. gives an error about un-used arguments in the function call. Also skimming the web does not produce any helpful examples. So please advise if there is a way to achieve my desired setting. Alternatively any hint how and where to make (moderate) modifications to the source code to get this functionality would be very welcome. It's rather tedious to make the changes. You need to change rgl.light, the rgl_light function in api.cpp that it calls, and the parts of light.hpp and light.cpp that are called by that. For completeness you'd also want to fix writeWebGL and scene3d and the functions they call. Duncan Murdoch Hi Duncan, thank you very much for your helpful advice. At the moment I have everything working except writeWebGL. I had a fair look into the matter and thought about implementing a slightly more sophisticated per-fragment shading and including all 8 light sources available in OpenGL. However I wonder that you probably had very good reason to limit the functionality to the simpler lighting model as is. Maybe you could share some thoughts on why making the lighting more complex might not be such a good idea, before I put too much effort into this topic. Thanks again, Alexander Senger __ 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] detecting entry into a recursive function
Thanks Jim -- I had considered this approach; is there any way to hide
such arguments from users?
Jim Lemon wrote:
On 02/17/2013 12:55 PM, Benjamin Tyner wrote:
Given a function that calls itself, what's the best way to detect the
entry point? The best I came up with is:
IsEntryPoint - function(){
par - sys.call(-1L)[[1]]
grandpar - sys.call(-2L)[[1]]
!identical(par, grandpar)
}
but this won't work for functions that don't directly call themselves;
for example, in this one the paste gets inserted in the call stack, so i
is always TRUE.
f - function(d){
i - IsEntryPoint()
if(d 1L) paste(d, f(d-1L))
}
Any ideas?
Hi Ben,
There are a number of ways to do this, from simply having a flag that
has one value at the initial call to the function and another value
when the function calls itself. for an example of this, see the code
for the sizetree function in the plotrix package, in particular the
firstcall argument.
Jim
__
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] Multidimensional correlation matrix question
Have you read An Introduction to R (ships with R) or made any other minimal effort with a beginning R tutorial to learn the language? Have you read the posting guide or any other previous posts to learn how to post coherently -- a small, reproducible example would be helpful, for instance. It certainly does not appear that you have done either. At present, your query appears unanswerable because we do not know the data structure of your data. Perhaps that is why you have received no replies? I suspect if you go through a tutorial, you will find that what you wish to do is simple. -- Cheers, Bert On Sat, Feb 16, 2013 at 11:50 PM, Tanushree Tunstall t...@gurukuli.co.uk wrote: , -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm __ 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] detecting entry into a recursive function
Make the flag an attribute of the function? Unless the user looks at
the attributes, it will be invisible.
(Not sure this does what you want, but maybe it's useful.)
-- Bert
On Sun, Feb 17, 2013 at 7:03 AM, Benjamin Tyner bty...@gmail.com wrote:
Thanks Jim -- I had considered this approach; is there any way to hide
such arguments from users?
Jim Lemon wrote:
On 02/17/2013 12:55 PM, Benjamin Tyner wrote:
Given a function that calls itself, what's the best way to detect the
entry point? The best I came up with is:
IsEntryPoint - function(){
par - sys.call(-1L)[[1]]
grandpar - sys.call(-2L)[[1]]
!identical(par, grandpar)
}
but this won't work for functions that don't directly call themselves;
for example, in this one the paste gets inserted in the call stack, so i
is always TRUE.
f - function(d){
i - IsEntryPoint()
if(d 1L) paste(d, f(d-1L))
}
Any ideas?
Hi Ben,
There are a number of ways to do this, from simply having a flag that has
one value at the initial call to the function and another value when the
function calls itself. for an example of this, see the code for the
sizetree function in the plotrix package, in particular the firstcall
argument.
Jim
__
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.
--
Bert Gunter
Genentech Nonclinical Biostatistics
Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm
__
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] Multidimensional correlation matrix question
https://github.com/hadley/devtools/wiki/Reproducibility John Kane Kingston ON Canada -Original Message- From: t...@gurukuli.co.uk Sent: Sun, 17 Feb 2013 07:50:17 + To: r-help@r-project.org Subject: [R] Multidimensional correlation matrix question Hello, I previously sent this message which was stripped off due to HTML mail. Sorry! Thanks. Hello All, I am new to the list. I have been learning to use R recently and its been great to see so much help available. However I must admit, I have stumbled upon a problem with correlation matrices and was hoping if someone could help. I have an excel workbook with a sheet per drug. Each sheet has 40 patients' drug response measured over different time points (several days). I have already been able to create a cor() matrix of all the drugs for all the patients on each day. However I have been asked to prepare a matrix of all days of all people for all drugs. So since there are 26 drugs, I am expected to form a 26x26 matrix for all this data. I am not even sure if this is possible. Can anyone point me in the right direction? PS: If I make individual matrix and combine them using cbind(), it is not a 26x26 sqaure matrix. Any help would be greatly appreciated. Thanks __ 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. FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks orcas on your desktop! __ 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] What value to put in range when I make kriging interpolation?
Hi, I am trying to make an automatic kriging interpolation algorithm.When I use the fit.variogram function what would it be a good startingvalue for the range? ThanksDimitris [[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] Hyperparameters in ARIMA models with dlm package
Hi, i'm beginner in Bayesian methods, I'm reading the documentation about dlm package and kalman filters, I'm looking for a example of transformation of ARIMA in a state space equivalent to use the dlm package and calcualte the hyperparameters. Someone can help me about it?. If it's possible with a arima(1,0,1) example, or more complex model. While I have more examples best for me. Thanks all [[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] How to do a backward calculation for each record in a dataset
Hi Experts, I have a dataset of 3 columns: customer.name product cost John Toothpaste 30 Mike Toothpaste 45 Peter Toothpaste 40 And I have a function of cost whereby cost = 3.40 + (1.20 * no.of.orders^2) I want to do a backward calculation for each records (each customer) to find his no.of.orders and create a new column named no.of.orders in that dataset but I don't know how to do. Please help me. Thank you everyone, Prakasit [[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] How to findout the name of a dataframe
Let'say we have a dataframe mydata with column v1. If mydata$v1 is passed to a function, is there way, then, to extract the name of the dataframe? What I now do is passing the name of the dataframe to the funcion, so passing two parameters. Maybe with mydata$v1 it is not possible, but with mydata['v1'] or mydata[,'v1'] it is? Thanks Frans --- Frans Marcelissen fransiepansiekever...@gmail.com [[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] Select components of a list
Sorry, I just noticed that in the first line of the solution 'l' got stripped
off, it should be read as 'lapply` instead of 'apply`.
lapply(1:length(models),function(i) lapply(models[[i]],function(x)
summary(x)$p.table[2,]))[[1]] #1st list component
Thanks,
A.K.
From: Gustav Sigtuna gsigt...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Sunday, February 17, 2013 5:49 AM
Subject: Re: Select components of a list
Dear Arun,
Thanks again. The script works perfectly for GLM. Strangely, it does not work
for GAM, although it has the same output for the linear part. I cannot figure
out the error message. I have attached the gam code and error message .
Thanks
On Sun, Feb 17, 2013 at 3:11 AM, arun smartpink...@yahoo.com wrote:
HI Gustav,
If you need the combined output:
res-lapply(1:length(models),function(i)
do.call(rbind,lapply(models[[i]],function(x)
summary(x)$coef[row.names(summary(x)$coef)%in%c(pm10,ozone,so2),c(1:2,4)])))
names(res)-1:length(res)
res1-do.call(rbind,lapply(res,function(i)
{row.names(i)-c(pm10,ozone,so2);data.frame(i)}))
names(res1)[2:3]- c(Std.Error,Pr(|z|))
res1
# Estimate Std.Error Pr(|z|)
#1.pm10 0.0005999185 0.0001486195 5.423004e-05
#1.ozone 0.0010117294 0.0003792739 7.640816e-03
#1.so2 0.0026595441 0.0009352046 4.457766e-03
#2.pm10 0.0005720549 0.0001740368 1.012696e-03
#2.ozone 0.0009128304 0.0004364390 3.647954e-02
#2.so2 0.0028256121 0.0010150314 5.373144e-03
#3.pm10 0.0005099552 0.0001559620 1.076462e-03
#3.ozone 0.0023896044 0.0004109854 6.087769e-09
#3.so2 0.0024097381 0.0009563814 1.174744e-02
#4.pm10 0.0009285593 0.0001766520 1.468764e-07
#4.ozone 0.0005455392 0.0004301502 2.047076e-01
#4.so2 0.0017251400 0.0011635156 1.381552e-01
A.K.
From: Gustav Sigtuna gsigt...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Saturday, February 16, 2013 7:44 PM
Subject: Re: Select components of a list
Hi Arun,
Thanks for taking your time to find a solution.
I have attached a R script that will recreate a comparable list from publicly
available data. My list is longer and created by various models than the one
created here. However, the final output is similar to the one produced by
script. My interest is to extract only the coefficients for pm10., ozone and
so2 ( Estimate, Std. Error and p value) .
Thanks
On Fri, Feb 15, 2013 at 9:04 PM, arun smartpink...@yahoo.com wrote:
Dear Gustav,
Thank you for the data. Could you select a smaller subset of the list and
dput() that subset? Your data is useful, but I would have to recreate list
of lists from that to test and sometimes that may not accurate represent the
format in your list as it is the summary().
Arun
From: Gustav Sigtuna gsigt...@gmail.com
To: smartpink...@yahoo.com
Sent: Friday, February 15, 2013 4:56 AM
Subject: Re: Select components of a list
Hi Arun,
Thanks for your help. Your mail landed in my spam folder and just saw it by
chance.
I have attached a text file that contains the list of my model. It was
extremely long, thus I took out the last part which is think is more
important.
In brief I have an output from GAM model which resulted from analysis of
ozone at three time points on 12 data sets
Thanks for your assistance
On Wed, Feb 13, 2013 at 8:21 PM, smartpink...@yahoo.com wrote:
Dear Lungo,
If you can email (smartpink...@yahoo.com) me the `list` (dput(list)), I can
take a look at it. Probably, you understand that my previous solution was
just guesswork. With regards to GLM, GAM, it is good to check the structure
of the list (str()). It gives information about whether a `generic` tool
could be applied to extract them or not.
Cheers.
Arun
quote author='Lungo'
Dear Arun, Your code and the example works fine. However my list is quite
different from the one showed in your example. As I have shown in my
question above I have 12 lists each having 3 lists underneath. I get the
lists by different models (GLM, GAM ) but the output I aim to have is the
estimates of the explanatory variable which is placed next to the intercept.
Thus I am looking for a “generic” tool that would extract these lists.
/quote
Quoted from:
http://r.789695.n4.nabble.com/Select-components-of-a-list-tp4658295p4658389.html
__
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] calculate classification error in test data set if training data set is not the same length
Dear R-Help Members, I have built a classification function using a baseline data set, that contains the group variable and have used it to classify the test data set. I am now trying to get the classification table for the training and test data set and classification success using: baseline.lda-lda(Stock ~ LTT + LF + LFM + LPO + LH + LPV + LPA + LD + LA + DAC + HH + HP + ML + OD + TV02 + TV03 + TV04 + TV05 + TH01 + TH06 + TH07 + TH08 + TD02 + TD03 + TD05 + TD06 + TD08 + WGHT + WDTH + PERI + CIRC + A02 + A03 + A04 + A05 + A08 + A12 + A14 + A15 + A16 + A17 + A18 + A26 + A27 + A30 + B02 + B03 + B08 + B09 + B10 + B11 + B14 + B15 + B20 + B22 + B24 + B25 + B26 + B29 + B30 + C02 + C03 + C05 + C06 + C08 + C09 + C11 + C12 + C13 + C14 + C16 + C18 + C19 + C23 + C25 + C28 + C29 + D01 + D02 + D03 + D04 + D06 + D07 + D08 + D10 + D14 + D15 + D18 + D19 + D20 + D23,data=baseline.data.scaled) ypred.train - predict(baseline.lda,baseline.data.scaled)$class ypred.test - predict(baseline.lda,mixed.data.scaled)$class # Training error table(ypred.train , baseline.data.scaled$Stock) mean(ypred.train == baseline.data.scaled$Stock) # Test error table(ypred.test , baseline.data.scaled$Stock) mean(ypred.test == baseline.data.scaled$Stock) This works for the training set, but not the test data set, as the baseline data set is the only one that contains the grouping variable, but is not the same length as the test data set (1161 samples in the training set, 236 in the test set). Is there a way to construct a classification table from the test data predictions and get the classification error? Thank you! Melanie Zoelck Melanie Zölck (Zoelck) PhD Candidate Galway-Mayo Institute of Technology Marine and Freshwater Research Centre Commercial Fisheries Research Group Department of Life Science Dublin Road Galway Republic of Ireland Mobile: +353 (0) 85 7246196 Skype: melaniezoelck E-mail: mzoe...@gmx.com or mzoe...@hotmaill.com __ 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] Select components of a list
Hi Gustav,
Just change `summary(x)$coef` to `summary(x)$p.table`
I am pasting the code from the attachment.
library(gamair)
library(mgcv)
data(chicago)
library(splines)
chicago$date-seq(from=as.Date(1987-01-01),
to=as.Date(2000-12-31),length=5114)
chicago$trend-seq(dim(chicago)[1])
names(chicago) [2] -pm10
names(chicago) [3] -pm25
names(chicago) [4] -ozone
names(chicago) [5] -so2
names(chicago) [7] -temp
chicago$trend-seq(dim(chicago)[1])
chicago$year-as.numeric(format(chicago$date,%Y))
chicago1-subset(chicago, as.Date(date) '1999-01-01')
year- matrix(1987:1998, ncol=3, byrow=TRUE)
fun -
function( y , x ){
a - gam(
death ~ pm10 + s(trend,k=35) , poisson , na.action = na.omit , data = x[
x$year %in% y , ]
)
b- gam(
death ~ ozone + s(trend,k=35), poisson , na.action = na.omit , data = x[
x$year %in% y , ]
)
c- gam (
death ~ so2 + ns(trend,k=35) , poisson , na.action = na.omit , data = x[
x$year %in% y , ]
)
list( a , b ,c)
}
models- apply(year, 1 , fun , x = chicago )
#solution
apply(1:length(models),function(i) lapply(models[[i]],function(x)
summary(x)$p.table[2,]))[[1]] #1st list component
#[[1]]
# Estimate Std. Error z value Pr(|z|)
#6.054413e-04 1.474943e-04 4.104845e+00 4.045864e-05
#
#[[2]]
# Estimate Std. Error z value Pr(|z|)
#0.000765 0.0003777224 2.6473846529 0.0081117027
#[[3]]
# Estimate Std. Error z value Pr(|z|)
#5.643234e-03 9.023766e-04 6.253746e+00 4.007234e-10
res-lapply(1:length(models),function(i)
do.call(rbind,lapply(models[[i]],function(x)
summary(x)$p.table[row.names(summary(x)$p.table)%in%c(pm10,ozone,so2),c(1:2,4)])))
names(res)-1:length(res)
res1- lapply(res,function(i)
{row.names(i)-c(pm10,ozone,so2);data.frame(i)})
library(abind)
res2-abind(res1,along=1,hier.names=T) #gives a matrix
colnames(res2)[2:3]- c(Std.Error,Pr(|z|))
res3- do.call(rbind,lapply(res,function(i)
{row.names(i)-c(pm10,ozone,so2);data.frame(i)}))
colnames(res3)[2:3]- c(Std.Error,Pr(|z|))
str(res2)
#num [1:12, 1:3] 0.000605 0.001 0.005643 0.00059 0.000839 ...
#- attr(*, dimnames)=List of 2
# ..$ : chr [1:12] 1.pm10 1.ozone 1.so2 2.pm10 ...
# ..$ : chr [1:3] Estimate Std.Error Pr(|z|)
str(res3)
#'data.frame': 12 obs. of 3 variables:
# $ Estimate : num 0.000605 0.001 0.005643 0.00059 0.000839 ...
# $ Std.Error: num 0.000147 0.000378 0.000902 0.000172 0.000427 ...
#$ Pr(|z|) : num 4.05e-05 8.11e-03 4.01e-10 6.23e-04 4.96e-02 ...
res2
# Estimate Std.Error Pr(|z|)
#1.pm10 0.0006054413 0.0001474943 4.045864e-05
#1.ozone 0.000765 0.0003777224 8.111703e-03
#1.so2 0.0056432338 0.0009023766 4.007234e-10
#2.pm10 0.0005899052 0.0001724085 6.226404e-04
#2.ozone 0.0008389801 0.0004272480 4.956676e-02
#2.so2 0.0032899751 0.0009475318 5.163027e-04
#3.pm10 0.0005398889 0.0001551911 5.035438e-04
#3.ozone 0.0023890220 0.0004082119 4.845107e-09
#3.so2 0.0049121476 0.0008818088 2.539574e-08
#4.pm10 0.0009341888 0.0001760271 1.113999e-07
#4.ozone 0.0005461742 0.0004253987 1.991731e-01
#4.so2 0.0055712219 0.0011123419 5.484117e-07
Hope it helps.
A.K.
From: Gustav Sigtuna gsigt...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Sunday, February 17, 2013 5:49 AM
Subject: Re: Select components of a list
Dear Arun,
Thanks again. The script works perfectly for GLM. Strangely, it does not work
for GAM, although it has the same output for the linear part. I cannot figure
out the error message. I have attached the gam code and error message .
Thanks
On Sun, Feb 17, 2013 at 3:11 AM, arun smartpink...@yahoo.com wrote:
HI Gustav,
If you need the combined output:
res-lapply(1:length(models),function(i)
do.call(rbind,lapply(models[[i]],function(x)
summary(x)$coef[row.names(summary(x)$coef)%in%c(pm10,ozone,so2),c(1:2,4)])))
names(res)-1:length(res)
res1-do.call(rbind,lapply(res,function(i)
{row.names(i)-c(pm10,ozone,so2);data.frame(i)}))
names(res1)[2:3]- c(Std.Error,Pr(|z|))
res1
# Estimate Std.Error Pr(|z|)
#1.pm10 0.0005999185 0.0001486195 5.423004e-05
#1.ozone 0.0010117294 0.0003792739 7.640816e-03
#1.so2 0.0026595441 0.0009352046 4.457766e-03
#2.pm10 0.0005720549 0.0001740368 1.012696e-03
#2.ozone 0.0009128304 0.0004364390 3.647954e-02
#2.so2 0.0028256121 0.0010150314 5.373144e-03
#3.pm10 0.0005099552 0.0001559620 1.076462e-03
#3.ozone 0.0023896044 0.0004109854 6.087769e-09
#3.so2 0.0024097381 0.0009563814 1.174744e-02
#4.pm10 0.0009285593 0.0001766520 1.468764e-07
#4.ozone 0.0005455392 0.0004301502 2.047076e-01
#4.so2 0.0017251400 0.0011635156 1.381552e-01
A.K.
From: Gustav Sigtuna gsigt...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Saturday, February 16, 2013 7:44 PM
Subject: Re: Select components of a list
Hi Arun,
Thanks for taking your time to find a solution.
I have attached a R script
[R] nested random factor using lme produces errors
Hi, I am running a mixed-effect model with a nested-random effect. I am interested in gut parasites in moose. I has three different type of treatment that I applied to moose which are from different families. My response variable is gut parasites and the factors are moose families which is nested within treatment. My data is balanced. To answer this question, I used the lme function like this : model=lme(parasite~drug,random=~1|drug/family) But doing a summary on this model gives me warning message : In pt(-abs(tTable[, t-value]), tTable[, DF]) : NaNs produced I don't understand why ?! I noticed that the p-values are not computed and have NAs values for drug2 and drug3 (from the summary of this model) Moreover, in the summary, I noticed that in the random effects line I have standard deviation for Formula: ~1 | drug and for Formula: ~1 | family %in% drug. Does R consider drug as a random factor as well ? And last question, how can I know if my random factor has a significant effect on the gut parasites ? Thank for your help, -- View this message in context: http://r.789695.n4.nabble.com/nested-random-factor-using-lme-produces-errors-tp4658837.html Sent from the R help mailing list archive at Nabble.com. __ 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 do a backward calculation for each record in a dataset
This sounds a bit too much like homework. And in any case https://github.com/hadley/devtools/wiki/Reproducibility John Kane Kingston ON Canada -Original Message- From: asltjoey.rs...@gmail.com Sent: Sun, 17 Feb 2013 20:10:13 +0700 To: r-help@r-project.org Subject: [R] How to do a backward calculation for each record in a dataset Hi Experts, I have a dataset of 3 columns: customer.name product cost John Toothpaste 30 Mike Toothpaste 45 Peter Toothpaste 40 And I have a function of cost whereby cost = 3.40 + (1.20 * no.of.orders^2) I want to do a backward calculation for each records (each customer) to find his no.of.orders and create a new column named no.of.orders in that dataset but I don't know how to do. Please help me. Thank you everyone, Prakasit [[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. FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks orcas on your desktop! __ 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 do a backward calculation for each record in a dataset
Hi, I am not sure I understand it correctly. dat1-read.table(text= customer.name product cost John Toothpaste 30 Mike Toothpaste 45 Peter Toothpaste 40 ,sep=,header=TRUE,stringsAsFactors=FALSE) dat1$no.of.orders- sqrt((dat1$cost-3.40)/1.20) dat1 # customer.name product cost no.of.orders #1 John Toothpaste 30 4.708149 #2 Mike Toothpaste 45 5.887841 #3 Peter Toothpaste 40 5.522681 A.K. - Original Message - From: Prakasit Singkateera asltjoey.rs...@gmail.com To: r-help@r-project.org Cc: Sent: Sunday, February 17, 2013 8:10 AM Subject: [R] How to do a backward calculation for each record in a dataset Hi Experts, I have a dataset of 3 columns: customer.name product cost John Toothpaste 30 Mike Toothpaste 45 Peter Toothpaste 40 And I have a function of cost whereby cost = 3.40 + (1.20 * no.of.orders^2) I want to do a backward calculation for each records (each customer) to find his no.of.orders and create a new column named no.of.orders in that dataset but I don't know how to do. Please help me. Thank you everyone, Prakasit [[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.
Re: [R] nested random factor using lme produces errors
melswed amelie.truchy at slu.se writes: Hi, I am running a mixed-effect model with a nested-random effect. I am interested in gut parasites in moose. I has three different type of treatment that I applied to moose which are from different families. My response variable is gut parasites and the factors are moose families which is nested within treatment. My data is balanced. To answer this question, I used the lme function like this : model=lme(parasite~drug,random=~1|drug/family) But doing a summary on this model gives me warning message : In pt(-abs(tTable[, t-value]), tTable[, DF]) : NaNs produced I don't understand why ?! I noticed that the p-values are not computed and have NAs values for drug2 and drug3 (from the summary of this model) Moreover, in the summary, I noticed that in the random effects line I have standard deviation for Formula: ~1 | drug and for Formula: ~1 | family %in% drug. Does R consider drug as a random factor as well ? And last question, how can I know if my random factor has a significant effect on the gut parasites ? This belongs on r-sig-mixed-mod...@r-project.org. Hint: it very rarely makes sense to include a categorical predictor such as drug as both a random and a fixed effect ... this model is overspecified. For computational and philosophical reasons, it seems unwise and odd (respectively) to treat drug as a random effect. I might have more to say but will say it (perhaps) if you repost on r-sig-mixed-models . __ 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 findout the name of a dataframe
Will this work for you:
myFunc - function(var){
+ # get the dataframe name
+ charName - deparse(substitute(var))
+ # parse out data.frame
+ dataFrame - sub(\\$.*, , charName)
+ cat(input:, charName, data.frame:, dataFrame, \n)
+ }
myFunc(mydata$V1)
input: mydata$V1 data.frame: mydata
On Sun, Feb 17, 2013 at 8:51 AM, Frans Marcelissen
frans.marcelis...@digipsy.nl wrote:
Let'say we have a dataframe mydata with column v1. If mydata$v1 is passed
to a function, is there way, then, to extract the name of the dataframe?
What I now do is passing the name of the dataframe to the funcion, so
passing two parameters. Maybe with mydata$v1 it is not possible, but with
mydata['v1'] or mydata[,'v1'] it is?
Thanks
Frans
---
Frans Marcelissen
fransiepansiekever...@gmail.com
[[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.
--
Jim Holtman
Data Munger Guru
What is the problem that you are trying to solve?
Tell me what you want to do, not how you want to do it.
__
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 do a backward calculation for each record in a dataset
Homework? We don't do homework here. -- Bert On Sun, Feb 17, 2013 at 5:10 AM, Prakasit Singkateera asltjoey.rs...@gmail.com wrote: Hi Experts, I have a dataset of 3 columns: customer.name product cost John Toothpaste 30 Mike Toothpaste 45 Peter Toothpaste 40 And I have a function of cost whereby cost = 3.40 + (1.20 * no.of.orders^2) I want to do a backward calculation for each records (each customer) to find his no.of.orders and create a new column named no.of.orders in that dataset but I don't know how to do. Please help me. Thank you everyone, Prakasit [[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. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm __ 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 findout the name of a dataframe
On Feb 17, 2013, at 5:51 AM, Frans Marcelissen wrote: Let'say we have a dataframe mydata with column v1. If mydata$v1 is passed to a function, is there way, then, to extract the name of the dataframe? What I now do is passing the name of the dataframe to the funcion, so passing two parameters. Maybe with mydata$v1 it is not possible, but with mydata['v1'] or mydata[,'v1'] it is? It will depend on the specifics. The usual way is with deparse(substitute(arg)) d - data.frame(a=a) gn - function(col) print(deparse(substitute(col))) gn(d) [1] d gn(d$a) [1] d$a You do realize that mydata$v1 is identical (after evaluation, anyway) to mydata[,'v1'] , but not to mydata['v1'], don't you? gn(d['a']) [1] d[\a\] Thanks Frans --- Frans Marcelissen fransiepansiekever...@gmail.com [[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. David Winsemius Alameda, CA, USA __ 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] list of matrices -- array
thanks to all! didn't know about simplify2array, nor about the abind package. they're exactly what i wanted. cheers, -m On Sun, Feb 17, 2013 at 9:41 AM, Tony Plate tapl...@gmail.com wrote: abind() (from package 'abind') can take a list of arrays as its first argument, so in general, no need for do.call() with abind(). As another poster pointed out, simplify2array() can also be used; while abind() gives more options regarding which dimension is created and how dimension names are constructed. x - list(A=cbind(X=c(a=1,b=2,c=3,d=4),Y=5:8,Z=9:12), B=cbind(X=c(a=13,b=14,c=15,d=16),Y=17:20,Z=21:24)) $A X Y Z a 1 5 9 b 2 6 10 c 3 7 11 d 4 8 12 $B X Y Z a 13 17 21 b 14 18 22 c 15 19 23 d 16 20 24 dim(abind(x, along=3)) [1] 4 3 2 dim(abind(x, along=1.5)) [1] 4 2 3 dim(abind(x, along=0.5)) [1] 2 4 3 dim(abind(x, along=1, hier.names=T)) # construct rownames in a hierarchical manner A.a, A.b, etc [1] 8 3 dim(abind(x, along=2, hier.names=T)) # construct colnames in a hierarchical manner [1] 4 6 abind(x, along=2, hier.names=T) A.X A.Y A.Z B.X B.Y B.Z a 1 5 9 13 17 21 b 2 6 10 14 18 22 c 3 7 11 15 19 23 d 4 8 12 16 20 24 On 2/14/2013 3:53 AM, Rolf Turner wrote: require(abind) do.call(abind,c(my_list,list(along=0))) # Gives 2 x 4 x 5 do.call(abind,c(my_list,list(along=3))) # Gives 4 x 5 x 2 The latter seems more natural to me. cheers, Rolf Turner On 02/14/2013 07:03 PM, Murat Tasan wrote: i'm somehow embarrassed to even ask this, but is there any built-in method for doing this: my_list - list() my_list[[1]] - matrix(1:20, ncol = 5) my_list[[2]] - matrix(20:1, ncol = 5) now, knowing that these matrices are identical in dimension, i'd like to unfold the list to a 2x4x5 (or some other permutation of the dim sizes) array. i know i can initialize the array, then loop through my_list to fill the array, but somehow this seems inelegant. i also know i can vectorize the matrices and unlist the list, then build the array from that single vector, but this also seems inelegant (and an easy place to introduce errors/bugs). i can't seem to find any built-in that handles this already... but maybe i just haven't looked hard enough :-/ __ 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] Hyperparameters in ARIMA models with dlm package
Hi: Like I said earlier, you really should read west and harrison first, especially if you're a beginner in bayesian methods. giovanni's book and package are both very nice ( thanks giovanni ) but the book is more of a summary of west and harrison and sort of assumes some familarity with the material. I checked the dlm book and section 3.2.5 gives an example of conversion of arma(1,1) to a DLM. there is also a function in the dlm package dlmModArma that does the general conversion to state space. But, when converting from arima to state space, one has to consider constraints on the arima parameters to in order maintain stationarity in the original arima space. I'm not sure if dlmModArma handles that. But it's a start. Maybe Giovanni will see your question and say more. On Sun, Feb 17, 2013 at 10:49 AM, Juan Manuel Becerra jmblue...@gmail.comwrote: Hi, i'm beginner in Bayesian methods, I'm reading the documentation about dlm package and kalman filters, I'm looking for a example of transformation of ARIMA in a state space equivalent to use the dlm package and calcualte the hyperparameters. Someone can help me about it?. If it's possible with a arima(1,0,1) example, or more complex model. While I have more examples best for me. Thanks all [[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. [[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] detecting entry into a recursive function
On 13-02-17 10:03 AM, Benjamin Tyner wrote:
Thanks Jim -- I had considered this approach; is there any way to hide
such arguments from users?
Don't export the recursive function, just a non-recursive function that
calls it.
Duncan Murdoch
Jim Lemon wrote:
On 02/17/2013 12:55 PM, Benjamin Tyner wrote:
Given a function that calls itself, what's the best way to detect the
entry point? The best I came up with is:
IsEntryPoint - function(){
par - sys.call(-1L)[[1]]
grandpar - sys.call(-2L)[[1]]
!identical(par, grandpar)
}
but this won't work for functions that don't directly call themselves;
for example, in this one the paste gets inserted in the call stack, so i
is always TRUE.
f - function(d){
i - IsEntryPoint()
if(d 1L) paste(d, f(d-1L))
}
Any ideas?
Hi Ben,
There are a number of ways to do this, from simply having a flag that
has one value at the initial call to the function and another value
when the function calls itself. for an example of this, see the code
for the sizetree function in the plotrix package, in particular the
firstcall argument.
Jim
__
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] histogram
HI Elisa, You could use ?cut() vec1-c(33,18,13,47,30,10,6,21,39,25,40,29,14,16,44,1,41,4,15,20,46,32,38,5,31,12,48,27,36,24,34,2,35,11,42,9,8,7,26,22,43,17,19,28,23,3,49,37,50,45) label1-unlist(lapply(mapply(c,lapply(seq(0,45,5),function(x) x),lapply(seq(5,50,5),function(x) x),SIMPLIFY=FALSE),function(i) paste(i[1],x=,i[2],sep=))) Count1-as.data.frame(table(cut(vec1,breaks=seq(0,50,5),labels=label1))) Count1 # Var1 Freq #1 0x=5 5 #2 5x=10 5 #3 10x=15 5 #4 15x=20 5 #5 20x=25 5 #6 25x=30 5 #7 30x=35 5 #8 35x=40 5 #9 40x=45 5 #10 45x=50 5 hist(vec1,breaks=seq(0,50,5),freq=T) hist(vec1,breaks=seq(0,50,5),prob=TRUE) lines(density(vec1)) label2-unlist(lapply(mapply(c,lapply(seq(0,45,5),function(x) x),lapply(seq(55,by=50,length.out=10),function(x) x),SIMPLIFY=FALSE),function(i) paste(i[1],x=,i[2],sep=))) count2-as.data.frame(table(cut(vec1,breaks=c(0,55,100,145,190,235,280,325,370,415,460),labels=label2))) count2 # Var1 Freq #1 0x=55 50 #2 5x=105 0 #3 10x=155 0 #4 15x=205 0 #5 20x=255 0 #6 25x=305 0 #7 30x=355 0 #8 35x=405 0 #9 40x=455 0 #10 45x=505 0 hist(vec1,breaks=c(0,55,100,145,190,235,280,325,370,415,460)) hist(vec1,breaks=c(0,55,100,145,190,235,280,325,370,415,460),prob=TRUE) lines(density(vec1)) A.K. From: eliza botto eliza_bo...@hotmail.com To: smartpink...@yahoo.com smartpink...@yahoo.com Sent: Sunday, February 17, 2013 2:35 PM Subject: histogram Dear Arun, [text file is attached in case format of email is changed] For the following data set 33 18 13 47 30 10 6 21 39 25 40 29 14 16 44 1 41 4 15 20 46 32 38 5 31 12 48 27 36 24 34 2 35 11 42 9 8 7 26 22 43 17 19 28 23 3 49 37 50 45 1. i first of all want to make classes in the following form class 0x=5 5x=10 10x=15 15x=20 . ... ... ... 45x=50 and then i want to count the number of elements in each class and ultimately i want to execute a table in the following form classnumber of elements in each class 0x=55 5x=105 10x=155 15x=205 . ... ... ... 45x=505 the command which i used is to count the number of elements in each class was length(which(x 45 x = 50)) 2. is there a loop command which can count the number of elements in each range instead of using an individual command for each range?? 3. then i want to make histogram out of that table. 4. lastly, i want to fit density curve on those histograms. I am greatful to you in advance. elisa __ 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] addition in the initial question
Dear Elisa, Try this: vec1-c(33,18,13,47,30,10,6,21,39,25,40,29,14,16,44,1,41,4,15,20,46,32,38,5,31,12,48,27,36,24,34,2,35,11,42,9,8,7,26,22,43,17,19,28,23,3,49,37,50,45) vec2-vec1[1:26] names(vec2)-LETTERS[1:26] label1-unlist(lapply(mapply(c,lapply(seq(0,45,5),function(x) x),lapply(seq(5,50,5),function(x) x),SIMPLIFY=FALSE),function(i) paste(i[1],x=,i[2],sep=))) dat1-data.frame(vec2,class=cut(vec2,breaks=seq(0,50,5),labels=label1),stringsAsFactors=FALSE) vecNew-cut(vec2,breaks=seq(0,50,5),labels=label1) res-data.frame(aggregate(row.names(dat1)~class,dat1,function(x) x),Frequency=as.data.frame(table(vecNew))[,2],stringsAsFactors=FALSE) names(res)[2]- header_elements_class res # class header_elements_class Frequency #1 0x=5 P, R, X 3 #2 5x=10 F, G 4 #3 10x=15 C, M, S, Z 3 #4 15x=20 B, N, T 2 #5 20x=25 H, J 2 #6 25x=30 E, L 3 #7 30x=35 A, V, Y 3 #8 35x=40 I, K, W 2 #9 40x=45 O, Q 2 #10 45x=50 D, U 2 str(res) #'data.frame': 10 obs. of 3 variables: # $ class : Factor w/ 10 levels 0x=5,5x=10,..: 1 2 3 4 5 6 7 8 9 10 # $ header_elements_class:List of 10 # ..$ 0: chr P R X # ..$ 1: chr F G # ..$ 2: chr C M S Z # ..$ 3: chr B N T # ..$ 4: chr H J # ..$ 5: chr E L # ..$ 6: chr A V Y # ..$ 7: chr I K W # ..$ 8: chr O Q # ..$ 9: chr D U # $ Frequency : int 3 4 3 2 2 3 3 2 2 2 A.K. From: eliza botto eliza_bo...@hotmail.com To: smartpink...@yahoo.com smartpink...@yahoo.com Sent: Sunday, February 17, 2013 3:27 PM Subject: addition in the initial question Dear Arun, just a small change in the initial question. what if, instead of counting number i want to get the headers of the numbers falling in each range. finally counting those numbers to get frequencies and then drawing histograms and fitting density curve. [HEADER] A B C D E F G H I J K L M N O P Q R S T U V W X Y Z [VECTOR] 33 18 13 47 30 10 6 21 39 25 40 29 14 16 44 1 41 4 15 20 46 32 38 5 31 12 classheader of elements in each class Frequency 0x=5P, R, X 3 5x=10F, G 2 10x=15 C, S 2 . ... ... ... 45x=50 U 1 thankyou very much in advance elisa __ 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] Multidimensional correlation matrix question
Thanks! I will do the needful. On Sun, Feb 17, 2013 at 3:32 PM, John Kane jrkrid...@inbox.com wrote: https://github.com/hadley/devtools/wiki/Reproducibility John Kane Kingston ON Canada -Original Message- From: t...@gurukuli.co.uk Sent: Sun, 17 Feb 2013 07:50:17 + To: r-help@r-project.org Subject: [R] Multidimensional correlation matrix question Hello, I previously sent this message which was stripped off due to HTML mail. Sorry! Thanks. Hello All, I am new to the list. I have been learning to use R recently and its been great to see so much help available. However I must admit, I have stumbled upon a problem with correlation matrices and was hoping if someone could help. I have an excel workbook with a sheet per drug. Each sheet has 40 patients' drug response measured over different time points (several days). I have already been able to create a cor() matrix of all the drugs for all the patients on each day. However I have been asked to prepare a matrix of all days of all people for all drugs. So since there are 26 drugs, I am expected to form a 26x26 matrix for all this data. I am not even sure if this is possible. Can anyone point me in the right direction? PS: If I make individual matrix and combine them using cbind(), it is not a 26x26 sqaure matrix. Any help would be greatly appreciated. Thanks __ 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. FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks orcas on your desktop! Check it out at http://www.inbox.com/marineaquarium [[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] Loop
Hi all,
I want to execute a loop of a program:
for (u in Timeframemin:Timeframe){}
Imagine that Timeframemin-10
Timefram-1
Is it posible to execute the loop but only proving from 10 to 1 but
jumping 10 each time, for example, execute for 10,20,30.to Timeframe.
Other question is, when a program is heavy and has a lot of loops to
execute (how can I know where loop is executing at the moment? There is
some print or something to see in wich step of the loop is the program to
see if the program is advancing?
Many thanks in advance.
[[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] strptime() with format %OS does not print millisecs in MacOS
Hi! I'm finding this on MacOS Lion 10.7.5 getOption(digits.secs) NULL a - 2012_10_01_14_13_32.445 strptime(a,format=%Y_%M_%d_%H_%M_%OS) [1] 2012-02-01 14:13:32 strptime(a,format=%Y_%M_%d_%H_%M_%OS3) [1] NA I can solve it with options(digits.secs=3) strptime(a,format=%Y_%M_%d_%H_%M_%OS) [1] 2012-02-01 14:13:32.445 But is this not in contradiction with the help page? if %OS is not followed by a digit, it uses the setting of getOption(digits.secs), or if that is unset, n = 3 Thanks Agus PS sessionInfo() R version 2.15.1 (2012-06-22) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] sp_0.9-99 loaded via a namespace (and not attached): [1] grid_2.15.1lattice_0.20-6 tools_2.15.1 -- Dr. Agustin Lobo Institut de Ciencies de la Terra Jaume Almera (CSIC) Lluis Sole Sabaris s/n 08028 Barcelona Spain Tel. 34 934095410 Fax. 34 934110012 e-mail agustin.l...@ictja.csic.es https://sites.google.com/site/aloboaleu/ __ 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] nested random factor using lme produces errors
I understand. I want to specify that drug is only a fixed factor and family should be the only random factor. So maybe, my R code is wrong If I specify random=~1|drug/family it is only because I wanted to specify that family is nested within drug. -- View this message in context: http://r.789695.n4.nabble.com/nested-random-factor-using-lme-produces-errors-tp4658837p4658878.html Sent from the R help mailing list archive at Nabble.com. __ 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] Loop
Hello,
As for the first question, you can do something like
Timeframemin-10
Timefram-1
uFrame - seq(Timeframemin, Timefram, by = 10)
for(u in uFrame) {}
As for the second question, the answer is yes, there is a print()
function, which can be used for your purpose.
Hope this helps,
Rui Barradas
Em 17-02-2013 18:53, Trying To learn again escreveu:
Hi all,
I want to execute a loop of a program:
for (u in Timeframemin:Timeframe){}
Imagine that Timeframemin-10
Timefram-1
Is it posible to execute the loop but only proving from 10 to 1 but
jumping 10 each time, for example, execute for 10,20,30.to Timeframe.
Other question is, when a program is heavy and has a lot of loops to
execute (how can I know where loop is executing at the moment? There is
some print or something to see in wich step of the loop is the program to
see if the program is advancing?
Many thanks in advance.
[[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.
Re: [R] detecting entry into a recursive function
On 02/18/2013 02:03 AM, Benjamin Tyner wrote:
Thanks Jim -- I had considered this approach; is there any way to hide
such arguments from users?
Hi Ben,
I played around with your solution for a while and if the first argument
to the function changes with each recursive call, you may be able to do
it like this:
factorial-function(x) {
cat(x,\n)
if(as.character(x) == as.character(sys.call())[2])
cat(I\'m the first!\n)
if(x1) x-x*factorial(x-1)
return(x)
}
factorial(3)
factorial(4)
...
Jim
__
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] questions hash functions
Hello R, could you explain to me how to resolve this question: If this is a matrix: Element S1 S2 S3 S4 001 01 101 00 210 01 300 10 400 11 510 00 1. How is possible to ompute the minhash signature for each column if we use the following three hash functions: h1(x) = 2x + 1 mod 6; h2(x) = 3x + 2 mod 6; h3(x) = 5x + 2 mod 6. 2. Which of these hash functions are true permutations? 3.How close are the estimated Jaccard similarities for the six pairs of columns to the true Jaccard similarities? Thank you! Tania __ 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] questions hash functions
I cannot imagine why you think this is the appropriate place to pose such a question. It certainly has nothing to do with R (which would be appropriate), and it does look like homework (which is identified as NOT appropriate in the Posting Guide for this mailing list). --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- Sent from my phone. Please excuse my brevity. Tania Patiño taniu...@gmail.com wrote: Hello R, could you explain to me how to resolve this question: If this is a matrix: Element S1 S2 S3 S4 001 01 101 00 210 01 300 10 400 11 510 00 1. How is possible to ompute the minhash signature for each column if we use the following three hash functions: h1(x) = 2x + 1 mod 6; h2(x) = 3x + 2 mod 6; h3(x) = 5x + 2 mod 6. 2. Which of these hash functions are true permutations? 3.How close are the estimated Jaccard similarities for the six pairs of columns to the true Jaccard similarities? Thank you! Tania __ 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] Computing Spectral Slope
Greetings,
I'm working on image classification and for that I want to use the spectral
slope as a feature for my classifier. For this I would prefer to calculate
this feature using R, so far I've read my image and converted it's RGB
representation into HSL. The spectral slope is computed over the Luminance
component, so at the moment what I have is a NxN matrix of Luminance values.
The following steps would be:
1- Calculating the FFT of the Luminance matrix.
2- Calculating the Power Spectra, which is computed as |FFT|^2.
3- Matching every Power Spectra value to its Spatial frequency and doing a
linear regression to compute the Spectral Slope.
My doubts start now: I know that there is a function called fft that
computes the fast fourier transform in R, but how does it exactly do it?
The formula for the DFT, for a NxN matrix is is F[u,v] = 1/N sum from {x=
0} to {N-1} sum from {y = 0} to {N-1} f(x,y) e^((-2*pi*j/N)*(xu+yv)),
however which coordinates origin does the R fft function assumes? Is it
possible to compute a centered fft, meaning that the origin of our
coordinate system is in the center of our matrix?
Then, is there any R function to calculate the power spectra of a matrix of
data without having to compute the FFT first?
Thank you in Advance :)
Marc
[[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] forecast ARMA(1,1)/GARCH(1,1) using fGarch library
Hi, i am working in the forecast of the daily price crude . The last prices of this data are the following: 100.60 101.47 100.20 100.06 98.68 101.28 101.05 102.13 101.70 98.27 101.00 100.50 100.03 102.23 102.68 103.32 102.67 102.23 102.14 101.25 101.11 99.90 98.53 96.76 96.12 96.54 96.30 95.92 95.92 93.45 93.71 96.42 93.99 93.76 95.24 95.63 95.95 95.83 95.65 96.61 91.30 91.66 96.23 94.44 94.50 96.52 97.07 97.37 95.31 96.10 94.35 93.34 93.68 93.65 95.16 94.32 94.82 94.93 95.72 96.41 96.70 95.87 95.46 96.83 96.49 96.70 99.61 100.84 99.90 99.65 99.22 98.84 99.08 97.53 98.51 99.17 100.07 101.49 102.40 103.24 102.36 100.70 100.93 104.43 105.67 106.23 109.98 108.80 109.10 108.86 108.68 109.59 110.41 The data consist of 2973 observations. For the analisys i considered the returns, the last ones are: 0.0066998270 0.0090753250 0.0141900670 0.0089664010 0.0082031250 -0.0085238280 -0.0162172720 0.0022840120 0.0346774990 0.0118739830 0.0052995170 0.0353007620 -0.0107292230 0.0027573530 -0.0021998170 -0.0016535000 0.0083732060 0.0074824350 For modelling the mean i fit an ARMA(1,1) and fot the volatility i fit a GARCH(1,1) , i used a t-student as conditional distribution, for this i used the fGarch librray, the code is the following: h-garchFit(~arma(1,1)+garch(2,2),data=R,cond.dist=std,TRACE=F) On the other hand, for the prediction i use the function predict. predict(h,10) meanForecast meanError standardDeviation 1 0.001451401 0.015316820.01531682 2 0.001265062 0.015400830.01539350 3 0.001263344 0.015496280.01548892 4 0.001263328 0.015573060.01556565 5 0.001263328 0.015664200.01565676 6 0.001263328 0.015740620.01573312 7 0.001263328 0.015828000.01582047 8 0.001263328 0.015903720.01589614 9 0.001263328 0.015987790.01598018 10 0.001263328 0.016062580.01605493 I am modelling this Y_t-mean=e_t=sigma_t*Z_t however, my question is ,the prediction for the return itself is the mean forecast? if this is the case my prediction for the price would be equal to (1+.001451401)*110.41 =110.57 but i think this is not a good prediction, because the volatility is not affecting so much , in addition the predicted prices are growing up nevertheless i would expect that at some point these ones decrease . So i would expect that the prediction for the return would be different. but i certanly dont know which is. I would apreciate if you could help with this. Greeting Gustavo [[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] cumulative sum by group and under some criteria
Hi,
Yes, I wanted to expand directly from d. If there are other variables,
says A, B, C that I need to keep in the final expanded data, how
to modify the code?
d-data.frame()
for (m1 in 2:3) {
for (n1 in 2:2) {
for (x1 in 0:(m1-1)) {
for (y1 in 0:(n1-1)) {
d-rbind(d,c(m1,n1,x1,y1))
}
}
}}
colnames(d)-c(m1,n1,x1,y1)
d
res1-do.call(rbind,lapply(unique(d$m1),function(m1)
do.call(rbind,lapply(unique(d$n1),function(n1)
do.call(rbind,lapply(0:(m1-1),function(x1)
do.call(rbind,lapply(0:(n1-1),function(y1)
do.call(rbind,lapply((m1+2):(7-n1),function(m)
do.call(rbind,lapply((n1+2):(9-m),function(n)
do.call(rbind,lapply(x1:(x1+m-m1), function(x)
expand.grid(m1,n1,x1,y1,m,n,x)) )
names(res1)- c(m1,n1,x1,y1,m,n,x)
set.seed(235)
d$A- sample(1:50,10,replace=TRUE)
set.seed(23)
d$B- sample(1:50,10,replace=TRUE)
d
# m1 n1 x1 y1 A B
#1 2 2 0 0 50 29
#2 2 2 0 1 40 12
#3 2 2 1 0 31 17
#4 2 2 1 1 7 36
#5 3 2 0 0 13 41
#6 3 2 0 1 27 22
#7 3 2 1 0 49 49
#8 3 2 1 1 47 49
#9 3 2 2 0 23 43
#10 3 2 2 1 4 50
library(plyr)
res2- join(res1,d,by=c(m1,n1,x1,y1),type=full)
res2
# m1 n1 x1 y1 m n x A B
#1 2 2 0 0 4 4 0 50 29
#2 2 2 0 0 4 4 1 50 29
#3 2 2 0 0 4 4 2 50 29
#4 2 2 0 0 4 5 0 50 29
#5 2 2 0 0 4 5 1 50 29
#6 2 2 0 0 4 5 2 50 29
#7 2 2 0 0 5 4 0 50 29
#8 2 2 0 0 5 4 1 50 29
#9 2 2 0 0 5 4 2 50 29
#10 2 2 0 0 5 4 3 50 29
#11 2 2 0 1 4 4 0 40 12
#12 2 2 0 1 4 4 1 40 12
#13 2 2 0 1 4 4 2 40 12
#14 2 2 0 1 4 5 0 40 12
#15 2 2 0 1 4 5 1 40 12
#16 2 2 0 1 4 5 2 40 12
#17 2 2 0 1 5 4 0 40 12
#18 2 2 0 1 5 4 1 40 12
#19 2 2 0 1 5 4 2 40 12
#20 2 2 0 1 5 4 3 40 12
#21 2 2 1 0 4 4 1 31 17
#22 2 2 1 0 4 4 2 31 17
#23 2 2 1 0 4 4 3 31 17
#24 2 2 1 0 4 5 1 31 17
#25 2 2 1 0 4 5 2 31 17
#26 2 2 1 0 4 5 3 31 17
#27 2 2 1 0 5 4 1 31 17
#28 2 2 1 0 5 4 2 31 17
#29 2 2 1 0 5 4 3 31 17
#30 2 2 1 0 5 4 4 31 17
#31 2 2 1 1 4 4 1 7 36
#32 2 2 1 1 4 4 2 7 36
#33 2 2 1 1 4 4 3 7 36
#34 2 2 1 1 4 5 1 7 36
#35 2 2 1 1 4 5 2 7 36
#36 2 2 1 1 4 5 3 7 36
#37 2 2 1 1 5 4 1 7 36
#38 2 2 1 1 5 4 2 7 36
#39 2 2 1 1 5 4 3 7 36
#40 2 2 1 1 5 4 4 7 36
#41 3 2 0 0 5 4 0 13 41
#42 3 2 0 0 5 4 1 13 41
#43 3 2 0 0 5 4 2 13 41
#44 3 2 0 1 5 4 0 27 22
#45 3 2 0 1 5 4 1 27 22
#46 3 2 0 1 5 4 2 27 22
#47 3 2 1 0 5 4 1 49 49
#48 3 2 1 0 5 4 2 49 49
#49 3 2 1 0 5 4 3 49 49
#50 3 2 1 1 5 4 1 47 49
#51 3 2 1 1 5 4 2 47 49
#52 3 2 1 1 5 4 3 47 49
#53 3 2 2 0 5 4 2 23 43
#54 3 2 2 0 5 4 3 23 43
#55 3 2 2 0 5 4 4 23 43
#56 3 2 2 1 5 4 2 4 50
#57 3 2 2 1 5 4 3 4 50
#58 3 2 2 1 5 4 4 4 50
A.K.
- Original Message -
From: arun smartpink...@yahoo.com
To: Joanna Zhang zjoanna2...@gmail.com
Cc: R help r-help@r-project.org
Sent: Saturday, February 16, 2013 11:49 PM
Subject: Re: [R] cumulative sum by group and under some criteria
d2- data.frame()
for (m1 in 2:3) {
for (n1 in 2:2) {
for (x1 in 0:(m1-1)) {
for (y1 in 0:(n1-1)) {
for (m in (m1+2): (7-n1)){
for (n in (n1+2):(9-m)){
for (x in x1:(x1+m-m1)){
d2- rbind(d2,c(m1,n1,x1,y1,m,n,x))
}}}
colnames(d2)-c(m1,n1,x1,y1,m,n,x)
res-do.call(rbind,lapply(2:3,function(m1)
do.call(rbind,lapply(2:2,function(n1)
do.call(rbind,lapply(0:(m1-1),function(x1)
do.call(rbind,lapply(0:(n1-1),function(y1)
do.call(rbind,lapply((m1+2):(7-n1),function(m)
do.call(rbind,lapply((n1+2):(9-m),function(n)
do.call(rbind,lapply(x1:(x1+m-m1), function(x) expand.grid(m1,n1,x1,y1,m,n,x))
)
names(res)- c(m1,n1,x1,y1,m,n,x)
attr(res,out.attrs)-NULL
identical(d2,res)
#[1] TRUE
A.K.
From: Joanna Zhang zjoanna2...@gmail.com
To: arun smartpink...@yahoo.com
Sent: Saturday, February 16, 2013 8:46 PM
Subject: Re: [R] cumulative sum by group and under some criteria
Hi,
What I need is to expand each row by adding several columns and . Let me
restate the question. I have a dataset d, I want to expand it to d2 showed
below.
d-data.frame()
for (m1 in 2:3) {
for (n1 in 2:2) {
for (x1 in 0:(m1-1)) {
for (y1 in 0:(n1-1)) {
d-rbind(d,c(m1,n1,x1,y1))
}
}
}}
colnames(d)-c(m1,n1,x1,y1)
d
m1 n1 x1 y1
1 2 2 0 0
2 2 2 0 1
3 2 2 1 0
4 2 2 1 1
5 3 2 0 0
6 3 2 0 1
7 3 2 1 0
8 3 2 1 1
9 3 2 2 0
10 3 2 2 1
I want to expand it as follows:
for (m in (m1+2): (7-n1){
for (n in (n1+2):(9-m){
for (x in x1:(x1+m-m1){
}}}
so for the first row,
m1 n1 x1 y1
1 2 2 0 0
it should be expanded as
m1 n1 x1 y1 m n x
2 2 0 0 4 4 0
2 2 0 0 4 4 1
2 2 0 0 4 4 2
2 2 0 0 4 5 0
2 2 0 0 4 5 1
2 2 0 0 4 5 2
On Tue, Feb 12, 2013 at 8:19 PM,
[R] system() stdout not recieved
I would like to execute a python script from R and receive the stdout in R. I
have windows xp and R 2.14.2. The test.py script should print hello, it does
work in cmd. Here is a couple tries, thanks for any suggestions!John
system('cmd test.py', intern = TRUE)[1] Microsoft Windows XP [Version
5.1.2600][2] (C)
Copyright 1985-2001 Microsoft Corp.
[3]
[4] C:\\Documents and
Settings\\HP_Administrator.-1\\My Documents\\Dropbox\\Cluster Re-Run
system('cmd test.py')Microsoft Windows XP [Version 5.1.2600](C) Copyright
1985-2001 Microsoft Corp.
C:\Documents and Settings\HP_Administrator.-1\My Documents\Dropbox\Cluster
Re-Run system('test.py')
[[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] compare posterior samples from R2OpenBugs and R function bugs{R2WinBUGS}
Hi all,
I used both OpenBugs and R function bugs{R2WinBUGS} to run a linear mixed
effects model based on the same data set and initial values. I got the same
summary statistics but different posterior samples. However, if I order
these two sets of samples, one is generated from OpenBugs and the other is
generated from R, they turn to be the same. And the samples from R do not
have any autocorrelation. I do not know why and how this R function destroy
the orders of posterior samples. Have anyone ever met this situation
before? Any idea is appreciated.
Thanks,
Jia
sessionInfo()R version 2.15.1 (2012-06-22)
Platform: x86_64-pc-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] R2WinBUGS_2.1-18 BRugs_0.8-0 coda_0.15-2 lattice_0.20-6
loaded via a namespace (and not attached):
[1] grid_2.15.1 tools_2.15.1
[[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.
