[R] R function for censored linear regression
Dear R-useRs, I'm looking for an R-function for censored linear regression. I have the following data x1 - rnorm(100) x2 - rnorm(100) y - x1 + 2*x2 + rnorm(100,0,0.5) stat - rep(1,100) stat[50:100] - 0 data - data.frame(y,x1,x2,stat) y is the dependent variable, x1 and x2 are the independent variables in a linear model. the variable y could be right-censored, this information is in the variable stat, where 1 denotes observed and 0 denotes censored. If stat is 0, then the value in y is the observed right-censored value and could be greater. Using the Tobit-model would not be the right thing here because the Tobit model assumes the same limit for all observations, in my data each value of y[50:100] could have a different limit. If i use linear regression lm1 - lm(y ~ x1 + x2, data=data) summary(lm1) the censoring is not incorporated, so my idea is to use survreg from the survival package library(survival) s1 - survreg(Surv(y, stat) ~ x1 + x2, data, dist='gaussian') summary(s1) my question is, is this the right approach for my aim? Is it right, that here each censored observations could have its own limit? Thanks and best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] survival package - calculating probability to survive a given time
Dear R users, i try to calculate the probabilty to survive a given time by using the estimated survival curve by kaplan meier. What is the right way to do that? as far as is see i cannot use the predict-methods from the survival package? library(survival) set.seed(1) time - cumsum(rexp(1000)/10) status - rbinom(1000, 1, 0.5) ## kaplan meier estimates fit - survfit(Surv(time, status) ~ 1) s - summary(fit) ## 1. possibility to get the probability for surviving 20 units of time ind - findInterval(20, s$time) cbind(s$surv[ind], s$time[ind]) ## 2. possibility to get the probability for surviving 20 units of time ind - s$time = 20 sum(ind) / length(ind) Thanks and best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] gamma glm - using of weights gives error
Dear R-users, i try to use the following code to do a gamma regression glm(x1 / x2 ~ x3 + x4 + x5 + x6 + x7 + x8, family=Gamma(link=log), weights=x2) but here i get the error Error: NA/NaN/Inf in foreign function call (arg 1) In addition: Warning message: step size truncated due to divergence x2 has integer values ranging from 1 to 6. If i do instead glm(x1 / x2 ~ x3 + x4 + x5 + x6 + x7 + x8, family=Gamma(link=log)) without using the weights-argument i get no error. So far i don't really understand what this argument realy does, what is the difference in the results? What can i do to locate the root of the error and how can i avoid this error? Can i to the regression without the weights argument? Thanks and best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Replacing double loop by apply
Dear R-users,
after trying and searching a long time i have the following question.
is it possible to replace to following double loop by some apply calls?
###
m1 - data.frame(v1=factor(letters[1:5]),
v2=factor(letters[2:6]),
v3=factor(letters[3:7]))
m2 - data.frame(v1=factor(letters[3:7]),
v2=factor(letters[1:5]),
v3=factor(letters[4:8]))
ind - matrix(logical(nrow(m2) * ncol(m2)),
nrow = nrow(m2), byrow = TRUE)
for (j in 1:ncol(m2))
{
for (i in 1:nrow(m2))
{
ind[i, j] - any(levels(m1[, j]) == m2[i, j])
}
ind[is.na(ind[, j]), j] - TRUE
m2[!ind[, j], j] - NA
}
###
thanks and best regards
Andreas
--
GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT!
__
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] Replacing double loop by apply
Dear Henrique Dallazuanna,
thank you very much, this really helped me.
by the way, do you thinks it is also possible to do the following loop in R?
thanks and best regards
Andreas
v1 - rnorm(10)
v2 - v1+rnorm(10)
v3 - v2+rnorm(10)
v4 - rnorm(10)
y - ifelse(v10,1,0)
dat - data.frame(v1=v1,v2=v2,v3=v3,v4=v4)
newdat-dat[9:10,]
dat - dat[1:8,]
newdat[2,1] - NA
pred - numeric(2)
vars - c(v1,v3,v4)
y - y[1:8]
for (i in 1:nrow(newdat))
{
ind - which(!is.na(newdat[i, ]) colnames(newdat) %in% vars)
vars.dat - colnames(dat[ind])
formu - as.formula(paste(cbind(y, 1-y) ~ , paste(vars.dat,
collapse = +)))
glm_tmp - glm(formu, family = binomial(link = logit), data = dat[, ind])
pred[i] - predict(glm_tmp, newdat = newdat[i, ind],
type = response)
}
On 14.05.2010 14:42, Henrique Dallazuanna wrote:
Try this:
mapply(function(x, y)x[match(x, y)], m2, m1)
On Fri, May 14, 2010 at 9:23 AM, Andreas Wittmann
andreas_wittm...@gmx.de mailto:andreas_wittm...@gmx.de wrote:
Dear R-users,
after trying and searching a long time i have the following question.
is it possible to replace to following double loop by some apply
calls?
###
m1 - data.frame(v1=factor(letters[1:5]),
v2=factor(letters[2:6]),
v3=factor(letters[3:7]))
m2 - data.frame(v1=factor(letters[3:7]),
v2=factor(letters[1:5]),
v3=factor(letters[4:8]))
ind - matrix(logical(nrow(m2) * ncol(m2)),
nrow = nrow(m2), byrow = TRUE)
for (j in 1:ncol(m2))
{
for (i in 1:nrow(m2))
{
ind[i, j] - any(levels(m1[, j]) == m2[i, j])
}
ind[is.na http://is.na(ind[, j]), j] - TRUE
m2[!ind[, j], j] - NA
}
###
thanks and best regards
Andreas
--
GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT!
__
R-help@r-project.org mailto: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.
--
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40 S 49° 16' 22 O
[[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 using ecm.mix
Dear R-users, i have the following exmple for which i want to use ecm.mix from the mix-package. with da.mix after using em.mix i get the error improper posterior--empty cells, which is not uncommen because of 17 * 5 * 3 = 255 cells. so the next attempt is to use the ecm.mix for the restricted model, but i get errors also. what can i do to get the imputations with the mix package working? maybe it could be reasonable to merge some levels before imputation? best regards Andreas y - matrix(0, nrow=100, ncol=4, byrow=TRUE) y[,1] - sample(17,100,replace=TRUE) y[,2] - sample(5,100,replace=TRUE) y[,3] - sample(3,100,replace=TRUE) y[,4] - rnorm(100) y[,1] - ifelse(rbinom(100,1,0.4)==1, NA, y[,1]) y[,2] - ifelse(rbinom(100,1,0.7)==1, NA, y[,2]) y[,3] - ifelse(rbinom(100,1,0.05)==1, NA, y[,3]) ## first attempt, imputation under unrestricted model s - prelim.mix(y,3) thetahat - em.mix(s) rngseed(1234567) # set random number generator seed newtheta - da.mix(s,thetahat,steps=100) # data augmentation ##ximp - imp.mix(s, newtheta, y) ## da.mix gives improper posterior--empty cells error ## second attempt, imputation under restricted model s - prelim.mix(y,3) margins1 - c(1,2,3) margins2 - c(1,2,0,2,3,0,1,3) design - diag(rep(1,255)) # identity matrix D=no of cells thetahat - ecm.mix(s,margins1,design) # find ML estimate thetahat - ecm.mix(s,margins2,design) # find ML estimate #rngseed(1234567) # random generator seed #newtheta - dabipf.mix(s,margins,design,thetahat,steps=200) #ximp - imp.mix(s,newtheta,stlouis) # impute under newtheta __ 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] Name for factor's levels with contr.sum
Hi R-useRs, after having read http://tolstoy.newcastle.edu.au/R/help/05/07/8498.html with the same topic but five years older. the solution for a contr.sum with names for factor levels for R version 2.10.1 will be to comment out the following line #colnames(cont) - NULL in contr.sum i guess? by the way, with contrasts=FALSE colnames are set, so i don't know what the aim is to avoid this with contrasts=TRUE? best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] enty-wise closest element
Dear R-users,
i have a simple problem maybe, but i don't see the solution. i want to
find the entry-wise closest element of an vector compared with another.
ind1-c(1,4,10)
ind2-c(3,5,11)
for (i in length(ind2):1)
{
print(which.min(abs(ind1-ind2[i])))
}
for ind2[3] it should be ind1[3] 10, for ind2[2] it should be ind1[2] 4
and for ind2[1] it should be ind1[1] 1. but with the for-loop above i
get ind1[3], ind1[2] and ind1[2].
any suggestions are quite welcome.
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] enty-wise closest element
Thank you very much for your quick answers.
Sorry, but i forgot to explain i want to find the closest and smaller
element for each entry, so ind1[3] is smaller as 11 and close to 11,
ind1[2] is smaller as 5 and close to 5 and ind1[1] is smaller than 3 and
close to 3.
best regards
Andreas
On Jan 17, 2010, at 11:00 AM, Andreas Wittmann wrote:
Dear R-users,
i have a simple problem maybe, but i don't see the solution. i want
to find the entry-wise closest element of an vector compared with
another.
ind1-c(1,4,10)
ind2-c(3,5,11)
for (i in length(ind2):1)
{
print(which.min(abs(ind1-ind2[i])))
}
for ind2[3] it should be ind1[3] 10, for ind2[2] it should be ind1[2]
4 and for ind2[1] it should be ind1[1] 1. but with the for-loop above
i get ind1[3], ind1[2] and ind1[2].
any suggestions are quite welcome.
You are failing to use the indices that which.min is providing. (And I
think the closest in ind1 to ind2[1] is not 1, but rather 4, so see if
this looks more responsive to your expectations:
for (i in length(ind2):1)
+ {
+ print( ind1[which.min(abs(ind1-ind2[i]))] )
+ }
[1] 10
[1] 4
[1] 4
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
__
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] parsing protocol of states
Dear R-users, actually i try to parse some state protocols for my work. i an easy stetting the code below works fine, if states are reached only once. in harder settings it could be possible that one state gets visited more times. in this case for me its interesting to see how much waiting time lies between to states on the whole. by the way i didn't use R as a parsing tool so far, so any advice for doing this more effectivly are quite welcome. str01 - 2007-10-12 11:50:05 state B. ,2007-10-12 11:50:05 state C. ,2007-10-12 13:23:24 state D. ,2007-10-12 13:23:43 state E. ,2007-10-14 15:43:19 state F. ,2007-10-14 15:43:20 state E. ,2007-10-14 15:43:25 state G. ,2007-10-14 15:43:32 state H. ,2007-10-14 15:43:41 state I. ,2007-10-14 15:43:47 state F. ,2007-10-14 15:43:47 state G. ,2007-10-14 15:48:08 state H. ,2007-10-16 10:10:20 state J. ,2007-10-19 11:12:54 state K ,2007-10-19 11:17:37 state D. ,2007-10-19 11:17:42 state E. ,2007-10-19 11:17:49 state F. ,2007-10-19 11:17:51 state E. ,2007-10-19 11:17:58 state H. ,2007-10-19 11:18:05 state J. ,2007-10-19 11:21:45 state L. str02 - unlist(strsplit(str01, \\,)) x1 - grep(state B, str02) x2 - grep(state C, str02) x3 - grep(state D, str02) x4 - grep(state E, str02) x5 - grep(state F, str02) x6 - grep(state G, str02) x7 - grep(state H, str02) x8 - grep(state I, str02) x9 - grep(state J, str02) x10 - grep(state K, str02) x11 - grep(state L, str02) t1 - substr(str02[x1], 1, 19) t1 - as.POSIXct(strptime(t1, %Y-%m-%d %H:%M:%S)) t2 - substr(str02[x2], 1, 19) t2 - as.POSIXct(strptime(t2, %Y-%m-%d %H:%M:%S)) t3 - substr(str02[x3], 1, 19) t3 - as.POSIXct(strptime(t3, %Y-%m-%d %H:%M:%S)) t4 - substr(str02[x4], 1, 19) t4 - as.POSIXct(strptime(t4, %Y-%m-%d %H:%M:%S)) t5 - substr(str02[x5], 1, 19) t5 - as.POSIXct(strptime(t5, %Y-%m-%d %H:%M:%S)) t6 - substr(str02[x6], 1, 19) t6 - as.POSIXct(strptime(t6, %Y-%m-%d %H:%M:%S)) t7 - substr(str02[x7], 1, 19) t7 - as.POSIXct(strptime(t7, %Y-%m-%d %H:%M:%S)) t8 - substr(str02[x8], 1, 19) t8 - as.POSIXct(strptime(t8, %Y-%m-%d %H:%M:%S)) t9 - substr(str02[x9], 1, 19) t9 - as.POSIXct(strptime(t9, %Y-%m-%d %H:%M:%S)) t10 - substr(str02[x10], 1, 19) t10 - as.POSIXct(strptime(t10, %Y-%m-%d %H:%M:%S)) t11 - substr(str02[x11], 1, 19) t11 - as.POSIXct(strptime(t11, %Y-%m-%d %H:%M:%S)) as.numeric(difftime(t11, t1, units=days)) ## waiting times between state E and F sum(as.numeric(difftime(t5, t4, units=days))) best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Avoid for-loop in creating data.frames
Dear R-users,
after several tries with lapply and searching the mailing list, i want
to ask, wheter and how it is possibly to avoid the for-loop in the
following piece of code?
set2-as.data.frame(matrix(rnorm(9),ncol=3))
set2[1,1] - NA
set2[3,2] - NA
set2[2,1] - NA
dimnames(set2)[1] - list(c(A,B,C))
r - !is.na(set2)
imp - vector(list, ncol(set2))
for (j in 1:dim(set2)[2])
{
imp[[j]] - as.data.frame(matrix(NA, nrow = sum(!r[,j]), ncol = 1))
dimnames(imp[[j]]) - list(row.names(set2)[r[,j] == FALSE], 1)
}
many thanks and best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Avoid for-loop in creating data.frames
Hi babtiste,
thank you very much for your fast answer. your solution is very good,
but i also need the dimnames as in the for-loop for further calculations.
best regards
Andreas
baptiste auguie wrote:
Hi,
Is the following close enough?
apply(set2, 2, function(x) x[is.na(x)])
HTH,
baptiste
2009/12/10 Andreas Wittmann andreas_wittm...@gmx.de:
Dear R-users,
after several tries with lapply and searching the mailing list, i want to
ask, wheter and how it is possibly to avoid the for-loop in the following
piece of code?
set2-as.data.frame(matrix(rnorm(9),ncol=3))
set2[1,1] - NA
set2[3,2] - NA
set2[2,1] - NA
dimnames(set2)[1] - list(c(A,B,C))
r - !is.na(set2)
imp - vector(list, ncol(set2))
for (j in 1:dim(set2)[2])
{
imp[[j]] - as.data.frame(matrix(NA, nrow = sum(!r[,j]), ncol = 1))
dimnames(imp[[j]]) - list(row.names(set2)[r[,j] == FALSE], 1)
}
many thanks and best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
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] Solve linear program without objective function
Dear R-users, i try to solve to following linear programm in R 0 * x_1 + 2/3 * x_2 + 1/3 * x_3 + 1/3 * x_4 = 0.3 x_1 + x_2 + x_3 + x_4 = 1 x_1, x_2, x_3, x_4 0, x_1, x_2, x_3, x_4 1 as you can see i have no objective function here besides that i use the following code. library(lpSolve) f.obj-c(1,1,1,1) f.con-matrix(c(0,2/3,1/3,1/3, 1,1,1,1, 1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1),nrow=6,byrow=TRUE) f.dir - c(=, =, , , , ) f.rhs - c(0.3, 1, 0, 0, 0, 0) lp (max, f.obj, f.con, f.dir, f.rhs)$solution the problem is, the condition x_1, x_2, x_3, x_4 0 is not fulfilled. Any advice would be very helpful. best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Generate missing data patterns
Dear R-users, i try to generate missing values in a matrix X according to a given missingnes pattern R with the probabilities p per row. X-matrix(rnorm(3*100),ncol=3) ## indicator matrix for missingnes (1 observed, 0 missing) R-matrix(c(1,1,1, 0,0,1, 1,1,0, 0,1,1),ncol=3,byrow=TRUE) ## probabilities for row 1, row 2, row 3 and row 4 p-c(0.375,0.25,0.25,0.125) ## does not exactly what i want, because i get rows ## of missinges pattern which are not in R X[rbinom(100,1,p[1])==1,R[1,]==1] - NA X[rbinom(100,1,p[2])==1,R[2,]==1] - NA X[rbinom(100,1,p[3])==1,R[3,]==1] - NA X[rbinom(100,1,p[4])==1,R[4,]==1] - NA So it would be great if i can get any advice how to do this. I also tried rmultinom or sample but without any success so far. Another question is what to to if want the missing pattern R and simulate a certain amount of missingnes mybe about 10 %? I guess i have to mix the probabilities for each row in such a way to get approximatly the wanted missingnes in percent. best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] predict: remove columns with new levels automatically
Thank you all for the good advice.
Now i did a fast hack, which does want i was looking for, maybe anyone
else finds this usefull
set.seed(0)
x - rnorm(9)
y - x + rnorm(9)
training - data.frame(x=x, y=y,
z1=c(rep(A, 3), rep(B, 3), rep(C, 3)),
z2=c(rep(F, 4), rep(G, 5)))
test - data.frame(x=t-rnorm(1), y=t+rnorm(1), z1=D, z2=F)
`predict.drop` - function(f, dat, newdat)
{
datlev - vector(list, ncol(dat))
newdatlev - vector(list, ncol(newdat))
`filllevs` - function(dat, veclev)
{
for (j in 1:ncol(dat))
{
if (is.factor(dat[,j]))
veclev[[j]] - levels(dat[,j])
else
veclev[[j]] - NULL
}
return(veclev)
}
datlev - filllevs(dat, datlev)
newdatlev - filllevs(newdat, newdatlev)
if (ncol(dat) == ncol(newdat))
{
drop - logical(ncol(dat))
names(drop) - colnames(dat)
for (j in 1:ncol(dat))
{
if (!is.null(datlev[[j]]))
{
if (!(newdatlev[[j]] %in% datlev[[j]]))
drop[j] - TRUE
}
}
}
else
stop(dat and newdat must have the same column length!)
m - lm(formula(f), data=dat[,(1:ncol(dat))[!drop]])
p - predict(m, newdat)
return(list(drop=drop, p=p))
}
predict.drop(x ~ ., training, test)
best regards
Andreas
David Winsemius wrote:
On Nov 25, 2009, at 1:48 AM, Andreas Wittmann wrote:
Sorry for my bad description, i don't want get a constructed
algorithm without own work. i only hoped to get some advice how to do
this. i don't want to predict any sort of data, i reference only to
newdata which variables are the same as in the model data. But if
factors in the data than i can by possibly that the newdata has a
level which doesn't exist in the original data.
So i have to compare each factor in the data and in the newdata and
if the newdata has a levels which is not in the original data and
drop this variable and do compute the model and prediction again.
I thought this problem is quite common and i can use an algorithm
somebody has already implemented.
best regards
Andreas
If you use str to look at the lm1 object you will find at the bottom a
list called x:
lm1$x will show you the factors that were present in variables at the
time of the model creation
lm1$x
$z
[1] A B C
New testing scenario good level and bad level:
test - data.frame(x=t-rnorm(2), y=t+rnorm(2), z=c(B, D) )
lm1 - lm(x ~ ., data=training)
predict(lm1, subset(test, z %in% lm1$x$z) ) # get prediction for
good level only
1
0.4225204
Original-Nachricht
Datum: Wed, 25 Nov 2009 00:48:59 -0500
Von: David Winsemius dwinsem...@comcast.net
An: Andreas Wittmann andreas_wittm...@gmx.de
CC: r-help@r-project.org
Betreff: Re: [R] predict: remove columns with new levels automatically
On Nov 24, 2009, at 2:24 PM, Andreas Wittmann wrote:
Dear R-users,
in the follwing thread
http://tolstoy.newcastle.edu.au/R/help/03b/3322.html
the problem how to remove rows for predict that contain levels which
are not in the model.
now i try to do this the other way round and want to remove columns
(variables) in the model which will be later problematic with new
levels for prediction.
## example:
set.seed(0)
x - rnorm(9)
y - x + rnorm(9)
training - data.frame(x=x, y=y, z=c(rep(A, 3), rep(B, 3),
rep(C, 3)))
test - data.frame(x=t-rnorm(1), y=t+rnorm(1), z=D)
lm1 - lm(x ~ ., data=training)
## prediction does not work because the variable z has the new level
D
predict(lm1, test)
## solution: the variable z is removed from the model
## the prediction happens without using the information of variable z
lm2 - lm(x ~ y, data=training)
predict(lm2, test)
How can i autmatically recognice this and calculate according to this?
Let me get this straight. You want us to predict in advance (or more
accurately design an algorithm that can see into the future and work
around) any sort of newdata you might later construct
--
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
--
Preisknaller: GMX DSL Flatrate für nur 16,99 Euro/mtl.!
http://portal.gmx.net/de/go/dsl02
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
__
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] predict: remove columns with new levels automatically
Dear R-users, in the follwing thread http://tolstoy.newcastle.edu.au/R/help/03b/3322.html the problem how to remove rows for predict that contain levels which are not in the model. now i try to do this the other way round and want to remove columns (variables) in the model which will be later problematic with new levels for prediction. ## example: set.seed(0) x - rnorm(9) y - x + rnorm(9) training - data.frame(x=x, y=y, z=c(rep(A, 3), rep(B, 3), rep(C, 3))) test - data.frame(x=t-rnorm(1), y=t+rnorm(1), z=D) lm1 - lm(x ~ ., data=training) ## prediction does not work because the variable z has the new level D predict(lm1, test) ## solution: the variable z is removed from the model ## the prediction happens without using the information of variable z lm2 - lm(x ~ y, data=training) predict(lm2, test) How can i autmatically recognice this and calculate according to this? Thanks Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] predict: remove columns with new levels automatically
Sorry for my bad description, i don't want get a constructed algorithm without own work. i only hoped to get some advice how to do this. i don't want to predict any sort of data, i reference only to newdata which variables are the same as in the model data. But if factors in the data than i can by possibly that the newdata has a level which doesn't exist in the original data. So i have to compare each factor in the data and in the newdata and if the newdata has a levels which is not in the original data and drop this variable and do compute the model and prediction again. I thought this problem is quite common and i can use an algorithm somebody has already implemented. best regards Andreas Original-Nachricht Datum: Wed, 25 Nov 2009 00:48:59 -0500 Von: David Winsemius dwinsem...@comcast.net An: Andreas Wittmann andreas_wittm...@gmx.de CC: r-help@r-project.org Betreff: Re: [R] predict: remove columns with new levels automatically On Nov 24, 2009, at 2:24 PM, Andreas Wittmann wrote: Dear R-users, in the follwing thread http://tolstoy.newcastle.edu.au/R/help/03b/3322.html the problem how to remove rows for predict that contain levels which are not in the model. now i try to do this the other way round and want to remove columns (variables) in the model which will be later problematic with new levels for prediction. ## example: set.seed(0) x - rnorm(9) y - x + rnorm(9) training - data.frame(x=x, y=y, z=c(rep(A, 3), rep(B, 3), rep(C, 3))) test - data.frame(x=t-rnorm(1), y=t+rnorm(1), z=D) lm1 - lm(x ~ ., data=training) ## prediction does not work because the variable z has the new level D predict(lm1, test) ## solution: the variable z is removed from the model ## the prediction happens without using the information of variable z lm2 - lm(x ~ y, data=training) predict(lm2, test) How can i autmatically recognice this and calculate according to this? Let me get this straight. You want us to predict in advance (or more accurately design an algorithm that can see into the future and work around) any sort of newdata you might later construct -- David Winsemius, MD Heritage Laboratories West Hartford, CT -- Preisknaller: GMX DSL Flatrate für nur 16,99 Euro/mtl.! http://portal.gmx.net/de/go/dsl02 __ 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] recode according to old levels
Dear R-users, i try to recode a factor according to old levels F - factor(sample(c(rep(A, 4), rep(B,2), rep(C,5 recode(F, levels(F)[c(1,3)]='X'; else='Y') i tried to work with eval or expression around levels(F)[c(1,3)], but nothing seems to work. Many thanks if anyone could tell me what i've missed and what's the problem here. best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] predict missing values with svm
Dear R-Users,
i want to use the function svm of the e1071 package to predict missing data
data(iris)
## create missing completely at random data
for (i in 1:5)
{
mcar - rbinom(dim(iris)[1], size=1, prob=0.1)
iris[mcar == 1, i] - NA
}
ok - complete.cases(iris)
model - svm(Species ~ ., data=iris[ok,])
## try to predict the missing values for Species
## neither
pred - predict(model, iris[5])
## nor
pred - predict(model, iris[!ok, -5])
## seems to work
Many thanks if anyone could tell me what i do wrong and what is the
problem here.
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] mice and Date-Time Classes
Dear R-Users, i want to use the function mice of the mice package with data, that contains dates, but this gives an error x - Sys.time() dat - data.frame(dates=(1:10)*60*60*24+x, size=rnorm(10)) dat$dates[c(3,7)]-NA mice(dat) Fehler in check.imputationMethod(imputationMethod, defaultImputationMethod, : The following functions were not found: mice.impute. I guess the mice package is not implemented for the POSIXt data class yet. A workaround could be to build time intervals and use mice with that and transform it back afterwards. Or is there a more efficient way? Many thanks if anyone could tell me what i do wrong and what is the problem here. best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Runtime Error with multinom
Dear R-users, i try to fit a multinomial model in order to get an imputation for a missing value in factor1. library(nnet) factor1 - factor(c(a,b,c,d)) factor2 - factor(c(e,f,g,h)) size - c(3,8,2,1) factor1[3] - NA Z-ifelse(is.na(factor1), 0, 1) assign(data, cbind.data.frame(factor1,factor2,size),pos=1) multinom(formula(data),data=data[Z,]) when entering the last line i get a runtime error and R crashes down. my system is windows xp and R 2.9.1, i tried it also with ubuntu 9.04 and R.2.8.1 but i get a quite similar error. Many thanks if anyone could tell me what i do wrong and what is the problem here. best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] integrate with large parameters
Dear R-users,
i have to integrate the following function
`fun1` - function (a, l1, l2)
{
exp(log(l1) * (a - 1) - l2 * lgamma(a))
}
but if l1 is large, i get the non-finite function value error, so my
idea is to rescale with exp(-l1)
`fun2` - function (a, l1, l2)
{
exp(log(l1) * (a - 1) - l2 * lgamma(a) - l1)
}
but it seems this doesn't solve the problem, when l1=1.0e80. maybe i
have to iterate between large values in integrate and large values in
exp(l1) and then use the normal or the rescaled version?
l1-10; l2-10;
integrate(Vectorize(function(a) {fun1(a, l1=l1, l2=l2)}), 1, Inf)$val
[1] 11.65311
integrate(Vectorize(function(a) {fun2(a, l1=l1, l2=l2)}), 1, Inf)$val
* exp(l1)
[1] 11.65312
l1-1.0e12; l2=50;
integrate(Vectorize(function(a) {fun1(a, l1=l1, l2=l2)}), 1, Inf)$val
[1] 929558588152
integrate(Vectorize(function(a) {fun2(a, l1=l1, l2=l2)}), 1, Inf)$val
* exp(l1)
[1] NaN
l1-1.0e20; l2=50;
integrate(Vectorize(function(a) {fun1(a, l1=l1, l2=l2)}), 1, Inf)$val
[1] 5.005192e+24
integrate(Vectorize(function(a) {fun2(a, l1=l1, l2=l2)}), 1, Inf)$val
* exp(l1)
[1] NaN
l1-1.0e80; l2=50;
integrate(Vectorize(function(a) {fun1(a, l1=l1, l2=l2)}), 1, Inf)$val
Fehler in integrate(Vectorize(function(a) { : non-finite function value
integrate(Vectorize(function(a) {fun2(a, l1=l1, l2=l2)}), 1, Inf)$val
* exp(l1)
[1] NaN
Many thanks if you have any advice for me!
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] integrate lgamma from 0 to Inf
Dear R users,
i try to integrate lgamma from 0 to Inf. But here i get the message roundoff
error is detected in the extrapolation table, if i use 1.0e120 instead of Inf
the computation works, but this is against the suggestion of integrates help
information to use Inf explicitly. Using stirlings approximation doesnt bring
the solution too.
## Stirlings approximation
lgammaApprox - function(x)
{
0.5*log(2*pi)+(x-(1/2))*log(x)-x
}
integrate(lgamma, lower = 0, upper = 1.0e120)
integrate(lgammaApprox, lower = 0, upper = 1.0e120)
integrate(lgamma, lower = 0, upper = 1.0e120)
1.374051e+242 with absolute error 3.2e+235
integrate(lgammaApprox, lower = 0, upper = 1.0e120)
1.374051e+242 with absolute error 3.2e+235
integrate(lgamma, lower = 0, upper = Inf)
integrate(lgammaApprox, lower = 0, upper = Inf)
integrate(lgamma, lower = 0, upper = Inf)
Fehler in integrate(lgamma, lower = 0, upper = Inf) :
roundoff error is detected in the extrapolation table
integrate(lgammaApprox, lower = 0, upper = Inf)
Fehler in integrate(lgammaApprox, lower = 0, upper = Inf) :
roundoff error is detected in the extrapolation table
Many thanks if you have any advice for me!
best regards
Andreas
--
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Continue R CMD BATCH on error
Dear R useRs,
after searching r-help and r-manuals for about one hour i have the
following, probably easy question for you.
i have the following R-code, in the file test01.R
`fun1` - function(x)
{
x - x + 2
if(x == 5)
stop(failure)
return(x)
}
`fun2` - function(x)
{
x - x + 4
return(x)
}
x - 1:10
val1 - val2 - numeric(10)
for(i in 1:10)
{
val1[i] - fun1(x[i])
}
for(i in 1:10)
{
val2[i] - fun2(x[i])
}
then i want to do
R CMD BATCH test01.R
the result file test01.Rout does not contain the computation of val2 and
it stops in the for loop of val1.
how can i avoid this and continue the computation on error?
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Different values for double integral
Dear R useRs,
i have the function f1(x, y, z) which i want to integrate for x and y.
On the one hand i do this by first integrating for x and then for y, on
the other hand i do this the other way round and i wondering why i
doesn't get the same result each way?
z - c(80, 20, 40, 30)
f1 - function(x, y, z) {dgamma(cumsum(z)[-length(z)], shape=x, rate=y)}
g1 - function(y, z) {integrate(function(x) {f1(x=x, y=y, z=z)}, 0.1,
0.5)$value}
g2 - function(x, z) {integrate(function(y) {f1(x=x, y=y, z=z)}, 0.1,
0.5)$value}
integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 5.909943e-09
integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 5.978334e-09
If you have any advice what is wrong here, i would be very thankful.
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Different values for double integral
Thank you all, for the very helpful advice.
i want to estimate the parameters omega and beta of the gamma-nhpp-model
with numerical integration.
so one step in order to do this is to compute the normalizing constant,
but as you see below i get different values
## some reliability data, theses are times between events like a
failures during software test.
`s` -
c(8100, 4800, 900, 450, 450, 6000, 2400, 2100, 2100,
1260, 1200, 300, 9000, 600, 2400, 600, 15060, 120, 360,
1200, 300, 1200, 1200, 3300, 19800, 3000, 600, 9600,
8400, 8100, 15600, 1800, 5400, 3900, 2400, 1200, 79500,
9000, 48900)
## likelihood function of the NHPP-gamma reliability model
`likelihood` - function(s, omega, beta, alpha=1)
{
me-length(s)-1
te-cumsum(s)[length(s)] ## endpoint of observation
omega^me*prod(dgamma(cumsum(s)[-length(s)], shape=alpha, rate=beta)) *
exp(-omega*pgamma(te, shape=alpha, rate=beta))
}
## normalizing constant first integrating omega then beta
`normConst1` - function(s, alpha=1, lowBeta=-Inf, uppBeta=Inf,
lowOmega=-Inf, uppOmega=Inf)
{
g - function(beta, ...)
{
integrate(function(omega) {likelihood(s=s, omega, beta,
alpha=alpha)}, lowOmega, uppOmega)$value
}
val - integrate(Vectorize(function(beta) {g(beta, lowOmega=lowOmega,
uppOmega=uppOmega, s=s,
alpha=alpha)}), lowBeta,
uppBeta)$value
return(val)
}
## normalizing constant first integrating beta then omega
`normConst2` - function(s, x, alpha=1, lowBeta=-Inf, uppBeta=Inf,
lowOmega=-Inf, uppOmega=Inf)
{
g - function(omega, ...)
{
integrate(function(beta) {likelihood(s=s, omega, beta,
alpha=alpha)}, lowBeta, uppBeta)$value
}
val - integrate(Vectorize(function(omega) {g(omega, lowBeta=lowBeta,
uppBeta=uppBeta, s=s,
alpha=alpha)}), lowOmega,
uppOmega)$value
return(val)
}
## integration limits were choosen by taking the quantiles of beta and
omega
## of a mcmc run
normConst1(s=s, lowBeta=2.8e-06, uppBeta=2.8e-05,
lowOmega=12.5, uppOmega=90)
[1] 5.131021e-162
normConst2(s=s, lowBeta=2.8e-06, uppBeta=2.8e-05,
lowOmega=12.5, uppOmega=90)
[1] 5.246008e-158
## i get normalizing constants with different values
best regards and many thanks
Andreas
Prof Brian Ripley wrote:
More generally, if you want to do a two-dimensional integral, you will
do better to us a 2D integration algorithm, such as those in package
'adapt'.
Also, these routines are somewhat sensitive to scaling, so if the
correct answer is around 5e-9, you ought to rescale. You seem to be
in the far right tail of your gamma distribution (but as you are not
integrating over z, pgamma is not appropriate).
More specifically, cumsum(z)[-length(z)] is constant ( c(80,100, 140)
) are can be done once. But without knowing the intention of f1, it
is impossible to show you better code.
On Sat, 24 Jan 2009, Duncan Murdoch wrote:
On 24/01/2009 5:23 AM, Andreas Wittmann wrote:
Dear R useRs,
i have the function f1(x, y, z) which i want to integrate for x and
y. On the one hand i do this by first integrating for x and then for
y, on the other hand i do this the other way round and i wondering
why i doesn't get the same result each way?
z - c(80, 20, 40, 30)
f1 - function(x, y, z) {dgamma(cumsum(z)[-length(z)], shape=x,
rate=y)}
g1 - function(y, z) {integrate(function(x) {f1(x=x, y=y, z=z)},
0.1, 0.5)$value}
g2 - function(x, z) {integrate(function(y) {f1(x=x, y=y, z=z)},
0.1, 0.5)$value}
integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 5.909943e-09
integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 5.978334e-09
If you have any advice what is wrong here, i would be very thankful.
It looks as though your f1 returns a vector result whose length will
be the max of length(z)-1, length(x), and length(y): that's not
good, when you don't have control over the lengths of x and y. I'd
guess that's your problem. I don't know what your intention was, but
if the lengths of x and y were 1, I think it should return a length 1
result.
More geneally, integrate() does approximate integration, and it may
be that you won't get identical results from switching the order even
if you fix the problems above.
Finally, if this is the real problem you're working on, you can use
the pgamma function to do your inner integral: it will be faster and
more accurate than integrate.
Duncan Murdoch
__
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
[R] transform R to C
Dear R users,
i would like to transform the following function from R-code to C-code and call
it from R in order to speed up the computation because in my other functions
this function is called many times.
`dgcpois` - function(z, lambda1, lambda2)
{
`f1` - function(alpha, lambda1, lambda2)
return(exp(log(lambda1) * (alpha - 1) - lambda2 * lgamma(alpha)))
`f2` - function(lambda1, lambda2)
return(integrate(f1, lower=1, upper=Inf, lambda1=lambda1,
lambda2=lambda2)$value)
return(exp(log(lambda1) * z - lambda2 * lgamma(z + 1) -
log(f2(lambda1=lambda1, lambda2=lambda2
}
In order to do this i read for example dgamma.c or dpois.c but for my it seems
rather cryptic and so i would like get some advice in writing the c-function.
First of all i think i cannot call the integrate r-function from c so i have to
use the internal c-function of integrate here?
Thanks,
Andreas
--
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Finding the first value without warning in a loop
Dear R useRs,
with the following piece of code i try to find the first value which can
be calculated without warnings
`test` - function(a)
{
repeat
{
## hide warnings
suppressWarnings(log(a))
if (exists(last.warning, envir = .GlobalEnv))
{
a - a + 0.1
## clear existing warnings
rm(last.warning, envir = .GlobalEnv)
}
if(a 5 || !exists(last.warning, envir = .GlobalEnv))
break
}
return(a)
}
if i run this with test(-3), i would expect a=0 as return value.
Is it also possible to hide warnings during my function, i guess i use
suppressWarnings in a wrong way here?
Thanks and best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] ESS Toolbar missing after Ubuntu Update
Dear R useRs, yesterday i updated my system from ubuntu 8.04 to 8.10. I use emacs- snapshot, this is emacs 23.0.60.1 and ess 5.3.8. Before the update i had when starting emacs with an R file an ess-toolbar with little icons to start R or to evaluate a line or a region of my R file, but no this toolbar is lost and i don't know how to get it again. I tried a lot of changes in the ess options but without any success. searching with google could not solve my problem. If you have any advice for me i would be very thankful best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] folded normal distribution in R
Dear R useRs, i wanted to ask if the folded normal destribution (Y = abs(X) with X normal distributed) with density and random number generator is implemented in R or in any R-related package so far? Maybe i can use the non-central chi-square distribution and rchisq(n, df=1, ncp0) here? Thanks and best regards Andreas __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R CMD SHLIB: file not recognized: File format not recognized
Hello Dirk,
thank you for your chick answer.
I tried another file and there it works. so i removed all files
which were created during of the compilation of add.c in windows
and so i could compile it under ubuntu too. During the windows
compilation there is some *.o file which is created during the
compilation, if i delete it and try the compilation under ubuntu
everything works fine, don't know why?
best regards
Andreas
Dirk Eddelbuettel wrote:
On Sun, Oct 19, 2008 at 01:27:06AM +0200, Andreas Wittmann wrote:
Dear R useRs,
on ubuntu 8.04 i try to create a shared object out of a c-file
this is
// add.c
#include Rinternals.h
SEXP addiere(SEXP a, SEXP b)
{
int i, n;
n = length(a);
for (i = 0; i n; i++)
REAL(a)[i] += REAL(b)[i];
return(a);
}
in terminal i type
R CMD SHLIB add.c
and get
gcc -std=gnu99 -shared -o add.so add.o -L/usr/lib/R/lib -lR
add.o: file not recognized: File format not recognized
my gcc version is 4.2.3 and my R version is 2.7.2.
I can't replicate that on Ubuntu 8.04:
[EMAIL PROTECTED]:/tmp$ R CMD SHLIB add.c
gcc -std=gnu99 -I/usr/share/R/include -fpic -g -O2 -c add.c -o add.o
gcc -std=gnu99 -shared -o add.so add.o -L/usr/lib/R/lib -lR
[EMAIL PROTECTED]:/tmp$
Works fine here. Have you compiled other files? Or are you maybe
missing some -dev packages?
Dirk
Searching R-help, Writing R Extensions and R Installation and
Administration guide
i don't get any idea whats wrong here?
Creating a dll-file on windows xp with the same c-file works fine.
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] R CMD SHLIB: file not recognized: File format not recognized
Dear R useRs,
on ubuntu 8.04 i try to create a shared object out of a c-file
this is
// add.c
#include Rinternals.h
SEXP addiere(SEXP a, SEXP b)
{
int i, n;
n = length(a);
for (i = 0; i n; i++)
REAL(a)[i] += REAL(b)[i];
return(a);
}
in terminal i type
R CMD SHLIB add.c
and get
gcc -std=gnu99 -shared -o add.so add.o -L/usr/lib/R/lib -lR
add.o: file not recognized: File format not recognized
my gcc version is 4.2.3 and my R version is 2.7.2.
Searching R-help, Writing R Extensions and R Installation and
Administration guide
i don't get any idea whats wrong here?
Creating a dll-file on windows xp with the same c-file works fine.
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] compute posterior mean by numerical integration
Dear R useRs,
i try to compute the posterior mean for the parameters omega and beta
for the following
posterior density. I have simulated data where i know that the true
values of omega=12
and beta=0.01. With the function postMeanOmega and postMeanBeta i wanted
to compute
the mean values of omega and beta by numerical integration, but instead
of omega=12
and beta=0.01 i get omega=11.49574 and beta=1.330105. I don't know what
is wrong
in my implementation, i guess the computation by numerical integration
strongly depends
on the integration limits low, upw, lob and upb, but what are good
choices for these to get
reasonable results?
###
## simulated data with omega=12, beta=0.01
data - c(8, 1, 6, 14, 1, 0, 16, 37, 15, 17, 6, 57)
## log likelihood function
loglike - function(t, omega, beta)
{
n - length(t)-1
res - n * log(omega) + n * log(beta) - beta *
sum(cumsum(t[-length(t)])) -
omega * (1-exp(-beta * cumsum(t)[n]))
return(res)
}
## log prior density
prior - function(omega, beta, o1, o2, b1, b2)
{
if(o1 o2 b1 b2)
res - dgamma(omega, o1, o2, log=T) + dgamma(beta, b1, b2, log=T)
else
res - 0 ## noninformative prior
return(res)
}
## log posterior density
logPost - function(t, omega, beta, o1, o2, b1, b2)
{
res - loglike(t, omega, beta) + prior(omega, beta, o1, o2, b1, b2)
return(res)
}
## posterior normalizing constant
PostNormConst - function(t, low, upw, lob, upb, o1, o2, b1, b2)
{
g - function(beta, ...)
{
integrate(function(omega) {logPost(t=t, omega, beta, o1=o1, o2=o2,
b1=b1, b2=b2)}, low, upw)$value
}
res - integrate(Vectorize(function(beta) {g(beta, low=low, upw=upw, t=t,
o1=o1, o2=o2, b1=b1,
b2=b2)}),
lob, upb)$value
return(res)
}
## posterior mean omega
postMeanOmega - function(t, norm, low, upw, lob, upb, o1, o2, b1, b2)
{
g - function(beta, ...)
{
integrate(function(omega) {logPost(t=t, omega, beta, o1=o1, o2=o2,
b1=b1, b2=b2) * omega}, low,
upw)$value
}
tmp - integrate(Vectorize(function(beta) {g(beta, low=low, upw=upw, t=t,
o1=o1, o2=o2, b1=b1,
b2=b2)}),
lob, upb)$value
res - tmp / norm
return(res)
}
## posterior mean beta
postMeanBeta - function(t, norm, low, upw, lob, upb, o1, o2, b1, b2)
{
g - function(beta, ...)
{
integrate(function(omega) {logPost(t=t, omega, beta, o1=o1, o2=o2,
b1=b1, b2=b2) * beta}, low,
upw)$value
}
tmp - integrate(Vectorize(function(beta) {g(beta, low=low, upw=upw, t=t,
o1=o1, o2=o2, b1=b1,
b2=b2)}),
lob, upb)$value
res - tmp / norm
return(res)
}
low - 3; upw - 20;
lob - 0; upb - 2;
norm - PostNormConst(t=data, low=low, upw=upw, lob=lob, upb=upb, o1=0,
o2=0, b1=0, b2=0)
postMeanOmega(t=data, norm=norm, low=low, upw=upw, lob=lob, upb=upb,
o1=0, o2=0, b1=0, b2=0)
postMeanBeta(t=data, norm=norm, low=low, upw=upw, lob=lob, upb=upb,
o1=0, o2=0, b1=0, b2=0)
###
If you have any advice, i would be very thankful,
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Computation of contour values - Speeding up computation
Thank you very much, yes maybe its worth working on it.
best regards
Andreas
Uwe Ligges wrote:
Andreas Wittmann wrote:
Dear R useRs,
i have the following code to compute values needed for a contour plot
myContour - function(a, b, plist, veca, vecb, dim)
{
tmpb - seq(0.5 * b, 1.5 * b, length=dim)
tmpa - seq(0.5 * a, 1.5 * a, length=dim)
z - matrix(0, nrow=dim, ncol=dim)
for(i in 1:dim)
{
for(j in 1:dim)
{
z[i, j] - posteriorPdf(a=tmpa[j], b=tmpb[i],
plist=plist, veca=veca, vecb=vecb) }
}
}
posteriorPdf - function(a, b, plist, veca, vecb)
{
res - sum(plist[, 1] *exp(vecb[, 1] *
log(vecb[, 2]) + (vecb[, 1] - 1.0) * log(b) - vecb[,
2] * b - lgamma(vecb[, 1])) *
exp(veca[, 1] * log(veca[, 2]) + (veca[, 1] - 1.0) *
log(a) - veca[, 2] * a - lgamma(veca[, 1])))
return(res) }
plist - matrix(0, 100, 3)
plist[, 1] - runif(100) veca - vecb - matrix(0, 100, 2)
veca[, 1] - seq(20, 50, len=100)
veca[, 2] - seq(10, 20, len=100)
vecb[, 1] - seq(50, 200, len=100)
vecb[, 2] - seq(1000, 40, len=100)
myContour(a=20, b=0.01, plist=plist, veca=veca, vecb=vecb, dim=50)
this is part of my other computations which i do with R. Here i
recognized, that my functions myContour and posteriorPdf took a long
time of my computations. The key to speed this up is to avoid the two
for-loops in myContour, i think. I tried a lot to do this with apply
or something like that, but i didn't get it.
If you have any advice how i can to this computations fast, i would
be very thankful, one idea is to use external c-code?
It takes 0.8 seconds on my machine. Not worth working on it, is it?
Your problem is that you are applying many calculations for all
iterations of the inner loop, even if the result won't change, example:
lgamma(veca[, 1]) will be calculated dim^2 times!
Hence you can improve your loop considerably:
myContour - function(a, b, plist, veca, vecb, dim)
{
tmpb - seq(0.5 * b, 1.5 * b, length=dim)
tmpa - seq(0.5 * a, 1.5 * a, length=dim)
z - matrix(0, nrow=dim, ncol=dim)
plist1 - plist[, 1]
vecb1l2 - vecb[, 1] * log(vecb[, 2])
vecb11 - vecb[, 1] - 1
vecb1lg - lgamma(vecb[, 1])
vecb2 - vecb[, 2]
veca1l2 - veca[, 1] * log(veca[, 2])
veca11 - veca[, 1] - 1
veca2 - veca[, 2]
veca1lg - lgamma(veca[, 1])
for(i in 1:dim)
{
for(j in 1:dim)
{
z[i, j] - sum(plist1 * exp(vecb1l2 + vecb11 * log(tmpb[i]) -
vecb2 * tmpb[i] - vecb1lg) * exp(veca1l2 + veca11 * log(tmpa[j]) -
veca2 * tmpa[j] - veca1lg))
}
}
z
}
Uwe Ligges
best regards
Andreas --
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
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] lower and upper limits in integrate as vectors
Dear R useRs,
i try to integrate the following function for many values
integrand - function(z)
{
return(z * z)
}
i do this with a for-loop
for(i in 2:4)
{
z - integrate(integrand, i-1, i)$value
cat(z, z, \n)
}
to speed up the computation for many values i tried vectors
in integrate to do this.
vec1-1:3
vec2-2:4
integrate(Vectorize(integrand), vec1, vec2)$value
but here it seems the integration works only for the first values of
vec1 and vec2.
If you have any advice how i can to this computations with vectors or
something like that, i would be very thankful,
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] lower and upper limits in integrate as vectors
dear baptiste,
thank you very much, that was excatly what i was looking for:-)
best regards
Andreas
baptiste auguie wrote:
Hi,
I think you want to Vectorize integrate(), not integrand(). Here's a
way using mapply,
integrand - function(z)
{
return(z * z)
}
vec1-1:3
vec2-2:4
mapply(integrate, lower=vec1, upper=vec2, MoreArgs=list(f=integrand) )
baptiste
On 20 Sep 2008, at 13:08, Andreas Wittmann wrote:
Dear R useRs,
i try to integrate the following function for many values
integrand - function(z)
{
return(z * z)
}
i do this with a for-loop
for(i in 2:4)
{
z - integrate(integrand, i-1, i)$value
cat(z, z, \n)
}
to speed up the computation for many values i tried vectors
in integrate to do this.
vec1-1:3
vec2-2:4
integrate(Vectorize(integrand), vec1, vec2)$value
but here it seems the integration works only for the first values of
vec1 and vec2.
If you have any advice how i can to this computations with vectors or
something like that, i would be very thankful,
best regards
Andreas
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
_
Baptiste Auguié
School of Physics
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK
Phone: +44 1392 264187
http://newton.ex.ac.uk/research/emag
__
__
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] certain number of equations in function depending on parameter
Hello everybody,
i have the following problem to write a function which recognizes depending
on the parameter-inputs how many equations for the calculation in the function
are needed.
Here is an example of my problem:
myfun - function(a, b, c, d)
{
k - length(a)
#here d = 3 for example, but how can i dynamically controll
#my function and tell her to build equations eq1 to eq5 if d = 5?
eq1 - function(a, b, y)
{
c[k-1] - a[k-1] + b * y
}
eq2 - function(a, b, y)
{
c[k-2] - a[k-2] + b * y
}
eq3 - function(a, b, y)
{
c[k-3] - a[k-3] + b * y
}
eq4 - function(a, b, z)
{
1 - sum(c(eq1(z), eq2(z), eq3(z), z))
}
sol - uniroot(eq4, lower=0, upper=1)
}
I hope my problems is explained clear enough. I would be very happy
if you can give me some advice.
best regards
Andreas
--
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
