Re: [R] Sensitivity and Specificity
MKclass::perfMeasures(predict_testing, truth = testing$case, namePos = 1) should also work and computes 80 performance measures. Best Matthias Am 25.10.22 um 06:42 schrieb Jin Li: Hi Greg, This can be done by: spm::pred.acc(testing$case, predict_testing) It will return both sensitivity and specificity, along with a few other commonly used measures. Hope this helps, Jin On Tue, Oct 25, 2022 at 6:01 AM Rui Barradas wrote: Às 16:50 de 24/10/2022, greg holly escreveu: Hi Michael, I appreciate your writing. Here are what I have after; predict_testing <- ifelse(predict > 0.5,1,0) head(predict) 1 2 3 5 7 8 0.29006984 0.28370507 0.10761993 0.02204224 0.12873872 0.08127920 # Sensitivity and Specificity sensitivity<-(predict_testing[2,2]/(predict_testing[2,2]+predict_testing[2,1]))*100 Error in predict_testing[2, 2] : incorrect number of dimensions sensitivity function (data, ...) { UseMethod("sensitivity") } specificity<-(predict_testing[1,1]/(predict_testing[1,1]+predict_testing[1,2]))*100 Error in predict_testing[1, 1] : incorrect number of dimensions specificity function (data, ...) { UseMethod("specificity") } On Mon, Oct 24, 2022 at 10:45 AM Michael Dewey wrote: Rather hard to know without seeing what output you expected and what error message you got if any but did you mean to summarise your variable predict before doing anything with it? Michael On 24/10/2022 16:17, greg holly wrote: Hi all R-Help , After partitioning my data to testing and training (please see below), I need to estimate the Sensitivity and Specificity. I failed. It would be appropriate to get your help. Best regards, Greg inTrain <- createDataPartition(y=data$case, p=0.7, list=FALSE) training <- data[ inTrain,] testing <- data[-inTrain,] attach(training) #model training and prediction data_training <- glm(case ~ age+BMI+Calcium+Albumin+meno_1, data = training, family = binomial(link="logit")) predict <- predict(data_training, data_predict = testing, type = "response") predict_testing <- ifelse(predict > 0.5,1,0) # Sensitivity and Specificity sensitivity<-(predict_testing[2,2]/(predict_testing[2,2]+predict_testing[2,1]))*100 sensitivity specificity<-(predict_testing[1,1]/(predict_testing[1,1]+predict_testing[1,2]))*100 specificity [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. -- Michael http://www.dewey.myzen.co.uk/home.html [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. Hello, Instead of computing by hand, why not use package caret? tbl <- table(predict_testing, testing$case) caret::sensitivity(tbl) caret::specificity(tbl) Hope this helps, Rui Barradas __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] confidence intervals for the difference between group means
you could try: library(MKinfer) meanDiffCI(a, b, boot = TRUE) Best Matthias Am 04.08.20 um 16:08 schrieb varin sacha via R-help: Dear R-experts, Using the bootES package I can easily calculate the bootstrap confidence intervals of the means like in the toy example here below. Now, I am looking for the confidence intervals for the difference between group means. In my case, the point estimate of the mean difference is 64.4. I am looking at the 95% confidence intervals around this point estimate (64.4). Many thanks for your response. library(bootES) a<-c(523,435,478,567,654) b<-c(423,523,421,467,501) bootES(a) bootES(b) __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] ncol() vs. length() on data.frames
should have added: dim(x)[2L] -> length(x) Am 31.03.20 um 16:21 schrieb Prof. Dr. Matthias Kohl: Dear Ivan, if I enter ncol in the console, I get function (x) dim(x)[2L] indicating that function dim is called. Function dim has a method for data.frame; see methods("dim"). The dim-method for data.frame is dim.data.frame function (x) c(.row_names_info(x, 2L), length(x)) Hence, it calls length on the provided data.frame. In addition, some "magic" with .row_names_info is performed, where base:::.row_names_info function (x, type = 1L) .Internal(shortRowNames(x, type)) Best Matthias Am 31.03.20 um 16:10 schrieb Ivan Calandra: Thanks Ivan for the answer. So it confirms my first thought that these two functions are equivalent when applied to a "simple" data.frame. The reason I was asking is because I have gotten used to use length() in my scripts. It works perfectly and I understand it easily. But to be honest, ncol() is more intuitive to most users (especially the novice) so I was thinking about switching to using this function instead (all my data.frames are created from read.csv() or similar functions so there should not be any issue). But before doing that, I want to be sure that it is not going to create unexpected results. Thank you, Ivan -- Dr. Ivan Calandra TraCEr, laboratory for Traceology and Controlled Experiments MONREPOS Archaeological Research Centre and Museum for Human Behavioural Evolution Schloss Monrepos 56567 Neuwied, Germany +49 (0) 2631 9772-243 https://www.researchgate.net/profile/Ivan_Calandra On 31/03/2020 16:00, Ivan Krylov wrote: On Tue, 31 Mar 2020 14:47:54 +0200 Ivan Calandra wrote: On a simple data.frame (i.e. each element is a vector), ncol() and length() will give the same result. Are they just equivalent on such objects, or are they differences in some cases? I am not aware of any exceptions to ncol(dataframe)==length(dataframe) (in fact, ncol(x) is dim(x)[2L] and ?dim says that dim(dataframe) returns c(length(attr(dataframe, 'row.names')), length(dataframe))), but watch out for AsIs columns which can have columns of their own: x <- data.frame(I(volcano)) dim(x) # [1] 87 1 length(x) # [1] 1 dim(x[,1]) # [1] 87 61 __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] ncol() vs. length() on data.frames
Dear Ivan, if I enter ncol in the console, I get function (x) dim(x)[2L] indicating that function dim is called. Function dim has a method for data.frame; see methods("dim"). The dim-method for data.frame is dim.data.frame function (x) c(.row_names_info(x, 2L), length(x)) Hence, it calls length on the provided data.frame. In addition, some "magic" with .row_names_info is performed, where base:::.row_names_info function (x, type = 1L) .Internal(shortRowNames(x, type)) Best Matthias Am 31.03.20 um 16:10 schrieb Ivan Calandra: Thanks Ivan for the answer. So it confirms my first thought that these two functions are equivalent when applied to a "simple" data.frame. The reason I was asking is because I have gotten used to use length() in my scripts. It works perfectly and I understand it easily. But to be honest, ncol() is more intuitive to most users (especially the novice) so I was thinking about switching to using this function instead (all my data.frames are created from read.csv() or similar functions so there should not be any issue). But before doing that, I want to be sure that it is not going to create unexpected results. Thank you, Ivan -- Dr. Ivan Calandra TraCEr, laboratory for Traceology and Controlled Experiments MONREPOS Archaeological Research Centre and Museum for Human Behavioural Evolution Schloss Monrepos 56567 Neuwied, Germany +49 (0) 2631 9772-243 https://www.researchgate.net/profile/Ivan_Calandra On 31/03/2020 16:00, Ivan Krylov wrote: On Tue, 31 Mar 2020 14:47:54 +0200 Ivan Calandra wrote: On a simple data.frame (i.e. each element is a vector), ncol() and length() will give the same result. Are they just equivalent on such objects, or are they differences in some cases? I am not aware of any exceptions to ncol(dataframe)==length(dataframe) (in fact, ncol(x) is dim(x)[2L] and ?dim says that dim(dataframe) returns c(length(attr(dataframe, 'row.names')), length(dataframe))), but watch out for AsIs columns which can have columns of their own: x <- data.frame(I(volcano)) dim(x) # [1] 87 1 length(x) # [1] 1 dim(x[,1]) # [1] 87 61 __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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 rum Multiple ANOVA and Multiple T-test between the same groups?
Have a look at Bioconductor package genefilter, especially functions colttests and colFtests. Best Matthias Am 10.02.19 um 10:35 schrieb AbouEl-Makarim Aboueissa: Dear All: good morning *Re:* How to rum Multiple ANOVA and Multiple T-test between the same groups. Your help will be highly appreciated. *1.* is there a way to run multiple t-tests on different variables between the same two groups. *Data for t-tests:* The data frame “dataTtest” has 5 variables (x1,x2,x3,x4,x5) and one factor (factor1) with 2 levels (group1, group2). x1<-rnorm(20,1,1) x2<-rnorm(20,2,1) x3<-rnorm(20,3,1) x4<-rnorm(20,4,1) x5<-rnorm(20,5,1) factor1<-rep(c("group1", "group2"), each = 10) dataTtest<-data.frame(x1,x2,x3,x4,x5,factor1) dataTtest *2.* is there a way to run *multiple ANOVA* and multiple comparisons *Tukey tests* on different variables between the same groups. *Data for ANOVA tests:* The data frame “dataANOVA” has 6 variables (x1,x2,x3,x4,x5,x6) and one factor (factor2) with 5 levels (group1, group2, group3, group4, group5). x1<-rnorm(40,1,1) x2<-rnorm(40,2,1) x3<-rnorm(40,3,1) x4<-rnorm(40,4,1) x5<-rnorm(40,5,1) x6<-rnorm(40,6,1) factor2<-rep(c("group1", "group2", "group3", "group4", "group5"), each = 8) dataANOVA<-data.frame(x1,x2,x3,x4,x5,x6,factor2) dataANOVA with many thanks abou __ *AbouEl-Makarim Aboueissa, PhD* *Professor, Statistics and Data Science* *Graduate Coordinator* *Department of Mathematics and Statistics* *University of Southern Maine* [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] Paired sample t-test with mi.t.test
Dear Joel, are you trying to apply function mi.t.test from my package MKmisc? Could you please try: mi.t.test(implist, "pre_test", "post_test", alternative = "greater", paired = TRUE, var.equal = TRUE, conf.level = 0.95) x and y are the names of the variables, not the variables themselves. Best Matthias Am 06.04.2017 um 18:32 schrieb Joel Gagnon: Dear all, It is my first time posting on this list so forgive me for any rookie mistakes I could make. I want to conduct t-tests on a dataset that has been imputed using the mice package: imput_pps <- mice(pps, m=20, maxit=20, meth='pmm') # pps is my dataset. It contains items from an 11-item questionnaire gather at pre and post test. So the data set has 22 columns. I then proceed to compute the total scores for the pre and post test on my imputed datasets: long_pps <- complete(imput_pps, action ="long", include = TRUE) long_pps$pre_test <- rowSums(long_pps[ ,c(3:13)]) long_pps$post_test <- rowSums(long_pps[ , c(14:24)]) I then used as.mids to convert back to mids object: mids_pps <- as.mids(long_pps) Next, I created an imputation list object using mitools: implist <- lapply(seq(mids_pps$m), function(im) complete(mids_pps, im)) implist <- imputationList(implist) Now, I want to conduct t-tests using the mi.t.test package. I tried the following code: mi.t.test(implist, implist$pre_test, implist$post_test, alternative = "greater", paired = TRUE, var.equal = TRUE, conf.level = 0.95) When I run this code, R tells me that Y is missing. I know this may sound stupid, but I thought that I specified Y with this line: implist$pre_test, implist$post_test - with implist$pre_test being X and implist$post_test being Y - like I usually do for a normal t-test using the t.test function. It seems I don't quite understand what the Y variable is supposed to represent. Could someone help me figure out what I am doing wrong? You help would be very much appreciated. Best regards, Joel Gagnon, Ph.D(c), Department of Psychology, Université du Québec à Trois-Rivières Québec, Canada [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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 install pkg when "not found" ?
use quotes! install.packages("ggplot2") Am 30.12.2015 um 06:41 schrieb Judson: Using "Install Packages" from CRAN, in RStudio on Windows 7, I downloaded (and supposedly installed) ggplot2 package to here: C:\Program Files\R\R-3.1.0\library\ggplot2_2.0.0\ggplot2 when I try this: require(ggplot2) I get the following: Loading required package: ggplot2 Warning message: In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, : there is no package called 'ggplot2' I also tried this: install.packages(ggplot2) Error in install.packages : object 'ggplot2' not found .. The above library has the usual things like MASS, stats, base, utils, graphics and so on . and if I do this: require(MASS) I get this: Loading required package: MASS so that works . Any idea what I'm doing wrong with ggplot2? ... judson blake [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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] C.D.F
Dear Diba, you could try package distr; eg. library(distr) G1 - Gammad(scale = 0.7, shape = 0.5) G2 - Gammad(scale = 2.1, shape = 1.7) G3 - G1+G2 # convolution G3 For the convolution exact formulas are applied if available, otherwise we use FFT; see also http://www.jstatsoft.org/v59/i04/ (will appear soon) resp. a previous version at http://arxiv.org/abs/1006.0764 hth Matthias Am 11.08.2014 um 10:17 schrieb pari hesabi: Hello everybody, Can anybody help me to write a program for the CDF of sum of two independent gamma random variables ( covolution of two gamma distributions) with different amounts of parameters( the shape parameters are the same)? Thank you Diba [[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. -- Prof. Dr. Matthias Kohl www.stamats.de __ 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] legendre quadrature
you could use package distrEx: library(distrEx) GLIntegrate(function(x) x^2, lower = -1, upper = 1, order = 50) hth Matthias On 01.05.2014 09:43, pari hesabi wrote: Hello everybody I need to approximate the amount of integral by using legendre quadrature. I have written a program which doesn't give me a logical answer; Can anybody help me and send the correct program? For example the approximated amount of integral of ( x ^2) on (-1,1) based on legendre quad rule. integrand-function(x) {x^2} rules - legendre.quadrature.rules( 50 ) order.rule - rules[[50]] chebyshev.c.quadrature(integrand, order.rule, lower = -1, upper = 1) Thank you Diba [[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. -- Prof. Dr. Matthias Kohl www.stamats.de __ 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] joint distribution
have a look at our package distr: library(distr) x1 - Norm(mean = 0, sd = 1) x2 - Binom(size = 1, prob = 0.75) x3 - x1 + x2 plot(x3) # to get density, cdf, quantile, and random numbers use d(x3)(5) p(x3)(0) q(x3)(0.975) r(x3)(20) # you can also have additonal coefficients; eg. x4 - 3*x1 + 2*x2 plot(x4) For more details on the computations see http://arxiv.org/abs/1006.0764 hth, Matthias On 22.04.2014 13:11, Ms khulood aljehani wrote: Hello i have two independent variablesx1 from normal (0,1)x2 from bernoulli (o.75) i need the density estimation of(b1*x1) + (b2*x2) where b1 and b2 are two fixed coefficients thank you [[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. -- Prof. Dr. Matthias Kohl www.stamats.de __ 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] Creating new instances from original ones
see function SMOTE in package DMwR hth Matthias On 31.03.2013 10:46, Nicolás Sánchez wrote: I have a question about data mining. I have a dataset of 70 instances with 14 features that belong to 4 classes. As the number of each class is not enough to obtain a good accuracy using some classifiers( svm, rna, knn) I need to oversampling the number of instances of each class. I have heard that there is a method to do this. It consists in generating these new instances as follows: new_instance original_instance + u(epsilon) U(epsilon) is a uniform number in the range [-epsilon,epsilon] and this number is applied to each feature of the dataset to obtain a new instance without modified the original class. Anybody has used this method to oversampling his data? Anybody has more information about it? 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] integrate function
use pmin instead of min. hth Matthias On 12.02.2013 16:41, dan wang wrote: Hi All, Can any one help to explain why min and max function couldn't work in the integrate function directly. For example, if issue following into R: integrand - function(x) {min(1-x, x^2)} integrate(integrand, lower = 0, upper = 1) it will return this: Error in integrate(integrand, lower = 0, upper = 1) : evaluation of function gave a result of wrong length However, as min(U,V) = (U+V)/2-abs(U-V)/2 Below code working; integrand - function(x){(1-x+x^2)/2-0.5*abs(1-x-x^2)} integrate(integrand, lower = 0, upper = 1) 0.151639 with absolute error 0.00011 I am confused... Dan [[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. -- Prof. Dr. Matthias Kohl www.stamats.de __ 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] Convolve 2 uniform distributed variables
take a look at the distr package library(distr) U - Unif() U2 - U + U # or U2 - convpow(U, 2) plot(U2) # or curve(d(U2)(x), from = 0, to = 2) Best Matthias On 24.05.2012 08:29, c...@mail.tu-berlin.de wrote: Hi, I want to convolve 2 uniform distributed [0,1] variables, so that I get this graphic: http://de.wikipedia.org/wiki/Benutzer:Alfred_Heiligenbrunner/Liste_von_Verteilungsdichten_der_Summe_gleichverteilter_Zufallsvariabler (second graphic) if I do u1-seq(0,1,length=100) fu1=dunif(u1,min=0,max=1) plot(u1,fu1,type=l,xlim=c(-2,2)) u2-seq(0,1,length=100) fu2=dunif(u2,min=0,max=1) u1u2-convolve(u1,rev(u1),typ=o) plot(u1u2) it does not work, could you help me please? The point is the convolve function I think, what do I have to type, to get the correct convolution of two uniform distributed [0,1] variables as to be seen in the second graphic in the given link? Thanks a lot __ 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ 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] Compare String Similarity
you should also look at Bioconductor Package Biostrings hth Matthias Am 19.04.2012 18:01, schrieb R. Michael Weylandt: Though if you do decide to use Levenstein, it's implemented here in R: http://finzi.psych.upenn.edu/R/library/RecordLinkage/html/strcmp.html I'm pretty sure this is in C code so it should be mighty fast. Michael On Thu, Apr 19, 2012 at 11:40 AM, Bert Guntergunter.ber...@gene.com wrote: Wrong list.This is R, not statistics (or linguistics?).Please post elsewhere. -- Bert On Thu, Apr 19, 2012 at 8:05 AM, Alekseiy Beloshitskiy abeloshits...@velti.com wrote: Dear All, I need to estimate the level of similarity of two strings. For example: string1- c(depending,audience,research, school); string2- c(audience,push,drama,button,depending); The words in string may occur in different order though. What function would you recommend to use to estimate similarity (e.g., levenstein, distance)? Appreciate for any advices. -Alex [[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. __ 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ 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] Sweave UFT8 problem
try: Sweave(sim_pi.Rnw, encoding = utf8) Best, Matthias On 15.04.2012 11:41, Philippe Grosjean wrote: Hello, Have you tried to put that command in a comment: %\usepackage[utf8]{inputenc} I haven't tested it in this particular case, but it works in some other situations. Best, Philippe ..¡})) ) ) ) ) ) ( ( ( ( ( Prof. Philippe Grosjean ) ) ) ) ) ( ( ( ( ( Numerical Ecology of Aquatic Systems ) ) ) ) ) Mons University, Belgium ( ( ( ( ( .. On 14/04/12 22:37, Mark Heckmann wrote: Hi, I work on MacOS, trying to Sweave an UFT8 document. AFAI remember R 2.14 used to render a warning when the encoding was not declared when using Sweave. With R 2.15 it seems to render an error. Sweave(sim_pi.Rnw) Error: 'sim_pi.Rnw' is not ASCII and does not declare an encoding Declaring an encoding by adding a line like \usepackage[utf8]{inputenc} in the preamble does the job. In my case though the .Rnw document does no have a preamble as it is just one chapter. All chapters are Sweaved separately (due to computation time). Hence I cannot inject the above line as LaTex will cause an error afterwards. (usepackage{} is only allowed in the preamble which only appears once in the main document, not in each chapter). How can I get around this not using the terminal for Sweaving, like e.g. R CMD Sweave --encoding=utf-8 sim_pi.Rnw ? Thanks Mark [[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. -- Prof. Dr. Matthias Kohl www.stamats.de __ 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] calculate quantiles of a custom function
Dear Gerhard, you could also use package distr; e.g. library(distr) ## use generating function AbscontDistribution D - AbscontDistribution(d = function(x) dbeta(x, 2, 6) + dbeta(x,6,2), low = 0, up = 1, withStand = TRUE) ## quantiles q(D)(seq(0,1,0.1)) Best Matthias On 03.01.2012 19:33, Albyn Jones wrote: What do quantiles mean here? If you have a mixture density, say myf - function(x,p0) p0*dbeta(x,2,6) + (1-p0)*dbeta(x,6,2) then I know what quantiles mean. To find the Pth quantile use uniroot to solve for the x such that myf(x,p0) - P =0. albyn Quoting VictorDelgado victor.m...@fjp.mg.gov.br: Gerhard wrote Suppose I create a custom function, consisting of two beta-distributions: myfunction - function(x) { dbeta(x,2,6) + dbeta(x,6,2) } How can I calculate the quantiles of myfunction? Thank you in advance, Gerhard Gehard, if do you want to know the quantiles of the new distribution created by myfunction. Maybe you can also do: x - seq(0,1,.01) # insert your 'x' q - myfunction(x) # And: quantile(x) 0% 25% 50% 75% 100% 0.00 1.476177 2.045389 2.581226 2.817425 # This gives the sample quantiles. You can also look foward to simulations (like Bert Gunter had suggested) to know better the properties of distributions quantiles obtained after 'myfunction'. - Victor Delgado cedeplar.ufmg.br P.H.D. student www.fjp.mg.gov.br reseacher -- View this message in context: http://r.789695.n4.nabble.com/calculate-quantiles-of-a-custom-function-tp4256887p4257551.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. __ 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. -- Prof. Dr. Matthias Kohl www.stamats.de __ 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] Defining functions inside loops
. __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Package distr and define your own distribution
Dear Gabriel, it is possible. You can define a new class or just an object of an already existing class. Did you look at: library(distr) vignette(newDistributions) as well as ?AbscontDistribution ?DiscreteDistribution Please let me know if you have further questions! Best Matthias On 11.02.2011 16:05, gabriel.ca...@ubs.com wrote: Hi all I am using the Package distr (and related) Do you know if it is possible to define your own distribution (object) GIVEN that you have an analytical form of the probability density function (pdf) ? I would then like to use the standard feature of the distr and related packages. Best regards Giuseppe Gabriel Cardi __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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 generate a random data from a empirical distribition
Hi Dennis, you should take a look at the CRAN task view for distributions http://cran.r-project.org/web/views/Distributions.html Beside that our distr-family of packages might be useful, see also http://www.jstatsoft.org/v35/i10/ http://cran.r-project.org/web/packages/distrDoc/vignettes/distr.pdf Best, Matthias On 27.07.2010 10:37, Dennis Murphy wrote: Hi: On Mon, Jul 26, 2010 at 11:36 AM, xin weixin...@stat.psu.edu wrote: hi, this is more a statistical question than a R question. but I do want to know how to implement this in R. I have 10,000 data points. Is there any way to generate a empirical probablity distribution from it (the problem is that I do not know what exactly this distribution follows, normal, beta?). My ultimate goal is to generate addition 20,000 data point from this empirical distribution created from the existing 10,000 data points. thank you all in advance. The problem, it seems to me, is the leap of faith you're taking that the empirical distribution of your manifest sample will serve as a useful data-generating mechanism for the 20,000 future observations you want to take. I would think that, if you intend to take a sample of 20,000 from ANY distribution, you would want some confidence in the specification of said distribution. Even if you don't know exactly what type of population distribution you're dealing with, there are ways to narrow down the set of possibilities. What is the domain/support of the distribution? For example, the Normal is defined on all of R (as in the real numbers, not our favorite statistical programming language), whereas the lognormal, Gamma and Weibull distributions, among others, are defined on the nonnegative reals. The beta distribution is defined on [0, 1]. Therefore, knowledge of the domain is useful in and of itself. Is it plausible that the distribution is symmetric, or should it have a distinct left or right skew? (Similar comments apply to discrete distributions.) Is censoring or truncation a relevant concern? If there is a random process that well describes how the data you observe are generated, that will certainly narrow down the class of potential data-generating mechanisms/distributions. Once you've narrowed down the class of possible distributions as much as possible, you could look into the fitdistr() function in MASS or the fitdistrplus package on CRAN to test out which candidates seem plausible wrt your existing sample and which are not. You are not likely to be able to narrow it down to one family of distributions, but you should have a much better idea about the characteristics of the distribution that gave rise to your sample of 10,000 (assuming, of course, that it is a *random* sample) after going through this exercise, which you can apply to the generation of the next 20,000 observations. OTOH, if your existing 10,000 observations were not produced by some random process, all bets are off. HTH, Dennis -- View this message in context: http://r.789695.n4.nabble.com/how-to-generate-a-random-data-from-a-empirical-distribition-tp2302716p2302716.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. [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] Robust Asymptotic Statistics (RobASt)
Dear Alex, after installation of ROptRegTS there should be a folder scripts in the package directory of ROptRegTS which includes some further examples. If you are assuming normal errors you should switch to package RobRex which is optimized for normal errors. There you will also find a folder scripts in the package directory. In general, multivariate regression should work. But, I have to admit that we didn't try many multivariate examples so far. Beside that, the user interface is still rudimentary at the moment and some ideas which will clearly reduce computation time (as in the case of package RobLox) are unfortunately not yet implemented. So, if you are willing to accept these drawbacks I will be glad to give you a hand to get our estimators to work for your data. Best, Matthias On 06.06.2010 06:47, Alex Weslowski wrote: Hi all, Other than: http://www.stamats.de/F2006.r Are there other good simple examples out there of using the ROptRegTS package (part of the RobASt project)? I'm hoping to plug it in for multivariate regression. Or is this not a good idea? Just trying to find out how it compares to rlm, lts, glm, etc. Hopefully this makes sense, I'm new to the world of statistics and R. Thanks! St0x __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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-pkgs] new package RFLPtools
The new package RFLPtools is available on CRAN. RFLPtools provides analysis functions for DNA fragment molecular weights (e.g.\ derived from RFLP-analysis) and nucleotide sequence similarities. It aims mainly at the identification of similar or identical fragment patterns to evaluate the amount of different genotypes gained from environmental samples during diversity studies and at further analysis of similarities of nucleotide sequences derived from pairwise sequence alignments (e.g.\ derived from standalone BLAST). Additionally, functions to check the data quality of molecular fingerprints are available. To identify the organisms represented by the extracted fragment patterns (e.g.\ RFLP genotypes), RFLPtools includes a function to compare samples with fragment patterns stored in a reference dataset, from which the taxonomic affiliations are already known. To identify unknown samples in scientific projects, DNA sequences are used, and the gained DNA sequences are compared to existing sequence data via alignments tools, as standalone BLAST. RFLPtools offers tools to generate groups based on tabular report files of sequence comparison of pairwise nucleotide sequence alignment. To get a first impression try: vignette(RFLPtools) as well as library(RFLPtools) example(RFLPplot) example(RFLPrefplot) Best regards, Fabienne Alexandra Matthias -- Dr. Matthias Kohl www.stamats.de ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ 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-pkgs] new version of RobASt-family of packages
The new version 0.7 of our RobASt-family of packages is available on CRAN for several days. As there were many changes we will only sketch the most important ones here. For more details see the corresponding NEWS files (e.g. news(package = RobAStBase) or using function NEWS from package startupmsg i.e. NEWS(RobAStBase)). First of all the new package RobLoxBioC was added to the family which includes S4 classes and methods for preprocessing omics data, in particular gene expression data. ## ## All packages (RandVar, RobAStBase, RobLox, RobLoxBioC, ## ROptEst, ROptEstOld, ROptRegTS, RobRex) ## - TOBEDONE file was added as a starting point for collaborations. The file can be displayed via function TOBEDONE from package startupmsg; e.g. TOBEDONE(distr) - tests/Examples folder for some automatic testing was introduced ## ## Package RandVar ## - mainly fixing of warnings and bugs ## ## Package RobAStBase ## - enhanced plotting, in particular methods for qqplot - unified treatment of NAs - extended implementation for total variation neighbourhoods - implementation of k-step estimator construction extended ## ## Package RobLox ## - na.rm argument added - introduction of finite-sample correction ## ## Package ROptEst ## - optional use of alternative algorithm to obtain Lagrange multipliers using duality based optimization - extended implementation for total variation neighbourhoods - solutions for general parameter transformations with nuisance components - several extensions to the examples in folder scripts - implementation of k-step estimator construction extended ## ## Package ROptEstOld ## - still needed for packages ROptRegTS and RobRex - removed Symmetry and DistributionSymmetry implementation to make ROptEstOld compatible with distr 2.2 ## ## Package ROptRegTS ## - still depends on ROptEstOld ## ## Packages RobRex ## - moved some of the examples in \dontrun{} to reduce check time ... - some minor corrections in ExamplesEstimation.R in folder scripts Best Peter Matthias -- Dr. Matthias Kohl www.stamats.de ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ 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-pkgs] new version of distr-family of packages
The new version 2.2 of our distr-family of packages has now been available on CRAN for several days. As there were many changes we will only sketch the most important ones here. For more details see the corresponding NEWS files (e.g. news(package = distr) or using function NEWS from package startupmsg i.e. NEWS(distr)). First of all a new package distrEllipse was added to the family which includes S4 classes and methods for elliptical symmetric distributions. ## ## All packages (distr, distrEx, distrSim, distrTEst, distrTeach, ## distrMod, distrEllipse, distrDoc as well as ## startupmsg and SweaveListingUtils) ## - TOBEDONE file was added as a starting point for collaborations. The file can be displayed via function TOBEDONE from package startupmsg; e.g. TOBEDONE(distr) - tests/Examples folder for some automatic testing was introduced ## ## Package distr ## - new diagnostic function (more specifically S4 method) qqplot() for quantile-quantile plots to work with distribution objects as argument; in addition to function qqplot() from package stats has functionality for point-wise and simultaneous confidence bands - under the hood: a new slot symmetry was introduced for distributions ## ## Package distrEx ## - enhanced E (expectation), m1df, m2df (clipped 1. and 2. moments) methods --- now you can specify lower and upper integration bounds in the E(.), var(.) functions - new class GPareto for generalized Pareto distribution (implemented by Nataliya Horbenko) ## ## Package distrMod ## - methods for qqplot() for parametric model objects as argument - unified treatment of NAs ## ## Package SweaveListingUtils ## - commands: * SweaveListingPreparations() is more flexible [see help NEWS file] * integrated Andrew Ellis's nice ideas into SweaveListingUtils to use \lstenvironment * individual markup style for base or recommended packages (checked for with new function isBaseOrRecommended()) is distinct now by default (extra color RRecomdcolor) * temporarily changes (like background color) made easier: by new functions + lstdefRstyle to change Rstyle + lstsetRall to change all R-like styles (i.e. Rstyle, Rinstyle, Routstyle, Rcodestyle) * now can specify markup styles for Sinput, Soutput, Scode, separately as Rin, Rout, Rcode * colors now have suffix color, i.e. Rcomment - Rcommentcolor, Rout - Routcolor - vignette: * included an example with escape sequences in vignette * included an example with framed code in vignette Best Peter Matthias -- Dr. Matthias Kohl www.stamats.de ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ 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] COURSE: Introduction to Metabolomics and biomarker research
Introduction to Metabolomics and biomarker research. The course will take place in Innsbruck, 16th to 17th November 2009, and will be held in German. For more details see: http://www.gdch.de/vas/fortbildung/kurse/fortbildung2009.htm#1175 -- Dr. Matthias Kohl www.stamats.de __ 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] help with median for each row
if your matrix has many rows you might want to consider rowMedians from Bioconductor package Biobase. hth, Matthias Greg Hirson schrieb: Edward, See ?apply x = matrix(rnorm(100), ncol = 10) apply(x, 1, median) Hope this helps, Greg Edward Chen wrote: Hi, I tried looking through google search on whether there's a way to computer the median for each row of a nxn matrix and return the medians for each row for further computation. And also if the number of columns in the matrix are even, how could I specify which median to use? Thank you very much! -- Dr. Matthias Kohl www.stamats.de __ 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 change the quantile method in bwplot
for teaching purposes I wrote a corresponding function; cf. qbxp.stats (as well as qboxplot ...) in package MKmisc. hth, Matthias Deepayan Sarkar schrieb: On Tue, Jul 21, 2009 at 7:47 AM, Jun Shenjun.shen...@gmail.com wrote: Uwe, Thank you for your reply. I am still not very clear about the meanings of the arguments in the stats function. To make it clearer, quantile() uses type=7 as default method. I believe this is the method bwplot() uses to calculate the quantiles. I want to use type=6 method for bwplot(). How do I achieve that? Thanks again. Maybe this will be clearer: bwplot() uses the boxplot.stats() function to compute the quantiles used, which in turn uses fivenum(), which has its own quantile calculation (and does not explicitly use quantile()). There is no easy way to allow for type=6 etc. here. bwplot() allows you to replace boxplot.stats() and provide your own alternative. So what you need to do is: (1) write a function, say, 'my.boxpot.stats', that takes the same arguments as boxplot.stats() and returns a similar result, but using your preferred calculation for the quantiles. There are many ways to do this. (2) plug in this function into the bwplot() call; e.g. bwplot(..., stats = my.boxplot.stats) -Deepayan Jun 2009/7/21 Uwe Ligges lig...@statistik.tu-dortmund.de Jun Shen wrote: Hi, everyone, Since quantile calculation has nine different methods in R, I wonder how I specify a method when calling the bwplot() in lattice. I couldn't find any information in the documentation. Thanks. bwplot() uses the panel function panel.bwplot() which allows to specify a function that calculates the statistics in its argument stats that defaults to boxplot.stats(). Hence you can change that function. Example with some fixed values: bwplot( ~ 1:10, stats = function(x, ...) return(list(stats=1:5, n=10, conf=1, 10, out=integer(0))) ) Uwe Ligges [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] trouble using optim for maximalisation of 2-parameter function
try: # first argument of llik.nor has to be the parameter llik.nor-function(theta,x){mu-theta[1];sig-theta[2];-length(x)*log(sqrt(2*pi))-length(x)*log(sig)-(1/(2*sig**2))*sum((x-mu)**2)} # optim by default does minimization # setting fnscale = -1 one obtains a maximization problem optim(c(1,5), llik.nor, x=x, control = list(fnscale = -1)) hth, Matthias Anjali Sing schrieb: Hello, I am having trouble using optim. I want to maximalise a function to its parameters [kind of like: univariate maximum likelihood estimation, but i wrote the likelihood function myself because of data issues ] When I try to optimize a function for only one parameter there is no problem: llik.expo-function(x,lam){(length(x)*log(lam))-(length(x)*log(1-exp(-1*lam* *cons*)))-lam*sum(x)} cons- data-c(.) expomx-optimize(llik.expo,c(0,100),maximum=TRUE,tol=0.0001,x=data) expomx To optimize to two parameters you can't use optimize, so I tried the following to test my input: llik.nor-function(x,theta){mu-theta[1];sig-theta[2];-length(x)*log(sqrt(2*pi))-length(x)*log(sig)-(1/(2*sig**2))*sum((x-mu)**2)} x-c(-1,4,6,4,2) normx-optim(c(1,20),llik.nor) the output should be close to parameters: mu=3 and sigma=2.366 [This I calculated by hand to compare with the output] but in stead of output I get an error: Error in fn(par, ...) : argument theta is missing, with no default I don't understand why this happened? I hope someone can help me with this for I am still a [R]ookie. Kind regards, Sonko Lady [Anjali] [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] implementing Maximum Likelihood with distrMod when only the PDF is known
Dear Guillaume, retrospectively, I'm not sure if the decision to have spezial initialize methods is the optimal way to go. In distrMod and our other packages on robust statistics we don't introduce special initialize methods, but use so-called generating functions. This approach has the advantage that the default initialize method can be used for programming where I find it very useful. I would say it is the canonical (and recommended) way to first define a function as generic (like mu) and then implement methods for it. Choosing names for which there is already an existing definition for a generic function may also not be what you want in general. By defining new methods you have to respect the definition of the generic function; that is, your method definition has to be with respect to the arguments in the given generic function and also dispatching will be on these arguments (cf. scale in the distrMod example). In defining our classes we decided that at least the slot r has to be filled (of class function) whereas d, p and q are of class OptionalFunction. Hence, there are functions to fill d, p and q for given r but, so far not for filling r, p and q given d. A way to avoid implementing r is given in http://www.stamats.de/distrModExample1.R I also do not fill the slots p and q in this example. This avoids the simulation of a large random sample and speeds up the computation of the MLE. However, this is rather a quick and dirty solution and it would of course be better to have a valid definition for r, d, p and q. Best, Matthias guillaume.martin schrieb: Dear Mathias, That's pretty amazing, thanks a lot ! I'll have to look all this through because I don't easily understand why each part has to be set up, in particular the initialize method part seems crucial and is not easy to intuit. From what I get, the actual name we give to a parameter (my original mu for example) is important in itself, and if we introduce new variable names we have to define a new generic, right? The simplest option then is to re-use an existing variable name that has the same properties/range, right? Another general question: my actual pdf is of the same type but not the exact same as the skew normal. In particular, I don't have a rule for building the slot r (eg the one borrowed from the sn package in your example); is it a problem? isn't it sufficient to give slot d, and then you have automatic methods implemented to get from d() to r() slots etc. is that right? Thanks a lot for your help and time ! Best, Guillaume Matthias Kohl a écrit : Dear Guillaume, thanks for your interest in the distrMod package. Regarding your question I took up your example and put a file under: http://www.stamats.de/distrModExample.R Hope that helps ... Don't hesitate to contact me if you have further questions! Best, Matthias guillaume.martin schrieb: Dear R users and Dear authors of the distr package and sequels I am trying to use the (very nice) package distrMod as I want to implement maximum likelihood (ML) fit of some univariate data for which I have derived a theoretical continuous density (pdf). As it is a parametric density, I guess that I should implement myself a new distribution of class AbscontDistributions (as stated in the pdf on creating new distributions in distr), and then use MLEstimator() from the distrMod package. Is that correct or is there a simpler way to go? Note that I want to use the distr package because it allows me to implement simply the convolution of my theoretical pdf with some noise distribution that I want to model in the data, this is more difficult with fitdistr or mle. It proved rather difficult for me to implement the new class following all the steps provided in the example, so I am asking if someone has an example of code he wrote to implement a parametric distribution from its pdf alone and then used it with MLEstimator(). I am sorry for the post is a bit long but it is a complicate question to me and I am not at all skillful in the handling of such notions as S4 - class, etc.. so I am a bit lost here.. As a simple example, suppose my theoretical pdf is the skew normal distribution (available in package sn): #skew normal pdf (default values = the standard normal N(0,1) fsn-function(x,mu=0,sd=1,d=0) {u = (x-mu)/sd; f = dnorm(u)*pnorm(d*u); return(f/sd)} # d = shape parameter (any real), mu = location (any real), sd = scale (positive real) #to see what it looks like try x-seq(-1,4,length=200);plot(fsn(x,d=3),type=l) #Now I tried to create the classes SkewNorm and SkewNormParameter copying the example for the binomial ##Class:parameters setClass(SkewNormParameter, representation=representation(mu=numeric,sd=numeric,d=numeric), prototype=prototype(mu=0,sd=1,d=0,name=gettext(Parameter of the Skew Normal distribution)), contains=Parameter ) ##Class: distribution (created using the pdf of the skew normal defined
Re: [R] implementing Maximum Likelihood with distrMod when only the PDF is known
Dear Guillaume, thanks for your interest in the distrMod package. Regarding your question I took up your example and put a file under: http://www.stamats.de/distrModExample.R Hope that helps ... Don't hesitate to contact me if you have further questions! Best, Matthias guillaume.martin schrieb: Dear R users and Dear authors of the distr package and sequels I am trying to use the (very nice) package distrMod as I want to implement maximum likelihood (ML) fit of some univariate data for which I have derived a theoretical continuous density (pdf). As it is a parametric density, I guess that I should implement myself a new distribution of class AbscontDistributions (as stated in the pdf on creating new distributions in distr), and then use MLEstimator() from the distrMod package. Is that correct or is there a simpler way to go? Note that I want to use the distr package because it allows me to implement simply the convolution of my theoretical pdf with some noise distribution that I want to model in the data, this is more difficult with fitdistr or mle. It proved rather difficult for me to implement the new class following all the steps provided in the example, so I am asking if someone has an example of code he wrote to implement a parametric distribution from its pdf alone and then used it with MLEstimator(). I am sorry for the post is a bit long but it is a complicate question to me and I am not at all skillful in the handling of such notions as S4 - class, etc.. so I am a bit lost here.. As a simple example, suppose my theoretical pdf is the skew normal distribution (available in package sn): #skew normal pdf (default values = the standard normal N(0,1) fsn-function(x,mu=0,sd=1,d=0) {u = (x-mu)/sd; f = dnorm(u)*pnorm(d*u); return(f/sd)} # d = shape parameter (any real), mu = location (any real), sd = scale (positive real) #to see what it looks like try x-seq(-1,4,length=200);plot(fsn(x,d=3),type=l) #Now I tried to create the classes SkewNorm and SkewNormParameter copying the example for the binomial ##Class:parameters setClass(SkewNormParameter, representation=representation(mu=numeric,sd=numeric,d=numeric), prototype=prototype(mu=0,sd=1,d=0,name=gettext(Parameter of the Skew Normal distribution)), contains=Parameter ) ##Class: distribution (created using the pdf of the skew normal defined above) setClass(SkewNorm,prototype = prototype( d = function(x, log = FALSE){fsn(x, mu=0, sd=1,d=0)}, param = new(SkewNormParameter), .logExact = TRUE,.lowerExact = TRUE), contains = AbscontDistribution ) #so far so good but then with setMethod(mu, SkewNormParameter, function(object) obj...@mu) #I get the following error message: Error in setMethod(mu, SkewNormParameter, function(object) obj...@mu) : no existing definition for function mu I don't understand because to me mu is a parameter not a function... maybe that is too complex programming for me and I should switch to implementing my likelihood by hand with numerical convolutions and optim() etc., but I would like to know how to use distr, so if there is anyone who had the same problem and solved it, I would be very grateful for the hint ! All the best, Guillaume -- Dr. Matthias Kohl www.stamats.de __ 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] Drawing dendrogram
maybe ?dendrapply is helpful - in particular see the example section. There is also some code at http://www.stamats.de/Heatmap.R (with commentary in German) hth, Matthias Seyit Ali Kayis schrieb: Dear all, I would like to draw a dendrogram and mark some parts/branches (by using segments) including their labels. If I draw it without specifying the length of x axix, I am able to do that (as in My dendrogram 1 of the following codes). However, if I want to specify the x axix, I am not able to draw marking line (by using segments) including labels (as in My dendrogram 2 of the following codes). Is there a way that I can include my labels into as well? Any help is deeply appreciated. Kind Regards Seyit Ali Codes: set.seed(201) x-matrix(rnorm(100, 0, 1), ncol=20, byrow=TRUE) myclust- hclust(dist(cor(x), method='euclidian'), method=ave) mydend-as.dendrogram(myclust) par(mfrow=c(1,2), mar=c(5,4,4,3), oma=c(1, 3, 1, 3)) plot(mydend, xlim=c(4, -0.2), horiz = TRUE, main=My Dendrogram 1) x0- -0.35 y0- 0.5 x1- 1.8 y1- 0.5 segments(x0, y0, x1, y1, col=red, lty=2) x0- -0.35 y0- 4.4 x1- 1.8 y1- 4.4 segments(x0, y0, x1, y1, col=red, lty=2) x0- -0.35 y0- 0.5 x1- -0.35 y1- 4.4 segments(x0, y0, x1, y1, col=red, lty=2) x0- 1.8 y0- 0.5 x1- 1.8 y1- 4.4 segments(x0, y0, x1, y1, col=red, lty=2) plot(mydend, xlim=c(4, 0.5), horiz = TRUE, main=My Dendrogram 2, dLeaf=0.35)#, yaxt=n) x0- 0.3 y0- 0.5 x1- 1.8 y1- 0.5 segments(x0, y0, x1, y1, col=red, lty=2) x0- 0.3 y0- 4.4 x1- 1.8 y1- 4.4 segments(x0, y0, x1, y1, col=red, lty=2) x0- 0.3 y0- 0.5 x1- 0.3 y1- 4.4 segments(x0, y0, x1, y1, col=red, lty=2) x0- 1.8 y0- 0.5 x1- 1.8 y1- 4.4 segments(x0, y0, x1, y1, col=red, lty=2) -- Dr. Seyit Ali KAYIS Selcuk University Faculty of Agriculture Kampus, Konya, TURKEY s_a_ka...@yahoo.com,s_a_ka...@hotmail.com Tell: +90 332 223 2830 Mobile: +90 535 587 1139 Fax: +90 332 241 0108 Greetings from Konya, TURKEY http://www.ziraat.selcuk.edu.tr/skayis/ -- _ [[elided Hotmail spam]] [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] Heatmap: draw horizontal line
Is heatmap.2(mat,Rowv=NA,trace = 'none',key=F, add.expr = abline(h = c(3.5,5.5), lwd = 3)) what your are looking for? hth, Matthias Paul Evans schrieb: Hi, I wanted to draw a heatmap with some horizontal lines. For example: #-- code -- mat - matrix(-1:1,7,9) heatmap.2(mat,Rowv=NA,trace = 'none',key=F) #end code - In this heatmap, I want to subgroup the rows. For instance, I would like to group rows 5,6 7 together, and would draw a black line between rows 4 5 to demarcate the group. Similarly, I want to group rows 3 4 in one group, and would like to put a black line between rows 2 3. I guess I'm looking for the 'hline' parameter or something, but couldn't find it in the documentation. Any help would be appreciated! [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] any other fast method for median calculation
there is function rowMedians in Bioconductor package Biobase which works for numeric matrices and might help. Matthias Dimitris Rizopoulos wrote: S Ellison wrote: Sorting with an appropriate algorithm is nlog(n), so it's very hard to get the 'exact' median any faster. However, if you can cope with a less precise median, you could use a binary search between max(x) and min(x) with low tolerance or comparatively few iterations. In native R, though, that isn;t going to be fast; interpreter overhead will likely more than wipe out any reduction in number of comparisons. In any case, it looks like you are not constrained by the median algorithm, but by the number of calls. You might do a lot better with apply, though apply(df,2,median) well, for data frames, I think sapply(...) or even unlist(lapply(...)) will be faster, e.g., mat - matrix(rnorm(50*2e05), 50, 2e05) DF - as.data.frame(mat) invisible({gc(); gc()}) system.time(apply(DF, 2, median)) invisible({gc(); gc()}) system.time(sapply(DF, median)) invisible({gc(); gc()}) system.time(unlist(lapply(DF, median), use.names = FALSE)) Best, Dimitris On my system 200k columns were processed in negligible time by apply and I'm still waiting for mapply. S Zheng, Xin (NIH) [C] zheng...@mail.nih.gov 14/04/2009 05:29:40 Hi there, I got a data frame with more than 200k columns. How could I get median of each column fast? mapply is the fastest function I know for that, it's not yet satisfied though. It seems function median in R calculates median by sort and mean. I am wondering if there is another function with better algorithm. Any hint? Thanks, Xin Zheng __ 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. *** This email and any attachments are confidential. Any use...{{dropped:8}} __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Finding Lambda in Poisson distribution
Hi, assuming the assumption of Poisson distribution is true (or at least approximately true), you can use MLE as mentioned by Roland using (fitdistr of MASS, mle of stats4, MLEstimator of distrMod, probably more). Using package distrMod you can compute minimum distance estimators (cf. function MDEstimator). Using package ROptEst you can compute (in some sense) optimally robust estimators (cf. function roptest). Best, Matthias Rau, Roland wrote: Hi, -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Saeed Ahmadi Sent: Monday, March 02, 2009 3:16 PM To: r-help@r-project.org Subject: [R] Finding Lambda in Poisson distribution Hi, I have a dataset. First of all, I know that my dataset shall follow the Poission distribution. Now I have two questions: 1) How can I check that my data follow the Poisson distribution? 2) How can I calculate Lambda of my data? is this maybe some homework? For 2): I simulated data and did some simple MLE. I don't know if this is the recommended way, but it worked fine for me. Best, Roland -- This mail has been sent through the MPI for Demographic Research. Should you receive a mail that is apparently from a MPI user without this text displayed, then the address has most likely been faked. If you are uncertain about the validity of this message, please check the mail header or ask your system administrator for assistance. __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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 normalize to a set of internal references
you should better ask this question on the Bioconductor mailing list. For qPCR normalisation strategies take a look at http://www.gene-quantification.info/ Best, Matthias Waverley wrote: Thanks for the advice. My question is more on how to do this? Let me use a biology gene analysis example to illustrate: In biology, there are always some house keeping genes which differ little even at pathological conditions. We know that at different batches, there are external factors affect the measurements. For example, overall signal intensity might be different due to lab reagents. A simplified picture: Day 1: Using control samples, I have measured #1 to #110 genes and get data. Day 2: Using disease samples, I have measured again #1 to #110 genes and get data. For those two data sets, I noticed the overall signal intensity in Day 1, for each gene, is more than Day 2. I know, from biological literature, gene 101 to 110, are house keeping genes, should not change much between disease and control. My questions arise, technically, how do I use gene 101 to 110 values to adjust the signals of gene 1 to 100 such that the batch effect can be corrected. The differences revealing from the comparative analysis of 1 ~ 100 genes between disease and control are due to biology rather than lab artifacts. So the question is how to do that mathematically? If I have only one house keeping gene, then I can divide every gene to that to normalize, then compare. But now I have 10 genes which can be utilized for normalization. I assume, the more reference genes to be used, the better, under this context. Can you help again? Thanks much in advance. Waverley wrote: Hi, I have a question of the method as how to normalize the data sets according to a set of the internal measurements. For example, I have performed two batches of experiments contrasting two different conditions (positive versus negative conditions): one at a time. 1. each experiment, I measure signals of variable v1 to v100. I want to understand v1 to v100 change under these two contrasting conditions 2. Also I know different variables v101 to v1110, total of 10 of them, although they are different from each other, but they would of the same or similar values under these two contrasting conditions 3. How do I do the internal normalization? How can I use the the variable v101 to v110 values to normalize the measures of v1 to v100 at either positive or negative condition to minimize batch effect? I hope the comparisons of values (v1 to v100) between two different conditions can be more accurate and robust to external noises. In general, I have a couple of matrices of the same dimensions and a reference matrix of values to be used as reference values to be normalize to. How should I do that? I don't understand your problem well, but in general internal normalization is by and large an attempt to avoid appropriate modeling (e.g., incorporating block effects or certain covariates in a regression model), and results in overstated confidence of the final estimates by not taking into account the imprecision in the normalizing factors. Frank -- Dr. Matthias Kohl www.stamats.de __ 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] Does The Randvar package contain a virus(Malware) ?
there were only minor changes in the latest version of RandVar (from 0.6.6 to 0.6.7). Might this be a mirror problem? Best Matthias Timthy Chang wrote: Today,I update the packages in R. but AntiVir Guard dectects the Randvar package as affected file. What happen ?? Thank you for your answer. http://www.nabble.com/file/p22281513/qq.gif -- Dr. Matthias Kohl www.stamats.de __ 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] heatmap.2 color issue
the heatmapCol function of package MKmisc might help ... try: library(MKmisc) example(heatmapCol) Best Matthias Liu, Hao [CNTUS] wrote: Dear All: I tried to use heatmap.2 to generate hierarchical clustering using the following command: heatmap.2(datamatrix, scale=row, trace=none, col=greenred(256), labRow=genelist[,1], margins=c(10,10), Rowv=TRUE, Colv=TRUE) datamatrix is subset of a RMA normalized data subset by a genelist. The problem is a lot of times, the z-score in key are from, like -5 to 15 or -15 to 5, as a result, the zero of z distribution are are either green region or red region of the key, the resulting heatmap are either generally greenish or redish. I wonder if there is a way to make the heatmap more balanced between red and green, I tried to read the heatmap.2 help but could not get a clear idea. Thanks Hao [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] Drawing from an empirical distribution
one could also use package distr; e.g., library(distr) x - 1:10 D - DiscreteDistribution(x) ## = r, d, p and q functions (also with log-argument) r(D)(5) p(D)(4) d(D)(1) q(D)(0.3) Best, Matthias roger koenker wrote: Sure, but it would be more 'fun' to modify ecdf() slightly to produce an ecqf() function -- essentially reversing the arguments to approxfun()-- and then use ecqf(runif(whatever)) no nit-picking about efficiency, please. url:www.econ.uiuc.edu/~rogerRoger Koenker emailrkoen...@uiuc.eduDepartment of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Jan 6, 2009, at 4:42 PM, Antonio, Fabio Di Narzo wrote: If the ecdf is 'ecdf(x)', do just: sample(x, size=whatever, replace=TRUE) HTH, Antonio. 2009/1/6 culpritNr1 ig2ar-s...@yahoo.co.uk: Hi All, Does anybody know if there is a simple way to draw numbers from an empirical distribution? I know that I can plot the empirical cumulative distribution function this easy: plot(ecdf(x)) Now I want to pick a number between 0 and 1 and go back to domain of x. Sounds simple to me. Any suggestion? Thank you, Your culprit (everybody needs a culprit) -- View this message in context: http://www.nabble.com/Drawing-from-an-empirical-distribution-tp21320810p21320810.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. -- Antonio, Fabio Di Narzo Ph.D. student at Department of Statistical Sciences University of Bologna, Italy __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] registering a generic method for a class
In case of S4 objects show is called. Hence, you should implement new S4-methods for show ... Matthias m.u.r. wrote: Sorry in advance if this is too simple a question, but I'm stuck with some odd behavior and I can't find the text to rid myself of this (admittedly somewhat trivial) problem. Note that I've done generic programming with S3 objects in the past, but I've never really dabbled in creating S4 objects until now. So, I've created a new S4 class, call it myclass, e.g. : setClass(myclass, representation(slot1 = list)); And I've also created an implementation of print for this class via: setMethod(print, myclass, function(x) { print(paste(a myclass instance with a list of length, length(x...@slot1), in slot1)); }); Now I can create a new object, say via: myobj - new(myclass, slot1 = list()); The question now is the difference in output when I type the object name alone on the command prompt vs when I type print(myobj) on the command prompt, as seen below: myobj An object of class “myclass” Slot slot1: list() print(myobj) [1] a myclass instance with a list of length 0 in slot1 I thought that the print method is called when an object is called from the command line, and to check this I checked the registered implementations of print via: methods(print) and in the long output I didn't see the print.myclass version. I had thought that the setMethod method would essentially create this function for me... but I must be missing a key step somewhere along the line. Can anyone help me figure out what I would need to do to get the generic version of the method (in this toy example, print) to see my class' method as an implementation of the generic? Thanks much for any help! __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Computational Probability
indeed one could use package distr ... library(distr) X - Unif(Min = 0, Max = 1) Y - convpow(X, 10) p(Y)(6) - p(Y)(4) Best Matthias Prof Brian Ripley wrote: Look at packages distr* : they can do your example and might do what your real applications. On Fri, 26 Dec 2008, Rory Winston wrote: Hi Firstly , happy Christmas to R-Help! Secondly, I wonder if anyone can help me with the following query: I am trying to reproduce some explicit probability calculations performed in APPL (a Maple extension for computational probability). For instance, in APPL, to compute the probability that the sum of 10 iid uniform variables [0,1] will be between 4 and 6, (i..e Pr( 4 \sum_{i=1}^{10}X_i 6)), I can type: X := UniformRV(0, 1); Y := ConvolutionIID(X, 10); CDF(Y,6) - CDF(Y,4); which gives the required probability .7222. Is there any way to perform these type of calcuations in R in a general way? I realise that a lot of the machinery behind these computations comes from Maple's symbolic engine, but are there any R extensions for these kind of calculation? Cheers Rory [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] Sweave/Rweave and results=verbatim
it should work without quotes; i.e., results = verbatim hth, Matthias Oliver Bandel wrote: Hello, in a Rnw-file I have this used stuff to try out tex-results... == = texme - function() cat( {\\bf Hallo, das ist voll fett!}\n ) @ results=verbatim= texme() @ == I used this command: R CMD Sweave example.Rnw and got this error: == [...] 7 : echo 8 : echo term verbatim (label=) Error in match.arg(options$results, c(verbatim, tex, hide)) : 'arg' should be one of verbatim, tex, hide Calls: Sweave - SweaveParseOptions - check - match.arg Execution halted == What did I do wrong here? Any hint is welcome, Oliver __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] check if a certain ... argument has been passed on to my user-defined function
and take a look at ?hasArg hth, Matthias Charles C. Berry wrote: See ?match.call and note the expand.dots arg. HTH, Chuck On Thu, 11 Dec 2008, Mark Heckmann wrote: Hi, How can I check if a certain ... argument has been passed on to my user-defined function or not? foo - function(data, ...) { ### here I want to check whether xlab was passed with the ... arguments ### or if the ... arguments did not contain an xlab argument } I tried missing(xlab) , exists(xlab) and several other things but did not find a solution. TIA, Mark __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:cbe...@tajo.ucsd.eduUC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Dr. Matthias Kohl www.stamats.de __ 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] S4 slot containing either aov or NULL
Dear Thomas, take a look at setOldClass ... ## register aov (an S3-class) as a formally defined class setOldClass(aov) ## list and some others are special cases; cf. class ? list ## now your code should work setClassUnion(aovOrNULL, c(aov, NULL)) setClass(c1, representation(value = aovOrNULL)) y1 - new(c1, value = NULL) #trying to assign an aov object to the slot doesn't work utils::data(npk, package=MASS) npk.aov - aov(yield ~ block + N*P*K, npk) y2 - new(c1, value = npk.aov ) Best, Matthias Thomas Kaliwe wrote: Dear listmembers, I would like to define a class with a slot that takes either an object of class aov or NULL. I have been reading S4 Classes in 15 pages more or less and Lecture: S4 classes and methods #First I tried with list and NULL setClass(listOrNULL) setIs(list, listOrNULL) setIs(NULL, listOrNULL) #doesn't work #defining a union class it works with list and null setClassUnion(listOrNULL, c(list, NULL)) setClass(c1, representation(value = listOrNULL)) y1 = new(c1, value = NULL) y2 = new(c1, value = list(a = 10)) #but it won't work with aov or null setClassUnion(aovOrNULL, c(aov, NULL)) setClass(c1, representation(value = aovOrNULL)) y1 = new(c1, value = NULL) #trying to assign an aov object to the slot doesn't work utils::data(npk, package=MASS) npk.aov - aov(yield ~ block + N*P*K, npk) y2 = new(c1, value = npk.aov ) Any ideas? Thank you Thomas Kaliwe __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Mallows' distance or Earth Mover's distance in R?
Hi Rainer, we have not implemented the Mallow's distance so far. Out package distrEx includes implementations of the Hellinger, the Kolmogorov, the total variation and the Cramer von Mises distance. Combined with our package distrMod one can compute minimum distance estimators. library(distrMod) ?MCEstimator ?MDEstimator hth, Matthias Rainer M Krug wrote: Hi I am looking for an implementation (or alternative to) Mallow's distance or the Earth Mover's distance to compare distributions or unnormalized distributions (signatures). Is there an implementation in R or can somebody recommend an alternative? Thanks Rainer -- Dr. Matthias Kohl www.stamats.de __ 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] folded normal distribution in R
Dear Andreas, try: library(distr) X - Norm() ## standard normal distribution Y - abs(X) plot(Y) Best regards, Matthias Andreas Wittmann wrote: 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. -- Dr. Matthias Kohl www.stamats.de __ 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] plot - central limit theorem
Hi Jörg, our package distrTeach has functions to visualize the central limit theorem and the law of large numbers for arbitrary univariate distributions; e.g., library(distrTeach) D - sin(Norm()) + Pois() plot(D) illustrateCLT(D, len = 10) illustrateLLN(D, m = 10) Best Matthias Jörg Groß wrote: Hi, Is there a way to simulate a population with R and pull out m samples, each with n values for calculating m means? I need that kind of data to plot a graphic, demonstrating the central limit theorem and I don't know how to begin. So, perhaps someone can give me some tips and hints how to start and which functions to use. thanks for any help, joerg __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] random normally distributed values within range
one could also use the Truncate-methods of package distr; cf. http://tolstoy.newcastle.edu.au/R/e5/help/08/09/1870.html If the situation is extreme Peter Daalgard gave some sophisticated code at http://tolstoy.newcastle.edu.au/R/e5/help/08/09/1892.html Best, Matthias Rolf Turner wrote: On 7/10/2008, at 11:54 AM, Achaz von Hardenberg wrote: Hi all, I need to create 100 normally distributed random values (X) which can not exceed a specific range (i.e. 0XY). With rnorm I cannot specify Max and min values among which values have to stay, like in runif so does some other simple way exist to do this with normally distributed random values? Presumably you want a truncated normal distribution. Duncan Murdoch posted some neat code on this list, back in the end of July this year, to generate samples from such a distribution. See: http://finzi.psych.upenn.edu/R/Rhelp02a/archive/137154.html cheers, Rolf Turner ## Attention:\ This e-mail message is privileged and confid...{{dropped:9}} __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] help on sampling from the truncated normal/gamma distribution on the far end (probability is very low)
you could use package distr and function Truncate; e.g. library(distr) N - Norm(mean = -4, sd = 1) NT - Truncate(N, lower = 0, upper = Inf) r(NT)(10) Unfortunatelly, your example using sd = 0.1 is very extreme and Truncate doesn't work; see also pnorm(0, mean = -4, sd = 0.1, lower.tail = FALSE) == 0 ## which on my system is TRUE Best, Matthias Moshe Olshansky wrote: Well, I made a mistake - your lambda should be 400 and not 40!!! --- On Thu, 18/9/08, Moshe Olshansky [EMAIL PROTECTED] wrote: From: Moshe Olshansky [EMAIL PROTECTED] Subject: Re: [R] help on sampling from the truncated normal/gamma distribution on the far end (probability is very low) To: r-help@r-project.org, Daniel Davis [EMAIL PROTECTED] Received: Thursday, 18 September, 2008, 5:00 PM Hi Sonia, If I did not make a mistake, the conditional distribution of X given that X 0 is very close to exponential distribution with parameter lambda = 40, so you can sample from this distribution. --- On Mon, 15/9/08, Daniel Davis [EMAIL PROTECTED] wrote: From: Daniel Davis [EMAIL PROTECTED] Subject: [R] help on sampling from the truncated normal/gamma distribution on the far end (probability is very low) To: r-help@r-project.org Received: Monday, 15 September, 2008, 2:28 PM Hi, guys, I am trying to sample from a truncated normal/gamma distribution. But only the far end of the distribution (where the probability is very low) is left. e.g. mu = - 4; sigma = 0.1; The distribution is Normal(mu,sigma^2) truncated on [0,+Inf]; How can I get a sample? I tried to use inverse CDF method, but got Inf as answers. Please help me out. Also, pls help me on the similar situation on gamma dist'n. Thanks, Sonia [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Dr. Matthias Kohl www.stamats.de __ 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] All possible pairs of two variables
take a look at the code of the function pairwise.table. It shows a way using function outer. hth, Matthias Yihui Xie wrote: Just simple loops as follows: z = NULL for (y in 0:10) { for (x in 0:y) { z = rbind(z, c(x, y)) } } z Yihui On Thu, Sep 11, 2008 at 12:34 PM, Ron Michael [EMAIL PROTECTED] wrote: I have two variables (x,y) : x : it takes all integer values from 0 to y and, y : takes all values from 0, 10 I am looking for some R code to find all possible pairs of (x,y). Can anyone please help me? -- Dr. Matthias Kohl www.stamats.de __ 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] fft: from characteristic funtion to density function
take a look at http://finzi.psych.upenn.edu/R/Rhelp02a/archive/130153.html hth, Matthias Jindan Zhou wrote: Hello all! I've posted the question before but I am still struggling for an answer, please help if you can;-) Suppose a discrete series of data is generated by the following equation: CF=exp(-(t^2)/2) which is the characteristic function of a random variable X with standard normal distribution, how do I *numerically* solve for its probability density function? i.e., obtain a series of data that plots a well know bell-shape curve? Or, find x such that Pr(Xx)=0.05? I have no problem to derive the forward and inverse Fourier equation pair, but just no clue on numerical method. A working code on this simplest exercise would clear off many, many of my confusions I have generated CF with t of 256 points equally spaced between -4 to 4 in Excel, but don't know how to proceed. Thanks for your help! Jindan __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Package for Tidying-up R-code
http://cran.r-project.org/doc/manuals/R-exts.html#Tidying-and-profiling-R-code hth, Matthias Don MacQueen wrote: I don't know what kind of tidying you want, but ESS within emacs is a possibility. -Don At 6:32 PM +0900 9/5/08, Gundala Viswanath wrote: Dear all, Is there any such package? I am thinking of something equivalent to Perl::Tidy. For example with VI editor one can visual highlight a portion of Perlcode and then issue the command: !perltidy then the code will be automatically arranged. - Gundala Viswanath Jakarta - Indonesia __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Plotting horizontal dendrograms
Dear Charles, the sample code you cite has an additional as.dendrogram compared to your code. Hence, dend1 is of class dendrogram (class(dend1)) and if you call plot(dend1) the method for dendrograms (i.e., plot.dendrogram) is executed which has a horiz argument; confer ?plot.dendrogram. Whereas in your code h is of class hclust (class(hclust)) and therefore plot(h) leads to a call of plot.hclust which has no horiz argument; cf. ?plot.hclust. So, adding as.dendrogram in your code should remove the warning messages and give what you want; e.g. plot(as.dendrogram(h), horiz = TRUE) Best, Matthias Charles K. Minns wrote: I am using hclust and plot to produce dendrograms. Using my input data I am able to complete an analysis and obtain a vertical plot. I want to be able to plot the dendrogram horizontally.I am using version 2.6 of R and have updated my packages recently. Using the sample script for dendrograms I can produce a horizontal plot using the instruction horiz = TRUE in plot(). When I use the same instruction in my own script I get a warning message that horiz, and horizontal, are not recognized commands in plot. I cannot find any documentation in R on how to change the orientation of plots and general searches for information on this subject. Can someone explain why my script does not work. Sample code follows: --- # Sample code extracted from the Dendrogram documentation in R # require(graphics); require(utils) hc - hclust(dist(USArrests), ave) (dend1 - as.dendrogram(hc)) # print() method str(dend1) # str() method str(dend1, max = 2) # only the first two sub-levels op - par(mfrow= c(2,2), mar = c(5,2,1,4)) plot(dend1, horiz = TRUE) plot(dend1) # # My code # x - read.table(Rclust2.data,header=TRUE) m - data.matrix(x) d - as.dist(m) h - hclust(d) plot(h, horiz = TRUE) -- The three plot commands execute. Plot 1 is horizontal, plot 2 is vertical and plot 3 is vertical. Plus plot 3 generates the console output: - # # My code # x - read.table(Rclust2.data,header=TRUE) m - data.matrix(x) d - as.dist(m) h - hclust(d) plot(h, horiz = TRUE) Warning messages: 1: In plot.hclust(h, horiz = TRUE) : horiz is not a graphical parameter 2: In plot.hclust(h, horiz = TRUE) : horiz is not a graphical parameter 3: In title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...) : horiz is not a graphical parameter --- I would appreciate help solving this problem. Doubtless there is a simple answer. [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] Eliminating [...] from print
Dear Ted, what about cat(\n) or message() Best, Matthias (Ted Harding) wrote: Hi Folks, This must be easy but I've not managed to locate the solution! Basically: I'm using sink() to save successively obtained results, e.g. I construct a set of regression coefficients etc. with names rows and columne (I want to see the names in the output) as an object (say Object), and then I emit it with print(Object) So far so good. But I also want to print a little header and separator prior to each such print(Object) which would identify it and explain what's going on. BUT: I can't get rid of the [1] which comes out each time I print a character string. Getting rid of the quotes is easy: just use quote=FALSE. For example: print (,quote=FALSE) gives [1] But how do I get rid of that [1]?? With thanks, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 01-Jun-08 Time: 18:53:00 -- XFMail -- __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Heatmap.2 - eliminate cluster and dendrogram
Dear Peter, heatmap.2(z, dendrogram = none, Rowv = FALSE, Colv = FALSE) should do what you want. Best, Matthias Peter Scacheri wrote: Using the heatmap.2 function, I am trying to generate a heatmap of a 2 column x 500 row matrix of numeric values. I would like the 1st column of the matrix sorted from the highest to the lowest values - so that the colors reflected in the first column of the heatmap (top to bottom) go from red to green. After sorting the matrix (z), I tried the following command, but the data remains clustered. heatmap.2((z),col=greenred(100),dendrogram=NULL,Rowv=FALSE) I also tried the following, but the data remained clustered. heatmap.2((z),col=greenred(100),Colv=1:ncol(z)) ANY IDEAS?? Thanks so much, Peter __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Function for subset of cases/lines
Dear Stefan, what about q1-c(4660,5621,5629,8030,8080,8180,8501,8190,8370,8200) x - data.frame(number = c(4660, 5010, 5621, 5629, 8030, 8080 , 8090, 8180, 8501, 8190, 8200, 8370, 8200), height = c(2.5, 1.4, 0.8, 2.3, 2.5, 2.4, 0.9, 1.4, 1.2, 1.9, 2.0, 2.1, 1.8)) ind - x[,number] %in% q1 mean(x[ind,height]) ## or shorter mean(x[x[,number] %in% q1,height]) Best, Matthias Stefan Uhmann wrote: Hi, I have a vector: q1-c(4660,5621,5629,8030,8080,8180,8501,8190,8370,8200) The following command gives me the mean of its elements: mean(q1) [1] 7346.1 What can I do to do the same for the variable 'height', but only for the cases/rows which have one of the elements of q1 as 'number': number height 1 46602.5 2 50101.4 3 56210.8 4 56292.3 5 80302.5 6 80802.4 7 80900.9 8 81801.4 9 85011.2 10 81901.9 11 82002.0 12 83702.1 13 82001.8 Thanks for any help, Stefan __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] Heatmap.2 - eliminate cluster and dendrogram
Dear Peter, the call works for me; e.g. x - matrix(rnorm(100), ncol = 4) heatmap.2(x, dendrogram=none, Colv = FALSE, Rowv = FALSE) There seems to be something wrong with your matrix z ... Best, Matthias Peter Scacheri wrote: Thanks Matthias, I tried that, but received the following error message: Error in image.default(1:nc, 1:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 + : dimensions of z are not length(x)(+1) times length(y)(+1) -Peter At 10:18 AM +0200 5/15/08, Matthias Kohl wrote: Dear Peter, heatmap.2(z, dendrogram = none, Rowv = FALSE, Colv = FALSE) should do what you want. Best, Matthias Peter Scacheri wrote: Using the heatmap.2 function, I am trying to generate a heatmap of a 2 column x 500 row matrix of numeric values. I would like the 1st column of the matrix sorted from the highest to the lowest values - so that the colors reflected in the first column of the heatmap (top to bottom) go from red to green. After sorting the matrix (z), I tried the following command, but the data remains clustered. heatmap.2((z),col=greenred(100),dendrogram=NULL,Rowv=FALSE) I also tried the following, but the data remained clustered. heatmap.2((z),col=greenred(100),Colv=1:ncol(z)) ANY IDEAS?? Thanks so much, Peter __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] dataframes to a list
Does mget(paste(data, 1:50, sep = ), envir = globalenv()) give what you want? hth, Matthias Shubha Vishwanath Karanth wrote: How do I proceed from this stage? paste(data,1:10,sep=,collapse=,) Shubha Karanth | Amba Research Ph +91 80 3980 8031 | Mob +91 94 4886 4510 Bangalore * Colombo * London * New York * San José * Singapore * www.ambaresearch.com -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Shubha Vishwanath Karanth Sent: Wednesday, May 14, 2008 3:59 PM To: [EMAIL PROTECTED] Subject: [R] dataframes to a list Hi R, I have the data frames, data1, data2data50. Now I want to put all of these in a single list. But, list(data1, data2,.data50) is very big to write. How do I then do it? Thanks, Shubha This e-mail may contain confidential and/or privileged i...{{dropped:13}} __ 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. This e-mail may contain confidential and/or privileged i...{{dropped:10}} __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] call function immediately before plot.new()
Hi Jake, are you looking for argument panel.first of plot.default? ?plot.default panel.first: an expression to be evaluated after the plot axes are set up but before any plotting takes place. This can be useful for drawing background grids or scatterplot smooths. hth, Matthias Jake Michaelson wrote: Hi all, I would like to be able to call a custom function automatically before plot.new() is called (more specifically, before a new plot is created on the current graphics device). Recently I've been poking around in the help files of some of the low(er) level plotting functions, and I seem to remember something saying that you could somehow add a setting to call a function immediately before plot.new() is called (something kind of in the spirit of .First). I've been looking everywhere, and I can't find the documentation I'm thinking of, which makes me wonder if I'm just making all this up. Can this be done? Thanks, Jake [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] Sum of random values
Hello Ricardo, no, the distribution needs not to be predefined. You can define your own distribution. For more details you may want to have a look at the vignette which is provided by package distrDoc. Do the random variables follow a discrete or an absolutely continuous distribution? To see how convolution is performed in package distr have a look at library(distr) getMethod(+, c(DiscreteDistribution, DiscreteDistribution)) ## respectively getMethod(+, c(AbscontDistribution, AbscontDistribution)) hth, Matthias Ricardo Bessa wrote: Thank you Gregory, I know the package distR but the the percentiles that I have are from an unkonwn distribuition, and I think that the package distR only works with pre-defined distribuions. Best Regards, Ricardo Bessa Subject: RE: [R] Sum of random values Date: Fri, 25 Apr 2008 09:16:10 -0600 From: [EMAIL PROTECTED] To: [EMAIL PROTECTED]; r-help@r-project.org You might want to look at the distr package (and its relatives). It provides methods for calculating the distribution function of combinations (the sum is one) of other distributions. I'm not sure how you would convert your percentiles to a distribution function, but there may be a way in the documentation for distr and friends. -- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare [EMAIL PROTECTED] (801) 408-8111 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Ricardo Bessa Sent: Thursday, April 24, 2008 5:12 PM To: r-help@r-project.org Subject: [R] Sum of random values Hello, I have two random variables with their percentiles which correspond to their ! pr! obability distribution function. My objective is to sum these two random variables. There exists any algorithm or procedure in R capable of converting the percentiles to a probability density function? is the fast Fourier transform function of R(fft) capable of doing the sum with a convolution?I'm just starting with this specific problem, so any help it will be very useful.Best regards, Ricardo Bessa _ [[elided Hotmail spam]][[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. _ [[elided Hotmail spam]] [[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. -- Dr. Matthias Kohl www.stamats.de __ 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] inverse F distribution in R?
just a comment to fstat... Using package distrEx one also has fstat functionality library(distrEx) E(Fd(df1 = 3, df2 = 3)) E(Fd(df1 = 4, df2 = 4)) E(Fd(df1 = 5, df2 = 5)) var(Fd(df1 = 5, df2 = 5)) for more functionals (skewness, kurtosis, IQR, ...) see for instance help(var, package = distrEx) Best Matthias Peter Dalgaard wrote: Jennifer Balch wrote: Dear all, I'm looking for a function that calls the inverse F-distribution. Something equivalent to FINV in matlab or excel. Does anyone know if such a function already exists for R? (I haven't been able to find one.) Thanks for any leads. I would guess that they are qf() or a variant thereof. I gather that _fpdf http://www.mathworks.com/access/helpdesk/help/toolbox/stats/fpdf.html__, fcdf http://www.mathworks.com/access/helpdesk/help/toolbox/stats/fcdf.html, finv http://www.mathworks.com/access/helpdesk/help/toolbox/stats/finv.html, frnd http://www.mathworks.com/access/helpdesk/help/toolbox/stats/frnd.html i_n matlab is df, pf, qf, and rf in R. We don't have an equivalent of fstat for mean and variance of the F distribution (this might actually be nice to have). Best, Jennifer __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] mean for each row
or rowMeans ... Matthias Hesen Peng wrote: How about rowSums(x)/ncol(x), where x is the matrix? On Mon, Mar 17, 2008 at 1:48 PM, Roslina Zakaria [EMAIL PROTECTED] wrote: Hi r-users, How do find the mean for each row? Thank you in advance for your help. 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 Day Totals 10.0 0.0 0.0 0.0 0.6 0.0 8.4 0.0 29.4 0.0 38.4 20.0 0.0 1.8 0.0 22.4 0.0 0.2 0.4 0.8 0.0 25.6 37.8 0.0 0.0 17.6 1.4 0.0 0.0 0.0 0.0 0.0 26.8 42.2 0.8 0.4 0.0 0.2 11.2 1.4 33.2 0.0 0.0 49.4 50.2 1.8 0.0 1.0 0.0 0.2 0.0 12.2 0.0 19.2 34.6 60.0 0.0 0.0 1.0 0.0 0.0 0.0 2.2 0.0 14.6 17.8 70.0 0.0 0.0 0.0 3.6 0.2 0.0 2.0 0.0 0.26.0 80.0 0.0 10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 10.0 90.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00.0 10 1.0 1.4 0.0 0.0 0.0 1.0 0.4 0.0 0.0 0.03.8 11 0.0 8.2 0.0 0.0 0.0 0.0 8.4 0.0 0.0 0.4 17.0 12 10.8 0.8 0.0 0.0 0.0 0.0 1.0 0.0 0.0 4.2 16.8 13 32.8 0.0 0.0 0.8 0.0 0.0 0.0 0.0 0.2 0.2 34.0 14 1.0 0.0 1.6 0.2 0.0 0.0 0.0 0.0 0.0 0.02.8 15 0.0 0.0 0.0 2.2 1.4 0.0 0.0 0.0 0.0 0.03.6 16 0.6 0.0 0.0 0.6 22.0 0.0 0.0 0.0 0.0 0.0 23.2 17 0.0 0.0 0.0 0.0 0.0 0.0 2.8 0.0 0.0 0.02.8 18 2.8 0.0 0.0 0.0 0.0 0.0 8.2 0.0 0.0 8.2 19.2 19 0.0 0.0 0.0 0.0 0.0 0.2 0.0 0.0 0.0 0.20.4 20 0.0 8.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.08.2 21 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00.0 22 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.40.4 23 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.2 0.0 0.01.2 24 0.0 0.0 0.0 0.0 1.0 0.0 0.0 2.0 8.2 0.0 11.2 25 3.2 0.0 0.0 0.6 0.6 0.0 0.0 0.0 11.8 0.0 16.2 26 0.0 0.0 26.2 0.0 12.6 0.0 0.0 2.2 0.0 0.0 41.0 27 0.2 0.0 10.6 0.0 1.2 0.0 0.0 1.8 0.0 0.0 13.8 28 0.0 4.0 0.0 5.8 0.0 0.0 0.0 0.0 0.0 0.09.8 29 0.2 12.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 12.6 30 0.0 2.6 0.0 0.0 2.2 0.0 0.0 0.0 0.0 0.04.8 31 0.0 0.0 0.6 0.0 0.8 0.0 0.0 0.0 4.8 0.06.2 Year Totals 62.8 40.2 51.2 29.8 70.0 12.8 30.8 57.2 55.2 47.6 457.6 Be a better friend, newshound, and __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] multiple boxplots using boxplot
Hi Marek, use a list instead of a data.frame; e.g. boxplot(list(x = rnorm(10), y = rnorm(20))) hth, Matthias Marek Bartkuhn wrote: Hi, I have data in 3 vectors a,b and c looking like eg a:1.2 3.4 1.4 .. I like to have the data from each vector as a boxplot, but the 3 boxplots within one graph. One example is the first graph at http://www.statmethods.net/graphs/boxplot.html Unfortunately I do not really understand the corresponding description. The vectors are of different length, therefore R complains when I try to build a data.frame of the 3 vectors: arguments imply differing number of rows: 6175, 20, 190 What do I have to do? Marek __ 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. -- Dr. Matthias Kohl www.stamats.de __ 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] likelihood from test result
Dear David, no, distrTEst won't help. It has a different intention. We are currently working on a new package distrMod (cf. https://r-forge.r-project.org/projects/distrmod/) which sometime might have such a functionality. Best, Matthias David Bickel wrote: Is there any automatic mechanism for extracting a likelihood or test statistic distribution (PDF or CDF) from an object of class htest or from another object of a general class encoding a hypothesis test result? I would like to have a function that takes x, an object of class htest, as its only argument and that returns the likelihood or test statistic distribution that was used to compute the p-value. It seems the only way to write such a function is to manually assign each test its statistic's distribution, e.g., like this: FUN - if(names(x$statistic) == t) dt else if(names(x$statistic) == X-squared) dchisq # etc. Is there a general S3 or S4 class other than htest that would better accommodate such extraction of distributions or likelihoods? I would also appreciate any suggestions for strategies or contributed packages that may facilitate automation. For example, would the distrTEst package help? David __ David R. Bickel Ottawa Institute of Systems Biology http://www.oisb.ca/members.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. __ 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] probability from different values
one addition ... Hello Ricardo, another solution could be using package distr: library(distr) A - c(18,18,18,19,20,21,22,23,24,25,26,27,28) DA - DiscreteDistribution(A) # maybe # support(DA) # plot (DA) B - c(82,83,84,85,85,86,87,88,89,90,91,91,92) DB - DiscreteDistribution(B) # support(DB) # plot(DB) DC - DB - DA# convolution with of DB with (-DA) # support(DC) # plot(DC) to get the probabilities for the C-values: p(DC)(support(DC)) hth, Matthias [EMAIL PROTECTED] wrote: i have two vectors (A and B) and i need create another vector (C) from the subtraction of A's values with the B's values. How can i estimate the probability of C's values if i have differents values combinations of A and B that can result in the same value? like 90 -20 = 70 and 91 - 21 = 70. example: B = {18 1818 19 20 21 22 23 24 25 26 27 28} A = {82 83 84 85 85 86 87 88 89 90 91 9192 } Try this: A = c(18,18,18,19,20,21,22,23,24,25,26,27,28) B = c(82,83,84,85,85,86,87,88,89,90,91,91,92) lenA = length(A) lenB = length(B) AA = rep(A, lenB) BB = rep(B, each=lenA) table(AA-BB)/(lenA*lenB) Regards, Richie. Mathematical Sciences Unit HSL ATTENTION: This message contains privileged and confidential inform...{{dropped:20}} __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] probability from different values
Hello Ricardo, another solution could be using package distr: library(distr) A - c(18,18,18,19,20,21,22,23,24,25,26,27,28) DA - DiscreteDistribution(A) # maybe # support(DA) # plot (DA) B - c(82,83,84,85,85,86,87,88,89,90,91,91,92) DB - DiscreteDistribution(B) # support(DB) # plot(DB) DC - DB - DA# convolution with of DB with (-DA) # support(DC) # plot(DC) hth, Matthias [EMAIL PROTECTED] wrote: i have two vectors (A and B) and i need create another vector (C) from the subtraction of A's values with the B's values. How can i estimate the probability of C's values if i have differents values combinations of A and B that can result in the same value? like 90 -20 = 70 and 91 - 21 = 70. example: B = {18 1818 19 20 21 22 23 24 25 26 27 28} A = {82 83 84 85 85 86 87 88 89 90 91 9192 } Try this: A = c(18,18,18,19,20,21,22,23,24,25,26,27,28) B = c(82,83,84,85,85,86,87,88,89,90,91,91,92) lenA = length(A) lenB = length(B) AA = rep(A, lenB) BB = rep(B, each=lenA) table(AA-BB)/(lenA*lenB) Regards, Richie. Mathematical Sciences Unit HSL ATTENTION: This message contains privileged and confidential inform...{{dropped:20}} __ 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] Simulate an AR(1) process via distributions? (without specifying a model specification)
Dear Pedro, you might be interested in the demo StationaryRegressorDistr of package distr. library(distr) demo(StationaryRegressorDistr) hth, Matthias [EMAIL PROTECTED] wrote: Thanks Prof. Ripley. My apologies for not including the code. Below I illustrate my point using the GLD package. Thank you very much for your time. Kind Regards, Pedro N. Rodriguez # Code begins # Simulate an ar(1) process # x = 0.05 + 0.64*x(t-1) + e # Create the vector x x - vector(length=1000) #simulate the own risk e - rnorm(1000) #Set the coefficient beta- 1.50 # set an initial value x[1]- 5 #Fill the vector x for(i in 2:length(x)) { x[i]- 0.05 + beta*x[i-1] + e[i] } #Check the AR(1) simulated_data_ar - arima(x,order=c(1,0,0)) simulated_data_ar #Using the G Lambda Distribution to fit the distribution. library(gld) resul1 - starship(x,optim.method=Nelder-Mead) lambdas1- resul1$lambda #Plot the Distribution plotgld(lambdas1[1],lambdas1[2],lambdas1[3],lambdas1[4]) #Random Deviates from GLD x_sim - rgl(1000,lambdas1[1],lambdas1[2],lambdas1[3],lambdas1[4]) #Fit an AR(1) gld_simulated - arima(x_sim,order=c(1,0,0)) gld_simulated #Code ends -Original Message- From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] Sent: Wednesday, November 28, 2007 11:37 AM To: Rodriguez, Pedro Cc: [EMAIL PROTECTED] Subject: Re: [R] Simulate an AR(1) process via distributions? (without specifying a model specification) On Wed, 28 Nov 2007, [EMAIL PROTECTED] wrote: Is it possible to simulate an AR(1) process via a distribution? Any distribution *of errors*, yes. Of the process values, not in general. I have simulated an AR(1) process the usual way (that is, using a model specification and using the random deviates in the error), and used the generated time series to estimate 3- and 4-parameter distributions (for instance, GLD). However, the random deviates generated from these distributions do not follow the specified AR process. How do you know that? Please give us the reproducible example we asked for (in the posting guide, at the bottom of every message), and we should be able to explain it to you. __ 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] skewness: generating distribution and estimating M3.
Dear Milton, package distrEx has methods to compute skewness. library(distrEx) ?skewness To generate new distributions you could use package distr. hth, Matthias Milton Cezar Ribeiro wrote: Dear all, How can I generate values which distrituions have left or rigth skewness. After that, which function compute the skewness of each distributions? Kind regards. miltinho para armazenamento! [[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] analytical solution to Sum of binominal distributed random numbers?
Hello Rainer, You could also get a very exact approximation using our package distr; e.g. library(distr) B1 - Binom(prob = 0.4, size = 20) B2 - Binom(prob = 0.2, size = 10) B3 - B1+B2 #convolution! plot(B3) hth, Matthias Rainer M. Krug schrieb: Frede Aakmann Tøgersen wrote: Perhaps http://stinet.dtic.mil/cgi-bin/GetTRDoc?AD=ADA266969Location=U2doc=GetTRDoc.pdf is something that you can use? Thanks a lot - that might help. Rainer Best regards Frede Aakmann Tøgersen Scientist UNIVERSITY OF AARHUS Faculty of Agricultural Sciences Dept. of Genetics and Biotechnology Blichers Allé 20, P.O. BOX 50 DK-8830 Tjele Phone: +45 8999 1900 Direct: +45 8999 1878 E-mail: [EMAIL PROTECTED] Web:http://www.agrsci.org This email may contain information that is confidential. Any use or publication of this email without written permission from Faculty of Agricultural Sciences is not allowed. If you are not the intended recipient, please notify Faculty of Agricultural Sciences immediately and delete this email. -Oprindelig meddelelse- Fra: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] På vegne af Rainer M Krug Sendt: 24. oktober 2007 09:11 Til: Charles C. Berry Cc: r-help Emne: Re: [R] analytical solution to Sum of binominal distributed random numbers? Hi Charles thanks for the pointing out that size and prob can be vectors as well - I tried it out but used 1 as the number of observations, assuming that and it only gave me one randon mumbewr (as it should be but not expected). But I was more looking at a analytical solution, as I have to sum up a huge number of random numbers. But I am going to try your solution as it should be much faster already. Thanks Rainer Charles C. Berry wrote: ?rbinom only says: size: number of trials (zero or more). prob: probability of success on each trial. But they can be vectors. BTW, you were aked to PLEASE ... provide minimal, self-contained, reproducible code. What you show cannot run without correction. Most likely, you intended size(n) to be the n-th element of the vector 'size', which in R is written 'size[ n ]' . In which case sum (rbinom( length(prob) , size, prob ) ) works. Chuck On Tue, 23 Oct 2007, Rainer M Krug wrote: Hi I have two vectors, prob and size, and I want to add the random deviates of these two, i.e. sum( sapply( 1:length(prob), function(n){ rbinom(1, size(n), prob(n) } ) ) My problem is that I have to do this for a large number of value combinations. Is there a faster way of doing this? Rainer __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Charles C. Berry(858) 534-2098 Dept of Family/Preventive Medicine E mailto:[EMAIL PROTECTED] UC San Diego http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ 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.