[R] Updating a worksheet in Excel file using RODBC
Hello! I have no problem reading Excel files (each worksheet in the file is a table which can be read - at least in my case). What I would like to do is to read such a table, change it (just the contents, not the format) and write it back, and this I can not do. I am getting the following error messages (3 slightly different attempts): sqlSave(con, x, tablename = Chimaera20_3years$, append = FALSE, + rownames = FALSE, colnames = TRUE, + verbose = TRUE, oldstyle = FALSE,safer=FALSE) Query: CREATE TABLE Chimaera20_3years$ (Date varchar(255), 000Tax varchar(255), 1500Tax varchar(255), 3000Tax varchar(255), 4650Tax varchar(255)) Error in sqlSave(con, x, tablename = Chimaera20_3years$, append = FALSE, : [RODBC] ERROR: Could not SQLExecDirect 37000 -3551 [Microsoft][ODBC Excel Driver] Syntax error in CREATE TABLE statement. sqlSave(con, x, tablename = [Chimaera20_3years$], append = FALSE, + rownames = FALSE, colnames = TRUE, + verbose = TRUE, oldstyle = FALSE,safer=FALSE) Query: CREATE TABLE [Chimaera20_3years$] (Date varchar(255), 000Tax varchar(255), 1500Tax varchar(255), 3000Tax varchar(255), 4650Tax varchar(255)) Error in sqlSave(con, x, tablename = [Chimaera20_3years$], append = FALSE, : [RODBC] ERROR: Could not SQLExecDirect 37000 -3553 [Microsoft][ODBC Excel Driver] Syntax error in field definition. sqlSave(con, x, tablename = [Chimaera20_3years], append = FALSE, + rownames = FALSE, colnames = TRUE, + verbose = TRUE, oldstyle = FALSE,safer=FALSE) Query: CREATE TABLE [Chimaera20_3years] (Date varchar(255), 000Tax varchar(255), 1500Tax varchar(255), 3000Tax varchar(255), 4650Tax varchar(255)) Error in sqlSave(con, x, tablename = [Chimaera20_3years], append = FALSE, : [RODBC] ERROR: Could not SQLExecDirect 37000 -3553 [Microsoft][ODBC Excel Driver] Syntax error in field definition. Am I doing it wrong way or is there a problem with the Excel driver? Thank you in advance, Moshe Olshansky Chimaera Capital Group Moshe Olshansky Chimaera Capital Limited Level 4 / 349 Collins Street Melbourne, Victoria 3000 Phone: +613 8614 8400 Fax: +613 8614 8410 Email: [EMAIL PROTECTED] Disclaimer: This message is intended only for the personal and confidential use of the designated recipient(s) named above. If you are not the intended recipient of this message you are hereby notified that any review, dissemination, distribution or copying of this message is strictly prohibited. This communication is for information purposes only and should not be regarded as an offer to sell or as a solicitation of an offer to buy any financial product, an official confirmation of any transaction, or as an official statement of Chimaera Capital Limited. E-mail transmissions cannot be guaranteed to be secure or error-free. Therefore, we do not represent that this information is complete or accurate and it should not be relied upon as such. All information is subject to change without notice. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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] Updating a worksheet in Excel file using RODBC
Dear Hans-Peter, Thank you very much for your note! I tried your package and it works all right (i.e. it indeed writes data to Excel files), however it creates a new Excel file and this is not what I really need. I need to update/create one sheet in the existing file. I am using R do compute some data but then it must be put into an Excel file and an Excel chart must be created. So I intended to (manually) create one excel file and to write a VBA macro which makes a desirable chart of a certain sheet (let's say Sheet1). Then I intended to autoomatically make many copies of this file with appropriate file names, write the right data to Sheet1 of each such file so that when it is opened a desired chart is automatically created. So creting a totally new file does not help me. Best regards, Moshe. -Original Message- From: Hans-Peter [mailto:[EMAIL PROTECTED] Sent: Saturday, 24 March 2007 3:44 AM To: Moshe Olshansky Cc: R Help Subject: Re: [R] Updating a worksheet in Excel file using RODBC Hi, 2007/3/23, Moshe Olshansky [EMAIL PROTECTED]: Hello! I have no problem reading Excel files (each worksheet in the file is a table which can be read - at least in my case). What I would like to do is to read such a table, change it (just the contents, not the format) and write it back, and this I can not do. I am getting the following error messages (3 slightly different attempts): [snip] As another option (if you work with Windows) you can check my xlsReadWrite package (- CRAN). It should work very well in your case (it's not suited if you want to use SQL (join) statements, but for plain data reading/writing it is nice). For both versions (free/pro) updates are pending. They should be released by end of next week (but no guarantees). -- Regards, Hans-Peter Moshe Olshansky Chimaera Capital Limited Level 4 / 349 Collins Street Melbourne, Victoria 3000 Phone: +613 8614 8400 Fax: +613 8614 8410 Email: [EMAIL PROTECTED] Disclaimer: This message is intended only for the personal and confidential use of the designated recipient(s) named above. If you are not the intended recipient of this message you are hereby notified that any review, dissemination, distribution or copying of this message is strictly prohibited. This communication is for information purposes only and should not be regarded as an offer to sell or as a solicitation of an offer to buy any financial product, an official confirmation of any transaction, or as an official statement of Chimaera Capital Limited. E-mail transmissions cannot be guaranteed to be secure or error-free. Therefore, we do not represent that this information is complete or accurate and it should not be relied upon as such. All information is subject to change without notice. __ R-help@stat.math.ethz.ch 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] Updating a worksheet in Excel file using RODBC
Dear Prof. Ripley, You seem not to have tried the simplest possible option. The following works for me (beware of wrapped lines from mailers) chan - odbcDriverConnect(DRIVER=Microsoft Excel Driver (*.xls);DBQ=C:\\bdr\\hills.xls; ReadOnly=False) sqlSave(chan, USArrests, tests, fast=TRUE) # or FALSE You are right - I have not. It does not work exactly as it should have but this solves my problem. I created a very small Excel file odbc1.xls containing 3 sheets (test, Sheet2 and Sheet3). Below is a short R session: chan - odbcDriverConnect(DRIVER=Microsoft Excel Driver (*.xls);DBQ=C:\\EFGraphs\\odbc1.xls; ReadOnly=False) x- c(1:6) x - matrix(x,nrow=3,ncol=2) x - data.frame(x) x X1 X2 1 1 4 2 2 5 3 3 6 sqlSave(chan, x, test, fast=FALSE) Error in sqlSave(chan, x, test, fast = FALSE) : table 'test' already exists sqlSave(chan, x, tests, fast=FALSE) As you see I was unable to overwrite an existing sheet (an attempt to drop this table also fails), but I was able to add a new sheet to an existing Excel file (after this action the file contains 4 sheets - the 3 it contained and the last sheet named tests). This allows me to do what I wanted, i.e. manually create an Excel file with a small VBA macro, make many copies of this file (under appropriate names), write an appropriate data to each file and then the macro will work on the right data (different for each file). Thanky you! Moshe. Moshe Olshansky Chimaera Capital Limited Level 4 / 349 Collins Street Melbourne, Victoria 3000 Phone: +613 8614 8400 Fax: +613 8614 8410 Email: [EMAIL PROTECTED] Disclaimer: This message is intended only for the personal and confidential use of the designated recipient(s) named above. If you are not the intended recipient of this message you are hereby notified that any review, dissemination, distribution or copying of this message is strictly prohibited. This communication is for information purposes only and should not be regarded as an offer to sell or as a solicitation of an offer to buy any financial product, an official confirmation of any transaction, or as an official statement of Chimaera Capital Limited. E-mail transmissions cannot be guaranteed to be secure or error-free. Therefore, we do not represent that this information is complete or accurate and it should not be relied upon as such. All information is subject to change without notice. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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] Updating a worksheet in Excel file using RODBC
OK. By the way, I only thought that I could do what I wanted! It worked once but then it failed. When I was trying to update an existing sheet I got an error message saying that it existed and when I was trying to make a new sheet (something that worked once) I got a message saying that there was no such table! -Original Message- From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] Sent: Tuesday, 27 March 2007 4:44 PM To: Moshe Olshansky Cc: [EMAIL PROTECTED] Subject: Re: Updating a worksheet in Excel file using RODBC Yes, sqlDrop does not work correctly for Excel worksheet names (and there are other quirks). As I said in another message, it is on my TODO list to make this work better, but in the absence of good documentation of what the Excel ODBC driver should do and several with known bugs it is largely a matter of trial-and-error. On Tue, 27 Mar 2007, Moshe Olshansky wrote: Dear Prof. Ripley, You seem not to have tried the simplest possible option. The following works for me (beware of wrapped lines from mailers) chan - odbcDriverConnect(DRIVER=Microsoft Excel Driver (*.xls);DBQ=C:\\bdr\\hills.xls; ReadOnly=False) sqlSave(chan, USArrests, tests, fast=TRUE) # or FALSE You are right - I have not. It does not work exactly as it should have but this solves my problem. I created a very small Excel file odbc1.xls containing 3 sheets (test, Sheet2 and Sheet3). Below is a short R session: chan - odbcDriverConnect(DRIVER=Microsoft Excel Driver (*.xls);DBQ=C:\\EFGraphs\\odbc1.xls; ReadOnly=False) x- c(1:6) x - matrix(x,nrow=3,ncol=2) x - data.frame(x) x X1 X2 1 1 4 2 2 5 3 3 6 sqlSave(chan, x, test, fast=FALSE) Error in sqlSave(chan, x, test, fast = FALSE) : table 'test' already exists sqlSave(chan, x, tests, fast=FALSE) As you see I was unable to overwrite an existing sheet (an attempt to drop this table also fails), but I was able to add a new sheet to an existing Excel file (after this action the file contains 4 sheets - the 3 it contained and the last sheet named tests). This allows me to do what I wanted, i.e. manually create an Excel file with a small VBA macro, make many copies of this file (under appropriate names), write an appropriate data to each file and then the macro will work on the right data (different for each file). Thanky you! Moshe. Moshe Olshansky Chimaera Capital Limited Level 4 / 349 Collins Street Melbourne, Victoria 3000 Phone: +613 8614 8400 Fax: +613 8614 8410 Email: [EMAIL PROTECTED] Disclaimer: This message is intended only for the personal and confidential use of the designated recipient(s) named above. If you are not the intended recipient of this message you are hereby notified that any review, dissemination, distribution or copying of this message is strictly prohibited. This communication is for information purposes only and should not be regarded as an offer to sell or as a solicitation of an offer to buy any financial product, an official confirmation of any transaction, or as an official statement of Chimaera Capital Limited. E-mail transmissions cannot be guaranteed to be secure or error-free. Therefore, we do not represent that this information is complete or accurate and it should not be relied upon as such. All information is subject to change without notice. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch 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] Creating an Access (.mdb) database using R
Hello! I have a short question: Is it possible to create a (non-existing) Access database using R (and if yes, how)? I need to create a new database and then insert a few tables into it. Thank you in advance, Moshe Olshansky [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch 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] Stock Price Correlation to Index Price Levels
Josh Kalish wrote: Thanks to all of the people who responded. What I was trying to do is to turn my matrix or frame containing index level returns and stock returns into a matrix of betas. I don't really need to worry about risk-free interest rates. I just need to be able to come up with a number that shows the expected index correlation. I was able finally to figure out how to use cor() to get what I think is an R^2 value. But, I'm trying to also figure out the ratio of correlation. For example, some stocks correlate very well and cor() returns a value of .92. But, how do you then figure out if the stock should have a 1.5:1 correlation? The way I would do it by hand is to turn the closes into daily returns and then get the mean() return for each stock against the index by day. I can't find an example as hard as I look, but this must be very common. If X is a vector of daily returns (today's close - yesterday's close) for your index (in dollars) and Y is the vector of daily returns for your stock (in dollars), then to hedge 1 share of your stock you need to hold u = -(1/r)*(Sy/Sx) units of the index, where r is the correlation coefficient, Sx is the standard deviation of daily returns of the index and Sy is the standard deviation of daily returns of Y. In R, u - -(1/cor(X,Y))*(sd(Y)/sd(X)) In this case your stock is fully hedged by the index (but the index of course is not fully hedged by the stock, unless the correlation coefficient is +/- 1). Hope this helps, Moshe Olshansky. __ R-help@stat.math.ethz.ch 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] Appropriate regression model for categorical variables
Tirtha wrote: Dear users, In my psychometric test i have applied logistic regression on my data. My data consists of 50 predictors (22 continuous and 28 categorical) plus a binary response. Using glm(), stepAIC() i didn't get satisfactory result as misclassification rate is too high. I think categorical variables are responsible for this debacle. Some of them have more than 6 level (one has 10 level). Please suggest some better regression model for this situation. If possible you can suggest some article. thanking you. Tirtha Hi Tirtha, Are your categorical variables really categorical? What I mean is if you variable is user's satisfaction level (0 for very unsatisfied, 1 for moderately unsatisfied, 2 for slightly unsatisfied, 4 for neutral, etc., finally 7 for very satisfied) then your variable is not really categorical (since 1 is closer to 3 than to 6) and then try what other people suggest. However, if your variable is, say, the 50-th amino acid in a certain gene (with values of 1 for the first amino acid, 2 for the second one,...,20 for the 20-th one) then your variable is really categorical (you generally can not say that amino acid 2 is much closer to amino acid 3 than to amino acid 17). In such a case I would have tried classification method which can treat categorical variables or, alternatively, may be regression trees (i.e. split on the values of categorical variables and at each node find regression coefficients of the continuous variables). Regards, Moshe Olshansky [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch 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] Responding to a posting in the digest
Is there a convenient way to respond to a particular posting which is a part of the digest? I mean something that will automatically quote the original message, subject, etc. Thank you! Moshe Olshansky [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch 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] if within a function
Hi Yuchen, In R, if you do not put an explicit return statement in the function, the value the function returns is the value of the last statement in the function. Unlike VB, it does not matter whether you assign this value to aaa (which is identical to your function name) or b or c or x etc. So either use an explicit return statement or make sure that the last statement in the function produces the right result. --- Yuchen Luo [EMAIL PROTECTED] wrote: Dear Friends. I found a puzzling phenomenon in R when you use 'if' within a function: # defining a function aaa aaa=function(a) {if (a==1) {aaa=1}; if (a!=1) {aaa=2} } # using the function: b=20 bbb=aaa(b) bbb [1] 2 typeof(bbb) [1] double c=1 ccc=aaa(c) ccc NULL typeof(ccc) [1] NULL It seems that only the last 'if' phrase works. Is it an instrinsic weakness of R? Is there a way to get around it? ( I use 'elseif' to get around this when there are only two cases to choose from, but what if there are more than two cases to choose from?) Best Yuchen [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Speed up R
Robert McFadden writes: Dear R Users, I hope that there is someone who has an experience with a problem that I describe below and will help me. I must buy new desktop computer and I'm wondering which processor to choose if my only aim is to speed up R. I would like to reduce a simulation time - sometimes it takes days. I consider buying one of them (I'm working under Win XP 32 bit): 1. Intel Core2 Duo E6700 2.67 GHz 2. Dual-Core Intel Xeon processor 3070 - 2,66 GHz 3. AMD Athlon 64 X2 6000+ Or simple Pentium 4? I'm very confused because I'm not sure whether R takes advantage dual-core or not. If not, probably Athlon would be better, wouldn't be? I would appreciate any help. Rob Hi Robert, Let me suggest you a dirty solution - if simulations take days and you must run them many times I would have rewriten them, let say, in C. I had a program in Matlab which took more than an hour to run and I had to run it many times, so I usually prepared a few runs, started them in the evening before leaving the office and got the results next morning. After a while I have re-written it in C (this took me a few days) and got a spead-up factor of about 100, so that now the run took just a few minutes! Languages like R and Matlab are extreemely convenient but if performance is a very important issue you shoul use C, Fortran, C++, etc. Regards, Moshe Olshansky. __ R-help@stat.math.ethz.ch 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] Binary decision problem
Hi Tolga, I do not see any problem with: max {a1*x1 + a2*x2 + ... + a800*x800} subject to: x1+x2+ ... + x800 = 100 b1*x1+b2*x2+ ... +b800*x800 = B c1*x1+c2*x2+ ... +c800*x800 = C and an additional condition that x1,x2,...,x800 are binary 0-1. Regards, Moshe Olshansky --- Tolga Uzuner [EMAIL PROTECTED] wrote: Dear R Users, I am trying to use LP_SOLVE and would appreciate any assistance with the following problem: - I am trying to choose a fixed number of items out of a batch of items: say 100 out of 800 - items have certain characteristics, say characteric A, B and C - I want to maximise the sum of A across all 100 items I choose while ensuring that the sum of B and C across the items do not exceed certain constraints How exactly do I set this up in lp_solve ? If I associate a boolean, 0 or 1, with each item, I can constrain the sum of the boolean to be equal to 100. But how do I then further specify the other constraints (on the sum of B and the sum of C) and the objective function (to maximise the sum of A) ? Thanks, Tolga __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] anova on data means
An ugly solution (but which works) would be something like this (assuming that d is a data.frame with 4 columns named Tank, Trt, Fish, Size): x -1:12 for (i in 1:12) x[i] - mean(d$Size[d$Tank == i]) Regards, Moshe Olshansky --- Ronaldo Reis Junior [EMAIL PROTECTED] wrote: Em Quinta 21 Junho 2007 16:56, Thomas Miller escreveu: I am transitioning from SAS to R and am struggling with a relatively simple analysis. Have tried Venables and Ripley and other guides but can't find a solution. I have an experiment with 12 tanks. Each tank holds 10 fish. The 12 tanks have randomly assigned one of 4 food treatments - S(tarve), L(ow), M(edium) and H(igh). There are 3 reps of each treatment. I collect data on size of each fish at the end of the experiment. So my data looks like Tank Trt Fish Size 1 S 1 3.4 1 S 2 3.6 1 S10 3.5 2 L 1 3.4 12M 10 2.1 To do the correct test of hypothesis using anova, I need to calculate the tank means and use those in the anova. I have tried using tapply() and by() functions, but when I do so I loose the treatment level because it is categorical. I have used Meandattapply(Size,list(Tank, Trt), mean) But that doesn't give me a dataframe that I can then use to do the actual aov analysis. So what is the most efficient way to accomplish the analysis Thanks Tom Miller Tom, try the aggregate funtion. Somethink like this meandat - aggregate(Size,list(Tank,Trt),mean) Inte Ronaldo -- Prof. Ronaldo Reis Júnior | .''`. UNIMONTES/Depto. Biologia Geral/Lab. de Ecologia | : :' : Campus Universitário Prof. Darcy Ribeiro, Vila Mauricéia | `. `'` CP: 126, CEP: 39401-089, Montes Claros - MG - Brasil | `- Fone: (38) 3229-8187 | [EMAIL PROTECTED] | [EMAIL PROTECTED] | http://www.ppgcb.unimontes.br/ | ICQ#: 5692561 | LinuxUser#: 205366 __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] sorting data
Try ReadWriteXls package from CRAN. --- dala [EMAIL PROTECTED] wrote: I have a 2 columns, Date and Number, in Excel. I copy and paste them into Notepad. I can use scan() to import the file but how do I plot this data with Date as the x-axis? -- View this message in context: http://www.nabble.com/sorting-data-tf3958889.html#a11233613 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R-excel
Also try xlsReadWrite package on CRAN. --- Erika Frigo [EMAIL PROTECTED] wrote: Good morning to everybody, I have a problem : how can I import excel files in R??? thank you very much Dr.sa. Erika Frigo Università degli Studi di Milano Facoltà di Medicina Veterinaria Dipartimento di Scienze e Tecnologie Veterinarie per la Sicurezza Alimentare (VSA) Via Grasselli, 7 20137 Milano Tel. 02/50318515 Fax 02/50318501 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] surprising difference in log()
This is not surprizing at all! The exact result is log(8,2) = 3, but the numerical procedure which calculates the logarithm may produce a result which is a few ULPs different from the exact one, i.e. you can get that log(8,2) = 2.99 and then floor(2.99) = 2. --- Fausto Galli [EMAIL PROTECTED] wrote: Hello everybody My collegue and I noticed a strange behaviour of R on different platforms. It's a simple computation, but results are rather different. On Windows XP: floor(log(8,2)) [1] 3 which is what one should expect. Here's instead the result with Mac OS X (same version, 2.5.0 (2007-04-23)) floor(log(8,2)) [1] 2 Is it a bug in R or in the operating system? Anyway, it's quite a surprising one. _ Fausto Galli Institute of Finance University of Lugano Via G. Buffi 13 CH-6904 Lugano, Switzerland. +41 (0)58 666 4497 http://www.people.lu.unisi.ch/gallif __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Gaussian elimination - singular matrix
All the nontrivial solutions to AX = 0 are the eigenvectors of A corresponding to eigenvalue 0 (try eigen function). The non-negative solution may or may not exist. For example, if A is a 2x2 matrix Aij = 1 for 1 =i,j =2 then the only non-trivial solution to AX = 0 is X = a*(1,-1), where a is any nonzero scalar. So in this case there is no non-negative solution. Let X1, X2,...,Xk be all the k independent eigenvectors corresponding to eigenvalue 0, i.e. AXi = 0 for i = 1,2,...,k. Any linear combination of them, X = X1,...,Xk, i.e. a1*X1 + ... + ak*Xk is also a solution of AX = 0. Let B be a matrix whose COLUMNS are the vectors X1,...,Xk (B = (X1,...,Xk). Then finding a1,...,ak for which all elements of X are non-negative is equivalent to finding a vector a = (a1,...,ak) such that B*a = 0. Of course a = (0,...,0) is a solution. The question whether there exists another one. Try to solve the following Linear Programming problem: max a1 subject to B*a = 0 (you can start with a = (0,...,0) which is a feasible solution). If you get a non-zero solution fine. If not try to solve min a1 subject to B*a = 0 if you still get 0 try this with max a2, then min a2, max a3, min a3, etc. If all the 2k problems have only 0 solution then there are no nontrivial nonnegative solutions. Otherwise you will find it. Instead of 2k LP (Linear Programming) problems you can look at one QP (Quadratic Programming) problem: max a1^2 + a2^2 + ... + ak^2 subject to B*a = 0 If there is a nonzero solution a = (a1,...,ak) then X = a1X1 +...+ak*Xk is what you are looking for. Otherwise there is no nontrivial nonnegative solution. --- Bruce Willy [EMAIL PROTECTED] wrote: I am sorry, there is just a mistake : the solution cannot be unique (because it is a vectorial space) (but then I might normalize it) can R find one anyway ? This is equivalent to finding an eigenvector in fact From: [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Date: Wed, 27 Jun 2007 22:51:41 + Subject: [R] Gaussian elimination - singular matrix Hello, I hope it is not a too stupid question. I have a singular matrix A (so not invertible). I want to find a nontrivial nonnegative solution to AX=0 (kernel of A) It is a special matrix A (so maybe this nonnegative solution is unique) The authors of the article suggest a Gaussian elimination method Do you know if R can do that automatically ? I have seen that solve has an option LAPACK but it does not seem to work with me :( Thank you very much _ Le blog Messenger de Michel, candidat de la Nouvelle Star : analyse, news, coulisses A découvrir ! [[alternative HTML version deleted]] _ Le blog Messenger de Michel, candidat de la Nouvelle Star : analyse, news, coulisses A découvrir ! [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] restructuring matrix
If your original matrix is A then unique(A$People) and unique(A$Desc) will produce a vector of different people and a vector of different descriptions. --- yoo [EMAIL PROTECTED] wrote: Hi all, let's say I have matrix PeopleDescValue Mary Height50 Mary Weight 100 FannyHeight 60 Fanny Height200 Is there a quick way to form the following matrix? People HeightWeight Mary 50 100 Fanny 60200 (Assuming I don't know the length of people/desc and let's say these are characters matrix.. I tried play with row(), col(), etc.. but I don't seem to find like a duplicate match function... I'm trying to write some one/two liner that convert my resulting matrix to vector and pick the appropriate fields.. etc ) Thanks! -- View this message in context: http://www.nabble.com/restructuring-matrix-tf3991741.html#a11334950 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Optimization
This is partially true since both the function to be maximized and the constraint are non-linear. One may substitute 1-x1-x2 for x3 and use (let say) Lagrange multipliers to get two non-linear equations with 2 unknowns for which there should be a function solving them. Then you must find the points where the constraint function intersects with the triangle {x1=0,x2=0,x1+x2=1}, which is easier (for each of the 3 edges you get a non-linear equation in one variable). So even though an (almost) analytical solution can be found it would be much more convenient to use an optimization function which (hopefully) does all this for you. Moshe. --- Ravi Varadhan [EMAIL PROTECTED] wrote: Hi, Your problem can be solved analytically. Eliminate one of the variables, say x3, from the problem by using the equality x1 + x2 + x3 = 1. Then solve for the intersection of the circle (in x1 and x2) defined by the radical constraint, with the straight line defined by the objective function. There will be, at most, two intersection points. The extremum has to be one of these two points, provided they also satisfy the other inequalities (To me, this sounds an awful lot like a homework problem). Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of massimiliano.talarico Sent: Monday, July 16, 2007 4:50 PM To: r-help Subject: [R] Optimization Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Date problem
Try the xlsReadWrite package from CRAN. --- Mr Natural [EMAIL PROTECTED] wrote: Dear Natalia: I have had the same problem. Solved it by keeping my dates in Excel files with the rest of my data. When I want to do something in R, highlight the date column and format it as general You will get a 5 digit number that represents the Julian date (more digits if you have minutes and seconds). Save the Excel file as a csv file. Then read it into R making the variables all numeric, with use of columnClasses, eg. from a big file that Im working on now, ss-read.csv(six.sites.epn.hist.csv,strip.white=T,na.strings=., colClasses=c(numeric,numeric,numeric, numeric,numeric,numeric, numeric,numeric,numeric, numeric,numeric,numeric, numeric,numeric,numeric, numeric,numeric,numeric,na.rm=T )) Then you make the Julian dates back into calendar dates that plot so nicely in R. In my file, I have one date column for every y variable because each has a different sequence of dates. Note that Excel and R have a different Julian day 1. Excel is Jan 1 1900, and R is Jan 1 1970. You set the day one Julian date in R with the origin feature in the dates(). The six interesting y variables (MP,C,D,UD,LD,and BS) each have two date columns, the most interesting of which is the calendar date, which are made from the Julian dates that I read in. dmp-dates(ss$MPdate,origin=c(month = 1, day = 1, year = 1900)) dc-dates(ss$Cdate,origin=c(month = 1, day = 1, year = 1900)) dd-dates(ss$Ddate,origin=c(month = 1, day = 1, year = 1900)) dld-dates(ss$LDdate,origin=c(month = 1, day = 1, year = 1900)) dud-dates(ss$UDdate,origin=c(month = 1, day = 1, year = 1900)) dbs-dates(ss$BSdate,origin=c(month = 1, day = 1, year = 1900)) Then, I plot a stack of 6 time series, one on top of each other. Note that I used arrow for error bars. #stack of plots, one for each site. attach(ss) par(mar=c(2,2,1,1), mfrow = c(6,1)) plot(MP~dmp, type=b,xlab=,ylim=c(0,1),xaxt = n) arrows(dmp,MP+semp,dmp,MP-semp,length=.02,angle=90,code=3) plot(C~dc, type=b,xlab=,ylim=c(0,1),xaxt = n) arrows(dc,C+sec,dc,C-sec,length=.02,angle=90,code=3) plot(D~dd, type=b,xlab=,ylim=c(0,1),xaxt = n) arrows(dd,D+sed,dd,D-sed,length=.02,angle=90,code=3) plot(LD~dld, type=b,xlab=,ylim=c(0,1),xaxt = n) arrows(dld,LD+seld,dld,LD-seld,length=.02,angle=90,code=3) plot(UD~dud, type=b,xlab=,ylim=c(0,1),xaxt = n) arrows(dud,UD+seud,dud,UD-seud,length=.02,angle=90,code=3) plot(BS~dbs, type=b,xlab=,ylim=c(0,1)) arrows(dbs,BS+sebs,dbs,BS-sebs,length=.02,angle=90,code=3) I hope this helps, Don natalia norden wrote: Hello, I have some stupid problems managing date data. I have a colomn date, which I converted from a character representation: for example: a=26/02/06 date=strptime(a,format=%d/%m/%y) For one part of the analysis, I'm interested only in the month and the year, so I did: m.y=strftime(date,format=%m/%y) This returns me 02/06, but this is an object of class character and I can't convert it into an object of class Date. Doing strptime(m.y,format=%m/%y) or as.Date(m.y,format=%m/%y) returns me NA How can a convert this colomn m.y in a Date class? Actually, I need it to plot fruiting data against time (month and year). Because I have many values of seeds in a month, I used the function tapply: seeds=tapply(data$seed,data$m.y,sum) But plot(x=names(seeds),y=seeds) doesn't work. Does anyone know an easier way? Thank you for your time, natalia __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- View this message in context: http://www.nabble.com/Date-problem-tf1280528.html#a11660658 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Optimization
You are right!!! For some strange reason I substituted ^ (exponentiation) for *, so the problem became Max x1^0.021986+x2^0.000964+x3^0.02913 with these conditions: x1+x2+x3=1; sqrt((x1^0.114434)^2+(x2^0.043966)^2+(x3^0.100031)^2)=0.04; which is clearly non-linear. --- Berwin A Turlach [EMAIL PROTECTED] wrote: G'day Moshe, On Tue, 17 Jul 2007 17:32:52 -0700 (PDT) Moshe Olshansky [EMAIL PROTECTED] wrote: This is partially true since both the function to be maximized and the constraint are non-linear. I am not sure what your definition of non-linear is, but in my book, and I believe by most mathematical/statistical definitions, the objective function is linear. The only non-linearity comes in through the second constraint. One may substitute 1-x1-x2 for x3 and use (let say) Lagrange multipliers to get two non-linear equations with 2 unknowns for which there should be a function solving them. Why would you want to use Lagrange multipliers? Isn't that a bit of an overkill? Once you substitute 1-x1-x2 for x3 in the second constraint, you have a quadratic equations in x1 and x2. So for any given value of x1 you can solve for x2 (or for any given value of x2 you can solve for x1). They still teach how to solve quadratic equations at school, don't they? ;-) Then you must find the points where the constraint function intersects with the triangle {x1=0,x2=0,x1+x2=1}, which is easier (for each of the 3 edges you get a non-linear equation in one variable). Even easier. Take an x1 between 0 and 1. If for that x1 the quadratic equation in x2 has no real solution, then x1 is not feasible. Otherwise find the values of x2 that solve the equation. Use each of these values together with x1 to calculate corresponding values of x3. Then check these tuples for feasibility. If they are feasible, evaluate the objective function and return the tuple with the larger function value. All the calculations outlined in the paragraph above are easily implemented in R, e.g. the function polyroot() returns the roots of a polynomial. Cheers, Berwin === Full address = Berwin A TurlachTel.: +65 6515 4416 (secr) Dept of Statistics and Applied Probability +65 6515 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED] Singapore 117546 http://www.stat.nus.edu.sg/~statba __ R-help@stat.math.ethz.ch 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] Linear programming question
Hi Tobias, Could you please explain the role of x and y - are they somehow related to S1,S2,S3,S4? Are they constant? Are they additional variable? What was your original problem (without the slack variables)? Regards, Moshe. --- Tobias Schlottmann [EMAIL PROTECTED] wrote: Hi everybody, consider please an optimization problem: minimize sum S1+S2 Subject to : y - x = A + S1 x - y = A + S2 and we want to add two more constraints: y - x = B - S3 x - y = B - S4 where A is a small constant value and B is a large constant value, S1 and S2 are surplus and S3 and S4 are slack variables. S3 and S4 have to be maximized in objective function. As objective function, is this correct? : minimize sum S1+ S2 - S3 -S4 where actually we want to minimize S1 and S2; and maximize S3 and S4. If it is not correct, what to do ? Thank you for any guide. Tobias - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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 heatmap - how to remove annoying X before numeric values?
Hi Suzanne, My solution (which I am sure is not the best) would be: heat - read.table('temp.txt') heat X1905 X1910 X1950 X1992 X2011 X2020 Gnat 0.08 0.29 0.29 0.37 0.39 0.43 Snake 0.16 0.34 0.32 0.40 0.41 0.53 Bat0.40 0.54 0.52 0.60 0.60 0.63 Cat0.16 0.27 0.29 0.39 0.37 0.41 Dog0.43 0.54 0.52 0.61 0.60 0.62 Lynx 0.50 0.57 0.54 0.59 0.50 0.59 a-names(heat) b-strsplit(a,split=X) w-unlist(b) w [1] 1905 1910 1950 1992 2011 2020 z - w[seq(2,length(w),by=2)] z [1] 1905 1910 1950 1992 2011 2020 names(heat) - z heat 1905 1910 1950 1992 2011 2020 Gnat 0.08 0.29 0.29 0.37 0.39 0.43 Snake 0.16 0.34 0.32 0.40 0.41 0.53 Bat 0.40 0.54 0.52 0.60 0.60 0.63 Cat 0.16 0.27 0.29 0.39 0.37 0.41 Dog 0.43 0.54 0.52 0.61 0.60 0.62 Lynx 0.50 0.57 0.54 0.59 0.50 0.59 Regards, Moshe. --- Suzanne Matthews [EMAIL PROTECTED] wrote: Hello All, I have a simple question based on how things are labeled on my heat map; particularly, there is this annoying X that appears before the numeric value of all the labels of my columns. Let's say I have the following silly data, stored in temp.txt 190519101950199220112020 Gnat0.080.290.290.370.390.43 Snake 0.160.340.320.400.410.53 Bat 0.400.540.520.600.600.63 Cat 0.160.270.290.390.370.41 Dog 0.430.540.520.610.600.62 Lynx0.500.570.540.590.5 0.59 I use the following commands to generate my heatmap: heat - read.table('temp.txt') x - as.matrix(heat) heatmap.2(x, keysize=1.2, dendrogram=none, trace=none, Colv = FALSE, main = Silly Data, labCol= NULL, margin=c(7,8)) This generates a very nice heatmap, but there is one thing I have an issue with: How do I get rid of the 'X' that seems to come automatically before my numeric column values? I just want those columns to be labeled 1905, 1910, 1950, and so on. I cannot find anything in the heatmap.2 documentation that suggests how I should do this. Thank you very much for your time, and patience in reading this! Sincerely, Suzanne [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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 heatmap - how to remove annoying X before numeric values?
I was right saying that my solution was not the best possible! --- Prof Brian Ripley [EMAIL PROTECTED] wrote: read.table('temp.txt', check.names = FALSE) would be easier (and more general, since make.names can do more than prepend an 'X'). On Wed, 18 Jul 2007, Moshe Olshansky wrote: Hi Suzanne, My solution (which I am sure is not the best) would be: heat - read.table('temp.txt') heat X1905 X1910 X1950 X1992 X2011 X2020 Gnat 0.08 0.29 0.29 0.37 0.39 0.43 Snake 0.16 0.34 0.32 0.40 0.41 0.53 Bat0.40 0.54 0.52 0.60 0.60 0.63 Cat0.16 0.27 0.29 0.39 0.37 0.41 Dog0.43 0.54 0.52 0.61 0.60 0.62 Lynx 0.50 0.57 0.54 0.59 0.50 0.59 a-names(heat) b-strsplit(a,split=X) w-unlist(b) w [1] 1905 1910 1950 1992 2011 2020 z - w[seq(2,length(w),by=2)] z [1] 1905 1910 1950 1992 2011 2020 names(heat) - z heat 1905 1910 1950 1992 2011 2020 Gnat 0.08 0.29 0.29 0.37 0.39 0.43 Snake 0.16 0.34 0.32 0.40 0.41 0.53 Bat 0.40 0.54 0.52 0.60 0.60 0.63 Cat 0.16 0.27 0.29 0.39 0.37 0.41 Dog 0.43 0.54 0.52 0.61 0.60 0.62 Lynx 0.50 0.57 0.54 0.59 0.50 0.59 Regards, Moshe. --- Suzanne Matthews [EMAIL PROTECTED] wrote: Hello All, I have a simple question based on how things are labeled on my heat map; particularly, there is this annoying X that appears before the numeric value of all the labels of my columns. Let's say I have the following silly data, stored in temp.txt 19051910195019922011 2020 Gnat0.080.290.290.370.39 0.43 Snake 0.160.340.320.400.41 0.53 Bat 0.400.540.520.600.60 0.63 Cat 0.160.270.290.390.37 0.41 Dog 0.430.540.520.610.60 0.62 Lynx0.500.570.540.590.5 0.59 I use the following commands to generate my heatmap: heat - read.table('temp.txt') x - as.matrix(heat) heatmap.2(x, keysize=1.2, dendrogram=none, trace=none, Colv = FALSE, main = Silly Data, labCol= NULL, margin=c(7,8)) This generates a very nice heatmap, but there is one thing I have an issue with: How do I get rid of the 'X' that seems to come automatically before my numeric column values? I just want those columns to be labeled 1905, 1910, 1950, and so on. I cannot find anything in the heatmap.2 documentation that suggests how I should do this. Thank you very much for your time, and patience in reading this! Sincerely, Suzanne [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch 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] Matrix Multiplication, Floating-Point, etc.
After multiplication by 10 you get 6*8 = 48 - the result is an exact machine number so there is no roundoff, while 0.6*0.8 = 0.48, where neither of the 3 numbers (0.6, 0.8, 0.48) is an exact machine mumber. However, (-0.6)*0.8 should be equal EXACTLY to -(0.6*0.8), and in fact you get that sum(ev1*ev2) is exactly 0. What is strange is that you are not getting this result from ev1 %*% ev2. This means that either %^% uses some non-straightforward algorithm or it somehow sets the rounding control to something different from round to nearest. In the later case (-0.6) does not necessarily equal to -(0.6) and the rounding after multiplication is not necessarily symetric. Regards, Moshe. --- Talbot Katz [EMAIL PROTECTED] wrote: Thank you for responding! I realize that floating point operations are often inexact, and indeed, the difference between the two answers is within the all.equal tolerance, as mentioned in FAQ 7.31 (cited by Charles): (as.numeric(ev1%*%ev2))==(sum(ev1*ev2)) [1] FALSE all.equal((as.numeric(ev1%*%ev2)),(sum(ev1*ev2))) [1] TRUE I suppose that's good enough for numerical computation. But I was still surprised to see that matrix multiplication (ev1%*%ev2) doesn't give the exact right answer, whereas sum(ev1*ev2) does give the exact answer. I would've expected them to perform the same two multiplications and one addition. But I guess that's not the case. However, I did find that if I multiplied the two vectors by 10, making the entries integers (although the class was still numeric rather than integer), both computations gave equal answers of 0: xf1-10*ev1 xf2-10*ev2 (as.numeric(xf1%*%xf2))==(sum(xf1*xf2)) [1] TRUE Perhaps the moral of the story is that one should exercise caution and keep track of significant digits. -- TMK -- 212-460-5430 home 917-656-5351 cell From: Charles C. Berry [EMAIL PROTECTED] To: Talbot Katz [EMAIL PROTECTED] CC: r-help@stat.math.ethz.ch Subject: Re: [R] Matrix Multiplication, Floating-Point, etc. Date: Mon, 30 Jul 2007 09:27:42 -0700 7.31 Why doesn't R think these numbers are equal? On Fri, 27 Jul 2007, Talbot Katz wrote: Hi. I recently tried the following in R 2.5.1 on Windows XP: ev2-c(0.8,-0.6) ev1-c(0.6,0.8) ev1%*%ev2 [,1] [1,] -2.664427e-17 sum(ev1*ev2) [1] 0 (I got the same result with R 2.4.1 on a different Windows XP machine.) I expect this issue is very familiar and probably has been discussed in this forum before. Can someone please point me to some documentation or discussion about this? Is there some standard way to get the correct answer from %*%? Thanks! -- TMK -- 212-460-5430home 917-656-5351cell __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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 plot a differential equation?
Hi Montserrat, What exactly would you like to plot? Your differential equation can be easily integrated so that you can get an implicit expression for F(x), i.e. expression like G(c,x,F(x)) = 0 where G is a known function and c is an arbitrary constant. For every value of c and each value of x such that x (-1/k2 + (k1+k2)*ln(k1+k2)/k2^2) [if my calculations are correct!] there exist 2 possible values of F(x): one with F(x) 1 and one with F(x) 1. So for each c you have two branches of the function F defined for x x(c). You will need a numerical procedure to find these two values of F(x). Regards, Moshe. --- Montse Rue [EMAIL PROTECTED] wrote: Hi, I would like to plot the following equation: dF(x)/dx=(k1+k2F(x))(1-F(x)) where k1 and k2 are parameters that I have estimated already. How can I plot the curve in R? Thanks! Montserrat Rue Universitat de Lleida (Spain) [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Solve equations
Hi Sebastian, Your equations can be easily solved - no programming is required! Let's number you equations: (1) 0.007= 2VZ (2) 0.03= W(Y+Z) (3) 0.034= X(y+Z) (4) 0.013 = (X+W)Y +(X-W)Z (5) X = W+V Substitute (5) into (3) and then divide (2) by (3) to get (W+V)/W = 0.034/0.03 so that (6) V = (2/15)*W (7) X = V+W = (17/15)*W Substitute (7) into (3) and (4) to get: (3') (17/15)*(WY + WZ) = 0.034 (4') (32/15)WY + (2/15)WZ = 0.013 Solving (3') and (4') for two new unknowns WY and WZ yields (8) WY = 0.0045 (9) WZ = 0.0255 So that (10) Z = (17/3)*Y Substitute (10) into (1) to get 2*Y*(17/3)*Y = 0.07 So (since Y 0), Y = sqrt(0.21/34) = 0.07859052 From (10), Z = (17/3)*Y = 0.4453463 From (8), W = 0.0045/Y = 0.05725881 From (6), V = (2/15)*W = 0.007634508 Finally, X = V + W = 0.06489332 So there is a unique solution (which luckily satisfies your constraints!). Regards, Moshe. --- sebastien puechmaille [EMAIL PROTECTED] wrote: Hello, I have a system of five equations to solve with five 'unknows'(V, W, X, Y and Z) and constraints. The equations are: 0.007= 2VZ 0.03= W(Y+Z) 0.034= X(y+Z) 0.013 = (X+W)Y +(X-W)Z X = W+V Constraints: 0VWX 0YZ1 Does anyone know a R-package to solve this system? Thanks, E-mail: [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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 a null expectation
Hi Ryan, I think that the only function you need is the one which generates random variable. The simplest one is to generate a uniform (0,1) variable U and then do not move if U 0.25, move to the first nearest point if 0.25 = U 0.5, move to the second nearest point if ).5 = U 0.75 and move to the third nearest point if U = 0.75. To make things work faster you may wish to prepare for which of your points the indexes of the 3 nearest points (i.e. if you have N points you prepare an Nx3 matrix of indexes). Regards, Moshe. --- Ryan Briscoe Runquist [EMAIL PROTECTED] wrote: Hello, I am trying to write a simulation to generate a null hypothesis so that I can test my data. Essentially we started out with a grid of points and I have a symmetrical matrix of the distance between all of these points. I then want to simulate random movement. So I want to be able to start at a point, search through the matrix for the three closest points and with equal probability move to one of them or stay at the same place and then record the distance moved as well as the relative elevation (a vector of data that I also have) and then starting at that grid-point repeat the process a number of times. Does anyone know of a function that could help me with this? Thanks, Ryan ~~ Ryan D. Briscoe Runquist Population Biology Graduate Group University of California, Davis [EMAIL PROTECTED] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Spatial sampling problem
Hi SK, If I understand you correctly you are trying to assign to each point of a 10x10 2D grid a normal variable with a certain dependency between the values at each grid point. Is this correct? If so, please explain the dependence structure. Do you need something like X(i,j) = f(i,j) + Z(i,j) for 1=i,j=10 where f is a deterministic function and Z(i,j) are independent normal variables? Or do you want to assume that Z(i,j) have a multivariate normal distribution (and can specify the covariance matrix)? Regards, Moshe. --- sk [EMAIL PROTECTED] wrote: Hi All, I am new in R and trying to simulate random normal 2D field with mean trend say north-south. My domain is 10x10 grid and I am trying to use mvnorm but do not know how to specify the domain and the mean field. I would appreciate any help. Cheers, SK - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] small sample techniques
Hi Nair, If the two populations are normal the t-test gives you the exact result for whatever the sample size is (the sample size will affect the number of degrees of freedom). When the populations are not normal and the sample size is large it is still OK to use t-test (because of the Central Limit Theorem) but this is not necessarily true for the small sample size. You could use simulation to find the relevant probabilities. --- Nair, Murlidharan T [EMAIL PROTECTED] wrote: If my sample size is small is there a particular switch option that I need to use with t.test so that it calculates the t ratio correctly? Here is a dummy example? á =0.05 Mean pain reduction for A =27; B =31 and SD are SDA=9 SDB=12 drgA.p-rnorm(5,27,9); drgB.p-rnorm(5,31,12) t.test(drgA.p,drgB.p) # what do I need to give as additional parameter here? I can do it manually but was looking for a switch option that I need to specify for t.test. Thanks ../Murli [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] XL files
Hi Animesh, Can you send an example of an Excel file you need to process (and the result you wish to get)? Regards, Moshe. --- Acharjee, Animesh [EMAIL PROTECTED] wrote: Dear All, I am new to R. I need to read XLs files, parse the data and convert them into txt format for further calculation . If any one have code for that , please send me code or suggestions for the same.Waiting for reply. Thanking you Animesh [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] small sample techniques
As Thomas Lumley noted, there exist several versions of t-test. If you use t1 - t.test(x,y) then no assumption is made of x and y having equal variance and of the two sample sizes being equal and then an approximate t-test is used with an approximate number of degrees of freedom (and this is what you got). If you use t2 - t.test(x,y,var.equal=TRUE) then equal variance is assumed and you get 8 degrees of freedom. If you use t3 - t.test(x,y,paired=TRUE) then equal sample sizes are assumed and the number of degrees of freedom is 4 (5-1). --- Nair, Murlidharan T [EMAIL PROTECTED] wrote: Indeed, I understand what you say. The df of freedom for the dummy example is n1+n2-2 = 8. But when I run the t.test I get it as 5.08, am I missing something? -Original Message- From: Moshe Olshansky [mailto:[EMAIL PROTECTED] Sent: Tuesday, August 07, 2007 9:05 PM To: Nair, Murlidharan T; r-help@stat.math.ethz.ch Subject: Re: [R] small sample techniques Hi Nair, If the two populations are normal the t-test gives you the exact result for whatever the sample size is (the sample size will affect the number of degrees of freedom). When the populations are not normal and the sample size is large it is still OK to use t-test (because of the Central Limit Theorem) but this is not necessarily true for the small sample size. You could use simulation to find the relevant probabilities. --- Nair, Murlidharan T [EMAIL PROTECTED] wrote: If my sample size is small is there a particular switch option that I need to use with t.test so that it calculates the t ratio correctly? Here is a dummy example? á =0.05 Mean pain reduction for A =27; B =31 and SD are SDA=9 SDB=12 drgA.p-rnorm(5,27,9); drgB.p-rnorm(5,31,12) t.test(drgA.p,drgB.p) # what do I need to give as additional parameter here? I can do it manually but was looking for a switch option that I need to specify for t.test. Thanks ../Murli [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] small sample techniques
Well, this an explanation of what is done in the paired t-test (and why the number of df is as it is). I was too lazy to write all this. It is nice that some list members are less lazy! --- Rolf Turner [EMAIL PROTECTED] wrote: On 9/08/2007, at 2:57 PM, Moshe Olshansky wrote: As Thomas Lumley noted, there exist several versions of t-test. snip If you use t3 - t.test(x,y,paired=TRUE) then equal sample sizes are assumed and the number of degrees of freedom is 4 (5-1). This is seriously misleading. The assumption is not that the sample sizes are equal, but rather that there is ***just one sample***, namely the sample of differences. More explicitly the assumptions are that x_i - y_i are i.i.d. Gaussian with mean mu and variance sigma^2. One is trying to conduct inference about mu, of course. It should also be noted that it is a crucial assumption for the ``non-paired'' t-test that the two samples be ***independent*** of each other, as well as being Gaussian. None of this is however germane to Nair's original question; it is clear that he is interested in a two-independent-sample t-test. cheers, Rolf Turner ## Attention: This e-mail message is privileged and confidential. If you are not the intended recipient please delete the message and notify the sender. Any views or opinions presented are solely those of the author. This e-mail has been scanned and cleared by MailMarshal www.marshalsoftware.com ## __ R-help@stat.math.ethz.ch 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] small sample techniques
Hi Murli, First of all, regarding prop.test, you made a typo: you should have used prop.test(c(69,90),c(300,300)) which gives you the squared value of 3.4228, and it's square root is 1.85 which is not too far from 1.94. I would use Fisher Exact Test (fisher.test). Two sided test has a p-value of 0.06411 so you do not reject H0, One sided test (i.e. H1 is that the first probability of success is smaller than the second) has a p-value of 0.03206, so you reject H0 (with 95% confidence level). You get similar results with two-sided and one-sided t-test. Moshe. P.S. if you use paired t-test you get nonsense since it uses pairwise differences, and in your case only 21 of 300 differences are non-zero! --- Nair, Murlidharan T [EMAIL PROTECTED] wrote: n=300 30% taking A relief from pain 23% taking B relief from pain Question; If there is no difference are we likely to get a 7% difference? Hypothesis H0: p1-p2=0 H1: p1-p2!=0 (not equal to) 1Weighed average of two sample proportion 300(0.30)+300(0.23) --- = 0.265 300+300 2Std Error estimate of the difference between two independent proportions sqrt((0.265 *0.735)*((1/300)+(1/300))) = 0.03603 3Evaluation of the difference between sample proportion as a deviation from the hypothesized difference of zero ((0.30-0.23)-(0))/0.03603 = 1.94 z did not approach 1.96 hence H0 is not rejected. This is what I was trying to do using prop.test. prop.test(c(30,23),c(300,300)) What function should I use? -Original Message- From: [EMAIL PROTECTED] on behalf of Nordlund, Dan (DSHS/RDA) Sent: Thu 8/9/2007 1:26 PM To: r-help@stat.math.ethz.ch Subject: Re: [R] small sample techniques -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Nair, Murlidharan T Sent: Thursday, August 09, 2007 9:19 AM To: Moshe Olshansky; Rolf Turner; r-help@stat.math.ethz.ch Subject: Re: [R] small sample techniques Thanks, that discussion was helpful. Well, I have another question I am comparing two proportions for its deviation from the hypothesized difference of zero. My manually calculated z ratio is 1.94. But, when I calculate it using prop.test, it uses Pearson's chi-squared test and the X-squared value that it gives it 0.74. Is there a function in R where I can calculate the z ratio? Which is ('p1-'p2)-(p1-p2) Z= S ('p1-'p2) Where S is the standard error estimate of the difference between two independent proportions Dummy example This is how I use it prop.test(c(30,23),c(300,300)) Cheers../Murli Murli, I think you need to recheck you computations. You can run a t-test on your data in a variety of ways. Here is one: x-c(rep(1,30),rep(0,270)) y-c(rep(1,23),rep(0,277)) t.test(x,y) Welch Two Sample t-test data: x and y t = 1.0062, df = 589.583, p-value = 0.3147 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.02221086 0.06887752 sample estimates: mean of x mean of y 0.1000 0.0767 Hope this is helpful, Dan Daniel J. Nordlund Research and Data Analysis Washington State Department of Social and Health Services Olympia, WA 98504-5204 __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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] binomial simulation
As I understand this, P(T+ | D-)=1-P(T+ | D+)=0.05 is the probability not to detect desease for a person at ICU who has the desease. Correct? What I asked was whether it is possible to mistakenly detect the desease for a person who does not have it? Assuming that this is impossible the formula is below: If there are N patients, each has a probability p to have the desease (p=0.6 in your case) and q is the probability to detect the desease for a person who has it (q = 0.95 for ICU and q = 0.8 for a regular unit), then P(k have the desease AND m are detected) = P(k have the desease)*P(m are detected / k have the desease) = C(N,k)*p^k*(1-p)^(N-k)*C(k,m)*q^m*(1-q)^(k-m) where C(a,b) is the Binomial coefficient a above b - the number of ways to choose b items out of a (when the order does not matter). You of course must assume that N = k = m = 0 (otherwise the probability is 0). To generate such pairs (k infected and m detected) you can do the following: k - rbinom(N,1,p) m - rbinom(k,1,q) Regards, Moshe. --- sigalit mangut-leiba [EMAIL PROTECTED] wrote: Hi, The probability of false detection is: P(T+ | D-)=1-P(T+ | D+)=0.05. and I want to find the joint probability P(T+,D+)=P(T+|D+)*P(D+) Thank you for your reply, Sigalit. On 8/13/07, Moshe Olshansky [EMAIL PROTECTED] wrote: Hi Sigalit, Do you want to find the probability P(T+ = t AND D+ = d) for all the combinations of t and d (for ICU and Reg.)? Is the probability of false detection (when there is no disease) always 0? Regards, Moshe. --- sigalit mangut-leiba [EMAIL PROTECTED] wrote: hello, I asked about this simulation a few days ago, but still i can't get what i need. I have 2 units: icu and regular. from icu I want to take 200 observations from binomial distribution, when probability for disease is: p=0.6. from regular I want to take 300 observation with the same probability: p=0.6 . the distribution to detect disease when disease occurred- *for someone from icu* - is: p(T+ | D+)=0.95. the distribution to detect disease when disease occurred- *for someone from reg.unit* - is: p(T+ | D+)=0.8. I want to compute the joint distribution for each unit: p(T+,D+) for icu, and the same for reg. I tried: pdeti - 0 pdetr - 0 picu - pdeti*.6 preg - pdetr*.6 dept - c(icu,reg) icu - rbinom(200, 1, .6) reg - rbinom(300, 1, .6) for(i in 1:300) { if(dept==icu) pdeti==0.95 if (dept==reg) pdetr==0.80 } print(picu) print(preg) and got 50 warnings: the condition has length 1 and only the first element will be used in: if (dept == icu) pdeti == 0.95 the condition has length 1 and only the first element will be used in: if (dept == reg) pdetr == 0.8 I would appreciate any suggestions, thank you, Sigalit. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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] invert 160000x160000 matrix
While inverting the matrix may be a problem, if you need to solve an equation A*x = b you do not need to invert A, there exist iterative methods which do need A or inv(A) - all you need to provide is a function that computes A*x for an arbitrary vector x. For such a large matrix this may be slow but possible. --- Paul Gilbert [EMAIL PROTECTED] wrote: I don't think you can define a matrix this large in R, even if you have the memory. Then, of course, inverting it there may be other programs that have limitations. Paul Jiao Yang wrote: Can R invert a 16x16 matrix with all positive numbers? Thanks a lot! __ R-help@stat.math.ethz.ch 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. La version française suit le texte anglais. This email may contain privileged and/or confidential inform...{{dropped}} __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] binomial simulation
No wonder that you are getting overflow, since gamma(N+1) = n! and 200! (200/e)^200 10^370. There exists another way to compute C(N,k). Let me know if you need this and I will explain to you how this can be done. But do you really need to compute the individual probabilities? May be you need something else and there is no need to compute the individual probabilities? Regards, Moshe. --- sigalit mangut-leiba [EMAIL PROTECTED] wrote: Thank you, I'm trying to run the joint probabilty: C(N,k)*p^k*(1-p)^(N-k)*C(k,m)*q^m*(1-q)^(k-m) and get the error: Error in C(N, k) : object not interpretable as a factor so I tried the long way: gamma(N+1)/(gamma(k+1)*(gamma(N-k))) and the same with k, and got the error: 1: value out of range in 'gammafn' in: gamma(N + 1) 2: value out of range in 'gammafn' in: gamma(N - k) Do you know why it's not working? Thanks again, Sigalit. On 8/14/07, Moshe Olshansky [EMAIL PROTECTED] wrote: As I understand this, P(T+ | D-)=1-P(T+ | D+)=0.05 is the probability not to detect desease for a person at ICU who has the desease. Correct? What I asked was whether it is possible to mistakenly detect the desease for a person who does not have it? Assuming that this is impossible the formula is below: If there are N patients, each has a probability p to have the desease (p=0.6 in your case) and q is the probability to detect the desease for a person who has it (q = 0.95 for ICU and q = 0.8 for a regular unit), then P(k have the desease AND m are detected) = P(k have the desease)*P(m are detected / k have the desease) = C(N,k)*p^k*(1-p)^(N-k)*C(k,m)*q^m*(1-q)^(k-m) where C(a,b) is the Binomial coefficient a above b - the number of ways to choose b items out of a (when the order does not matter). You of course must assume that N = k = m = 0 (otherwise the probability is 0). To generate such pairs (k infected and m detected) you can do the following: k - rbinom(N,1,p) m - rbinom(k,1,q) Regards, Moshe. --- sigalit mangut-leiba [EMAIL PROTECTED] wrote: Hi, The probability of false detection is: P(T+ | D-)=1-P(T+ | D+)=0.05. and I want to find the joint probability P(T+,D+)=P(T+|D+)*P(D+) Thank you for your reply, Sigalit. On 8/13/07, Moshe Olshansky [EMAIL PROTECTED] wrote: Hi Sigalit, Do you want to find the probability P(T+ = t AND D+ = d) for all the combinations of t and d (for ICU and Reg.)? Is the probability of false detection (when there is no disease) always 0? Regards, Moshe. --- sigalit mangut-leiba [EMAIL PROTECTED] wrote: hello, I asked about this simulation a few days ago, but still i can't get what i need. I have 2 units: icu and regular. from icu I want to take 200 observations from binomial distribution, when probability for disease is: p=0.6. from regular I want to take 300 observation with the same probability: p=0.6 . the distribution to detect disease when disease occurred- *for someone from icu* - is: p(T+ | D+)=0.95. the distribution to detect disease when disease occurred- *for someone from reg.unit* - is: p(T+ | D+)=0.8. I want to compute the joint distribution for each unit: p(T+,D+) for icu, and the same for reg. I tried: pdeti - 0 pdetr - 0 picu - pdeti*.6 preg - pdetr*.6 dept - c(icu,reg) icu - rbinom(200, 1, .6) reg - rbinom(300, 1, .6) for(i in 1:300) { if(dept==icu) pdeti==0.95 if (dept==reg) pdetr== 0.80 } print(picu) print(preg) and got 50 warnings: the condition has length 1 and only the first element will be used in: if (dept == icu) pdeti == 0.95 the condition has length 1 and only the first element will be used in: if (dept == reg) pdetr == 0.8 I would appreciate any suggestions, thank you, Sigalit. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.htmlhttp://www.r-project.org/posting-guide.html and provide commented, minimal, self-contained, === message truncated === __ R-help@stat.math.ethz.ch 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 write to a table column by column?
Hi Yuchen, First of all please notice that you may not have more than 2^8 = 256 columns in Excel, so if you have more than 256 time-series you can not put them all in one Excel sheet. I also believe that you can not have more than 2^16 = 65536 rows. If you do not have more than 256 time-series, I see two possibilities. The first one is to create a list which will hold all your time-series (each list member holds one time-series) and then make a double loop: one from 1 to as long as there is more data and if the first variable is i you create a string by looping through the list and adding to your string either element i (and a tab) of time-series j if it's length is = j or a tab (and a tab) otherwise. Then you can write each row to the file (using 'cat'). When you are done import the resulting file into Excel. Another possibility is to use xlsReadWritePro which allows you to write to a specific column of a specific sheet of an Excel file (and to append this). Regards, Moshe. --- Yuchen Luo [EMAIL PROTECTED] wrote: Dear friends. Every loop of my program will result in a list that is very long, with a structure similar to the one below: Lst - list(name=Fred, wife=Mary, daily.incomes=c(1:850)) Please notice the large size of daily.incomes. I need to store all such lists in a csv file so that I can easily view them in Excel. Excel cannot display a row of more than 300 elements, therefore, I have to store the lists as columns. It is not hard to store one list as a column in the csv file. The problem is how to store the second list as a second column, so that the two columns will lie side by side to each other and I can easily compare their elements. ( If I use 'appened=TRUE', the second time series will be stored in the same column. ) Thank you for your tine and your help will be highly appreciated!! Best Yuchen Luo [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] binomial simulation
Thank you - I wasn't aware of this function. One can even use lchoose which allows really huge arguments (more than 2^1000)! --- Lucke, Joseph F [EMAIL PROTECTED] wrote: C is an R function for setting contrasts in a factor. Hence the funky error message. ?C Use choose() for your C(N,k) ?choose choose(200,2) 19900 choose(200,100) 9.054851e+58 N=200; k=100; m=50; p=.6; q=.95 choose(N,k)*p^k*(1-p)^(N-k)*choose(k,m)*q^m*(1-q)^(k-m) 6.554505e-41 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Moshe Olshansky Sent: Wednesday, August 15, 2007 2:06 AM To: sigalit mangut-leiba; r-help Subject: Re: [R] binomial simulation No wonder that you are getting overflow, since gamma(N+1) = n! and 200! (200/e)^200 10^370. There exists another way to compute C(N,k). Let me know if you need this and I will explain to you how this can be done. But do you really need to compute the individual probabilities? May be you need something else and there is no need to compute the individual probabilities? Regards, Moshe. --- sigalit mangut-leiba [EMAIL PROTECTED] wrote: Thank you, I'm trying to run the joint probabilty: C(N,k)*p^k*(1-p)^(N-k)*C(k,m)*q^m*(1-q)^(k-m) and get the error: Error in C(N, k) : object not interpretable as a factor so I tried the long way: gamma(N+1)/(gamma(k+1)*(gamma(N-k))) and the same with k, and got the error: 1: value out of range in 'gammafn' in: gamma(N + 1) 2: value out of range in 'gammafn' in: gamma(N - k) Do you know why it's not working? Thanks again, Sigalit. On 8/14/07, Moshe Olshansky [EMAIL PROTECTED] wrote: As I understand this, P(T+ | D-)=1-P(T+ | D+)=0.05 is the probability not to detect desease for a person at ICU who has the desease. Correct? What I asked was whether it is possible to mistakenly detect the desease for a person who does not have it? Assuming that this is impossible the formula is below: If there are N patients, each has a probability p to have the desease (p=0.6 in your case) and q is the probability to detect the desease for a person who has it (q = 0.95 for ICU and q = 0.8 for a regular unit), then P(k have the desease AND m are detected) = P(k have the desease)*P(m are detected / k have the desease) = C(N,k)*p^k*(1-p)^(N-k)*C(k,m)*q^m*(1-q)^(k-m) where C(a,b) is the Binomial coefficient a above b - the number of ways to choose b items out of a (when the order does not matter). You of course must assume that N = k = m = 0 (otherwise the probability is 0). To generate such pairs (k infected and m detected) you can do the following: k - rbinom(N,1,p) m - rbinom(k,1,q) Regards, Moshe. --- sigalit mangut-leiba [EMAIL PROTECTED] wrote: Hi, The probability of false detection is: P(T+ | D-)=1-P(T+ | D+)=0.05. and I want to find the joint probability P(T+,D+)=P(T+|D+)*P(D+) Thank you for your reply, Sigalit. On 8/13/07, Moshe Olshansky [EMAIL PROTECTED] wrote: Hi Sigalit, Do you want to find the probability P(T+ = t AND D+ = d) for all the combinations of t and d (for ICU and Reg.)? Is the probability of false detection (when there is no disease) always 0? Regards, Moshe. --- sigalit mangut-leiba [EMAIL PROTECTED] wrote: hello, I asked about this simulation a few days ago, but still i can't get what i need. I have 2 units: icu and regular. from icu I want to take 200 observations from binomial distribution, when probability for disease is: p=0.6. from regular I want to take 300 observation with the same probability: p=0.6 . the distribution to detect disease when disease occurred- *for someone from icu* - is: p(T+ | D+)=0.95. the distribution to detect disease when disease occurred- *for someone from reg.unit* - is: p(T+ | D+)=0.8. I want to compute the joint distribution for each === message truncated === __ R-help@stat.math.ethz.ch 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] Combine matrix
This does not work in the general case. To see this, try: rownames(a}[5] - a and see what happens then. --- François Pinard [EMAIL PROTECTED] wrote: [Gianni Burgin] let say something like this a=matrix(1:25, nrow=5) rownames(a)=letters[1:5] colnames(a)=rep(A, 5) a A A A A A a 1 6 11 16 21 b 2 7 12 17 22 c 3 8 13 18 23 d 4 9 14 19 24 e 5 10 15 20 25 b=matrix(1:40, nrow=8) rownames(b)=c(rep(a,4),rep(b,4)) colnames(b)=rep(B, 5) b B B B B B a 1 9 17 25 33 a 2 10 18 26 34 a 3 11 19 27 35 a 4 12 20 28 36 b 5 13 21 29 37 b 6 14 22 30 38 b 7 15 23 31 39 b 8 16 24 32 40 as a results I wold like something like A A A A A B B B B B a 1 6 11 16 21 1 9 17 25 33 a 1 6 11 16 21 2 10 18 26 34 a 1 6 11 16 21 3 11 19 27 35 a 1 6 11 16 21 4 12 20 28 36 b 2 7 12 17 22 5 13 21 29 37 b 2 7 12 17 22 6 14 22 30 38 b 2 7 12 17 22 7 15 23 31 39 b 2 7 12 17 22 8 16 24 32 40 does it is clear? is there a function that automate this operation? Like, maybe: cbind(a[rownames(b),], b) -- François Pinard http://pinard.progiciels-bpi.ca __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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 to find coodinates in an array
A not very good solution is as below: If your array's dimensions were KxMxN and the linear index is i then n - ceiling(i/(K*M)) i1 - i - (n-1)*(K*M) m - ceiling(i1/K) k - i1 - (m-1)*K and your index is (k,m,n) I am almost sure that there is a function in R which does this (it exists in Matlab). Regards, Moshe. --- Ana Conesa [EMAIL PROTECTED] wrote: Dear list, I am looking for a function/way to get the array coordinates of given elements in an array. What I mean is the following: - Let X be a 3D array - I find the ordering of the elements of X by ord - order(X) (this returns me a vector) - I now want to find the x,y,z coordinates of each element of ord Can anyone help me? Thanks! Ana __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Fwd: Re: Combine matrix
Note: forwarded message attached. Mistakenly I did not send this to the list! ---BeginMessage--- Here is a solution which (I think) works, but it contains a loop which most likely can be eliminated by more knowledgable people): Let a and b be like in your example but to make things more interesting do the following: rownames(a)[3:4]-c(b,a) and now: a1 - a[rownames(a) %in% rownames(b),] b1 - b[rownames(b) %in% rownames(a),] c - unique(rownames(a1)) ia - numeric(0) ib - numeric(0) for (ronam in c) { i - which(rownames(a1) == ronam) j - which(rownames(b1) == ronam) i1 - rep(i,length(j)) ia - c(ia,t(matrix(i1,nrow=length(i ib - c(ib,rep(j,length(i))) } result - cbind(a1[ia,],b1[ib,]) --- Gianni Burgin [EMAIL PROTECTED] wrote: let say something like this a=matrix(1:25, nrow=5) rownames(a)=letters[1:5] colnames(a)=rep(A, 5) a A A A A A a 1 6 11 16 21 b 2 7 12 17 22 c 3 8 13 18 23 d 4 9 14 19 24 e 5 10 15 20 25 b=matrix(1:40, nrow=8) rownames(b)=c(rep(a,4),rep(b,4)) colnames(b)=rep(B, 5) b B B B B B a 1 9 17 25 33 a 2 10 18 26 34 a 3 11 19 27 35 a 4 12 20 28 36 b 5 13 21 29 37 b 6 14 22 30 38 b 7 15 23 31 39 b 8 16 24 32 40 as a results I wold like something like A A A A A B B B B B a 1 6 11 16 21 1 9 17 25 33 a 1 6 11 16 21 2 10 18 26 34 a 1 6 11 16 21 3 11 19 27 35 a 1 6 11 16 21 4 12 20 28 36 b 2 7 12 17 22 5 13 21 29 37 b 2 7 12 17 22 6 14 22 30 38 b 2 7 12 17 22 7 15 23 31 39 b 2 7 12 17 22 8 16 24 32 40 does it is clear? is there a function that automate this operation? thank you very much! On 8/16/07, jim holtman [EMAIL PROTECTED] wrote: Can you provide an example of what you mean; e.g., the two input matrices and the desired output. On 8/16/07, Gianni Burgin [EMAIL PROTECTED] wrote: Hi R user, I am new to R, and I have a very simple question for you. I have two matrix A and B, with internally redundant rownames (but variables are different). Some, but not all the rownames are shared among the two matrix. I want to create a greater matrix that combines the previuos two, and has all the possible combinations of matching rownames lines among matrix A and B. looking for the solution I bumped in merge but actually works on data.frame, and in dataframe there could be no redundancy in names. can you help me?? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem you are trying to solve? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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. ---End Message--- __ R-help@stat.math.ethz.ch 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] Partial comparison in string vector
Use grep(^e,x) if you are looking for the entries where e is the first character. --- Vladimir Eremeev [EMAIL PROTECTED] wrote: Hi! seq(along=x) %in% grep(e,x) Steve Powell-4 wrote: I have a vector of strings x=c(w,ex,ee) And I want to get a logical vector showing the positions where my search string e matches the elements partially, i.e. is at least the left-hand part of the target strings, i.e. I want to get a vector FALSE TRUE TRUE. -- View this message in context: http://www.nabble.com/Partial-comparison-in-string-vector-tf4304145.html#a12251593 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Calculating diameters of cirkels in a picture.
Hi Bart, I have never used image processing software in R (I was doing this with Matlab), but here is what I would have done algorithmically: 1) convert the picture to gray-scale 2) find a threshold value which separates the circles from the background and convert your image to black and white 3) if the circles are far apart use morphological closing to fill in small holes inside the circles (may be do this several times) 4) use labeling to split the image into connected components 5) for each connected component get it's area (the number of pixels) and use the formula S = Pi*R^2 to find the approximate radii. Regards, Moshe. --- Julian Burgos [EMAIL PROTECTED] wrote: Hi Bart, If you only have 36 circles, the fastest way would be to use some image processing software and measure the circles by hand. One option is to use ImageJ, which you can download here http://rsb.info.nih.gov/ij/ Julian Bart Joosen wrote: Hi, Maybe this is more a programming questions than a specific R-project question, but maybe there is someone who can point me in the right direction. I have a picture of cirkels which I took with a digital camera. Now I want to use the diameter of the cirkels on the picture for analysis in R. I can use pixmap to import the picture, but how do I find the outside cirkels and calculate the diameter? I pointed out that I can use the edci package, but then I need to preprocess the data to reduce the points, otherwise it takes a long time, and my computer crashes. If you want to see such a picture, I cropped a larger one, and highlighted the cirkel which is of interest. In a real world, this is a plate with 36 cirkels, which all should be measured. www.users.skynet.be/fa244930/fotos/outlined.jpg Thanks for your time Bart [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Choice of data.frame column by index?
This won't work since it produces a matrix (try this). What should work is x[(1:nrow(x)) + nrow(x)*(v-1)] --- Johannes Graumann [EMAIL PROTECTED] wrote: Thanks! Joh On Thursday 23 August 2007 12:01:50 you wrote: x[cbind(1:nrow(x), the.vector)] Patrick Burns [EMAIL PROTECTED] +44 (0)20 8525 0696 http://www.burns-stat.com (home of S Poetry and A Guide for the Unwilling S User) Johannes Graumann wrote: Hello, Imagine a data frame like so: Intensity0 Intensity1 1 767432.1 451743.4 2 3998988.0 4642145.0 3 818974.6 552315.8 and a vector like so: [1] 1 2 1 How can I get R to produce a vector that contains the value in one column or the other depending on the vector? The result should look like [1] 767432.1 4642145.0 818974.6 Thanks for any hints! Joh __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Choice of data.frame column by index?
Correct, I didn't notice that the coma was inside the cbind(). Sorry... --- Rolf Turner [EMAIL PROTECTED] wrote: On 24/08/2007, at 12:51 PM, Moshe Olshansky wrote: This won't work since it produces a matrix (try this). On the contrary, Patrick's solution is correct. I tried it. It works just fine. cheers, Rolf Turner What should work is x[(1:nrow(x)) + nrow(x)*(v-1)] --- Johannes Graumann [EMAIL PROTECTED] wrote: Thanks! Joh On Thursday 23 August 2007 12:01:50 you wrote: x[cbind(1:nrow(x), the.vector)] Patrick Burns [EMAIL PROTECTED] +44 (0)20 8525 0696 http://www.burns-stat.com (home of S Poetry and A Guide for the Unwilling S User) Johannes Graumann wrote: Hello, Imagine a data frame like so: Intensity0 Intensity1 1 767432.1 451743.4 2 3998988.0 4642145.0 3 818974.6 552315.8 and a vector like so: [1] 1 2 1 How can I get R to produce a vector that contains the value in one column or the other depending on the vector? The result should look like [1] 767432.1 4642145.0 818974.6 Thanks for any hints! Joh __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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. ## Attention: This e-mail message is privileged and confidential. If you are not the intended recipient please delete the message and notify the sender. Any views or opinions presented are solely those of the author. This e-mail has been scanned and cleared by MailMarshal www.marshalsoftware.com ## __ R-help@stat.math.ethz.ch 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] Calculating diameters of cirkels in a picture.
Hi Bart, One more comment: You do not really need the morphological closing to close the holes inside the circles. Another possibility is to reverse the black-and-withe picture, i.e. make the holes and background be 1 and the circles 0, label the connected components and then only the component which touches the boundaries is the background while all other components are holes and you can make them white (1) in the original black-and-white image. --- Moshe Olshansky [EMAIL PROTECTED] wrote: Hi Bart, I have never used image processing software in R (I was doing this with Matlab), but here is what I would have done algorithmically: 1) convert the picture to gray-scale 2) find a threshold value which separates the circles from the background and convert your image to black and white 3) if the circles are far apart use morphological closing to fill in small holes inside the circles (may be do this several times) 4) use labeling to split the image into connected components 5) for each connected component get it's area (the number of pixels) and use the formula S = Pi*R^2 to find the approximate radii. Regards, Moshe. --- Julian Burgos [EMAIL PROTECTED] wrote: Hi Bart, If you only have 36 circles, the fastest way would be to use some image processing software and measure the circles by hand. One option is to use ImageJ, which you can download here http://rsb.info.nih.gov/ij/ Julian Bart Joosen wrote: Hi, Maybe this is more a programming questions than a specific R-project question, but maybe there is someone who can point me in the right direction. I have a picture of cirkels which I took with a digital camera. Now I want to use the diameter of the cirkels on the picture for analysis in R. I can use pixmap to import the picture, but how do I find the outside cirkels and calculate the diameter? I pointed out that I can use the edci package, but then I need to preprocess the data to reduce the points, otherwise it takes a long time, and my computer crashes. If you want to see such a picture, I cropped a larger one, and highlighted the cirkel which is of interest. In a real world, this is a plate with 36 cirkels, which all should be measured. www.users.skynet.be/fa244930/fotos/outlined.jpg Thanks for your time Bart [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Calculating diameters of cirkels in a picture.
Hi Bart, Let's assume that you situation was simpler - you have a BW (Black and White) image containing circles (in white) and you need to find the diameter of each circle (and of course to know how many circles you have). This can be done with labeling of connected components. You say that two pixels are neighbors if they have common edge (4-connectivity) or at least a common vertex (8-connectivity). So now you can treat your image (white pixels) as a graph (with edges connecting any two neighbors). Then each connected component of that graph corresponds to a circle. There exists a well know algorithm to do this. It takes the original BW image (where every image pixel has the value of 1 and background pixel the value of 0) and produces an image where every background pixel still has the value of 0, every pixel of the first connected component has the value of 1, every pixel of the second connected component has the value of 2, etc. So no you can process each connected component (circle in your case) separately. Basically this is all you need. You can either count the number of pixels having the value of k to find the area (and then the diameter) or just take (maximal x value) - (minimal x value) + 1. In your case it can happen that after you convert your image into BW image some circles will have holes inside with some small objects inside these holes, and you do not want to consider these small objects as additional circles. So I thought of using morphological closing to get rid of small holes, but as I wrote in the following note you do not need this. When you get the BW image take the complimentary one (i.e. background pixels have the value of 1 and image pixels the value of 0). Label the connected components of the background. Only one of them is real background - all others are inside circles. Real background touches the image boundaries. Now go to the original BW image and give all the pixels outside the real background the value of 1. Now all your circles are full (no holes) and you can proceed as above. Best regards, Moshe. --- Bartjoosen [EMAIL PROTECTED] wrote: Hi All, I really like to thank you for the answers, while I was searching for some edge detection and clustering algorithms, Moshe came with a simple but effective solution: use the area to find the diameter! But I tried Moshe's solution, but I couldn't figure out what you mean with morphological closing and the labeling to split the images. Could you please clarify this a bit? Thanks for your support Bart Moshe Olshansky-2 wrote: Hi Bart, One more comment: You do not really need the morphological closing to close the holes inside the circles. Another possibility is to reverse the black-and-withe picture, i.e. make the holes and background be 1 and the circles 0, label the connected components and then only the component which touches the boundaries is the background while all other components are holes and you can make them white (1) in the original black-and-white image. --- Moshe Olshansky [EMAIL PROTECTED] wrote: Hi Bart, I have never used image processing software in R (I was doing this with Matlab), but here is what I would have done algorithmically: 1) convert the picture to gray-scale 2) find a threshold value which separates the circles from the background and convert your image to black and white 3) if the circles are far apart use morphological closing to fill in small holes inside the circles (may be do this several times) 4) use labeling to split the image into connected components 5) for each connected component get it's area (the number of pixels) and use the formula S = Pi*R^2 to find the approximate radii. Regards, Moshe. --- Julian Burgos [EMAIL PROTECTED] wrote: Hi Bart, If you only have 36 circles, the fastest way would be to use some image processing software and measure the circles by hand. One option is to use ImageJ, which you can download here http://rsb.info.nih.gov/ij/ Julian Bart Joosen wrote: Hi, Maybe this is more a programming questions than a specific R-project question, but maybe there is someone who can point me in the right direction. I have a picture of cirkels which I took with a digital camera. Now I want to use the diameter of the cirkels on the picture for analysis in R. I can use pixmap to import the picture, but how do I find the outside cirkels and calculate the diameter? I pointed out that I can use the edci package, but then I need to preprocess the data to reduce the points, otherwise it takes a long time, and my computer crashes. If you want to see such a picture, I cropped a larger one, and highlighted the cirkel which is of interest. In a real world, this is a plate with 36 cirkels, which all should
Re: [R] Excel
As far as I understand, changing the format changes the way data is displayed by Excel but this does not change the data itself - if while reading the data Excel decided that it was a date, it is being converted to an integer (the number of days since January 1, 1900 - and they mistakenly think that 1900 was a leap year) and it is stored this way. --- David Scott [EMAIL PROTECTED] wrote: On Tue, 28 Aug 2007, Robert A LaBudde wrote: If you format the column as Text, you won't have this problem. By leaving the cells as General, you leave it up to Excel to guess at the correct interpretation. Not true actually. I had converted the column to Text because I saw the interpretation as a date in the .xls file. I saved the .csv file *after* the column had been converted to Text. Looking at the .csv file in a text editor, the entry is correct. I have just rechecked this. On reopening the .csv using Excel, the entry AUG2699 had been interpreted as a date, and was showing as Aug-99. Most bizarre is that the NHI value of AUG1838 has *not* been interpreted as a date. David Scott You will note that the conversion to a date occurs immediately in Excel when you enter the value. There are many formats to enter dates. Either pre-format the column as Text, or prefix the individual entry with an ' to indicate text. A similar problem occurs in R's read.table() function when a factor has levels that can be interpreted as numbers. At 10:11 PM 8/27/2007, David wrote: A common process when data is obtained in an Excel spreadsheet is to save the spreadsheet as a .csv file then read it into R. Experienced users might have learned to be wary of dates (as I have) but possibly have not experienced what just happened to me. I thought I might just share it with r-help as a cautionary tale. I received an Excel file giving patient details. Each patient had an ID code in the form of three letters followed by four digits. (Actually a New Zealand National Health Identification.) I saved the .xls file as .csv. Then I opened up the .csv (with Excel) to look at it. In the column of ID codes I saw: Aug-99. Clicking on that entry it showed 1/08/2699. In a column of character data, Excel had interpreted AUG2699 as a date. The .csv did not actually have a date in that cell, but if I had saved the .csv file it would have. David Scott Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: [EMAIL PROTECTED] Least Cost Formulations, Ltd.URL: http://lcfltd.com/ 824 Timberlake Drive Tel: 757-467-0954 Virginia Beach, VA 23464-3239Fax: 757-467-2947 Vere scire est per causas scire __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. _ David Scott Department of Statistics, Tamaki Campus The University of Auckland, PB 92019 Auckland 1142,NEW ZEALAND Phone: +64 9 373 7599 ext 86830 Fax: +64 9 373 7000 Email:[EMAIL PROTECTED] Graduate Officer, Department of Statistics Director of Consulting, Department of Statistics __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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 multiply all dataframe rows by another dataframe's columns
Below is one way to do this: Loc1.alleles - c(1,5,6,7,8) Loc1.Freq- c(0.35, 0.15, 0.05, 0.10, 0.35) Loc1 - cbind( Loc1.alleles,Loc1.Freq) X- data.frame(Loc1) Loc2.alleles - c(1,4,6,8) Loc2.Freq - c(0.35, 0.35, 0.10, 0.20) Loc2 - cbind(Loc2.alleles, Loc2.Freq) Y - data.frame (Loc2) Z-Y[,2] %o% X[,2] Z [,1] [,2] [,3] [,4] [,5] [1,] 0.1225 0.0525 0.0175 0.035 0.1225 [2,] 0.1225 0.0525 0.0175 0.035 0.1225 [3,] 0.0350 0.0150 0.0050 0.010 0.0350 [4,] 0.0700 0.0300 0.0100 0.020 0.0700 dim(Z) - c(prod(dim(Z)),1) Z [,1] [1,] 0.1225 [2,] 0.1225 [3,] 0.0350 [4,] 0.0700 [5,] 0.0525 [6,] 0.0525 [7,] 0.0150 [8,] 0.0300 [9,] 0.0175 [10,] 0.0175 [11,] 0.0050 [12,] 0.0100 [13,] 0.0350 [14,] 0.0350 [15,] 0.0100 [16,] 0.0200 [17,] 0.1225 [18,] 0.1225 [19,] 0.0350 [20,] 0.0700 Regards, Moshe. --- Luke Neraas [EMAIL PROTECTED] wrote: Hello, I have two data frames, X and Y, with two columns each and different numbers of rows. # creation of data frame X Loc1.alleles - c(1,5,6,7,8) Loc1.Freq- c(0.35, 0.15, 0.05, 0.10, 0.35) Loc1 - cbind( Loc1.alleles,Loc1.Freq) X- data.frame(Loc1) #creation of data frame Y Loc2.alleles - c(1,4,6,8) Loc2.Freq - c(0.35, 0.35, 0.10, 0.20) Loc2 - cbind(Loc2.alleles, Loc2.Freq) Y - data.frame (Loc2) I would like a flexible way to multiply each element of second column of the X data frame by each element of the second column of the Y data frame. example of what the operation need to do: X[1,2]*Y[,2] ; X[2,2]*Y[,2].X[4,2]*Y[,2] I have worked on a variety of for loops to get this to work without success. The final result should look like a column like this all_X[,2] * all_Y[,2] 0.1225 0.1225 0.0350 0.0700 0.0525 0.0525 0.0150 0.0300 0.0175 0.0175 0.0050 0.0100 0.0350 0.0100 0.0200 0.1225 0.1225 0.0350 0.0700 any help would be greatly appreciated Luke Neraas [EMAIL PROTECTED] [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] element wise opertation between a vector and a list
Hi, Getting your e is not a problem: a- list(c(1,3),c(1,2),c(2,3)) b-c(10,20,30) aa-matrix(unlist(a),nrow=2) ee-aa+rbind(b,b) e-apply(ee,2,sum) e [1] 24 43 65 But this will not work if different list members have different number of elements. --- Yongwan Chun [EMAIL PROTECTED] wrote: I want to try to get a result of element wise addition between a vector and a list. It can be done with for statement. In order to reducing computing time, I have tried to avoid for state. If anybody give me an idea, I would apprecite it much. for example, with a b as below lines, a- list(c(1,3),c(1,2),c(2,3)) b-c(10,20,30) I would like to have a list (like d) or a vector (like e) as below. d-list(c((1+10),(3+10)),c((1+20),(2+20)),c((2+30),(3+30))) e- c((1+10)+(3+10),(1+20)+(2+20),(2+30)+(3+30)) Thanks, Yongwan Chun __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] how to compute cross correlation
The simlest answer is that X and Y are two vectors of length n then (un-normalized) cross correlation between X and Y is sum(i=1 to n) {X(i)-Xbar)*(Y(i)-Ybar)} where Xbar and Ybar are the means of X and Y respectively. You may need something different, so could you be more specific? What do you mean by ASCII - they are in an ASCII files? Regards, Moshe. --- Yogesh Tiwari [EMAIL PROTECTED] wrote: Hello R Users, How to compute cross correlation between two time series. Data is in ASCII format. I am using R on windows. Many thanks, Regards, Yogesh [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] The quadprog package
Hi Thomas, On my computer the solution is (0,0,1,0,0,0,0) (within machine accuracy) and it satisfies the constraints. --- [EMAIL PROTECTED] wrote: Hi everybody, I'm using Windows XP Prof, R 2.5.1 and a Pentium 4 Processor. Now, I want to solve a quadratic optimization program (Portfolio Selection) with the quadprog package I want to minimize (\omega'%*%\Sigma%*%\omega) Subject to (1) \iota' %*% \omega = 1 (full investment) (2) R'%*%\omega = \mu (predefined expectation value) (3) \omega \ge 0 (no short sales). Where \omega is a Nx1 vector of the weights invested in each asset \Sigma is the NxN variance-covariance matrix \iota is a 1xN vector of 1's R' is a Nx1 vector of the expected returns and \mu is a number, the postualted return of the investor I've done the following code but it doesn't make what I want it to do. The weights aren't all positive and the \mu isn't reached. What's wrong with my code? Require(quadprog) Dmat-diag(1,7,7) # muss als quadratische Matrix eingegeben werden Dmat dvec-matrix(0,7,1) # muss als Spaltenvektor eingegeben werden dvec mu-0 # (in Mio. â¬) bvec-c(1,mu,matrix(0,7,1)) # muss als Spaltenvektor eingegeben werden bvec mu_r-c(19.7,33.0,0.0,49.7, 82.5, 39.0,11.8) Amat-matrix(c(matrix(1,1,7),7*mu_r,diag(1,7,7)),9,7,byrow=T) # muss als Matrix angegeben werden, wie sie wirklich ist Amat meq-2 loesung-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2) loesung # Ãberprüfen, ob System richtig gelöst wurde loesung$solution %*% mu_r sum(loesung$solution) for (i in 1:7){ a-loesung$solution[i]=0 print(a) } Thanks in advance for your answers. __ Thomas Schwander MVV Energie Konzern-Risikocontrolling Telefon 0621 - 290-3115 Telefax 0621 - 290-3664 E-Mail: [EMAIL PROTECTED] . Internet: www.mvv.de MVV Energie AG . Luisenring 49 . 68159 Mannheim Handelsregister-Nr. HRB 1780, Amtsgericht Mannheim Vorsitzender des Aufsichtsrates: Herr Dr. Peter Kurz Vorstand: Dr. Rudolf Schulten (Vorsitzender) . Matthias Brückmann . Dr. Werner Dub . Hans-Jürgen Farrenkopf [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] interpolation
If your data is in data.txt file you can do the following: x - read.table(file=data.txt,header=TRUE) t-ISOdatetime(x[,1],x[,2],x[,3],x[,4],x[,5],x[,6]) secs - as.numeric(t) So now secs represents your time in seconds and you can use any type of interpolation you wish to interpolate co2obs. --- Yogesh Tiwari [EMAIL PROTECTED] wrote: Hello R Users, I am new to R and I have simple problem for R users. I have CO2 observations defined on time axis(yr,mo,day,hr,min,sec). (DATA ATTACHED HERE) First I want to convert time axis as one axis as 'hour' on regular interval as 1 hour. Say 00 hrs to 24hrs(jan1), 25hrs to 48hrs(jan2) and so on. Then I want to interpolate CO2 at every hour. Kindly anybody can help, Many thanks, Regards, Yogesh yr mo dy hr min sec co2obs 1998 1 9 0 35 24 365.19 1998 1 9 1 17 39 363.54 1998 1 9 1 58 41 364.24 1998 1 9 2 39 42 364.88 1998 1 9 3 20 43 365.06 1998 1 9 4 1 44 364.75 1998 1 9 4 42 45 364.77 1998 1 9 5 23 46 364.87 1998 1 9 6 4 51 364.77 1998 1 9 6 45 52 364.73 1998 1 9 7 26 54 364.76 1998 1 9 8 7 55 363.49 1998 1 23 1 6 16 364.31 1998 1 23 1 48 32 364.38 1998 1 23 2 29 33 364.67 1998 1 23 3 10 35 365.53 1998 1 23 3 51 36 365.16 1998 1 23 4 32 37 364.56 1998 1 23 5 13 40 364.62 1998 1 23 5 54 41 365.05 1998 1 23 6 35 42 365.13 1998 1 23 7 16 44 365.45 1998 1 23 7 57 49 364.77 1998 1 23 8 38 49 364.65 1998 2 3 0 44 21 362.43 1998 2 3 1 26 36 363.96 1998 2 3 2 7 38 364.59 1998 2 3 2 48 40 364.62 1998 2 3 3 29 40 366.22 1998 2 3 4 10 41 365.4 1998 2 3 4 51 41 365.34 1998 2 3 5 32 42 365.31 1998 2 3 6 13 44 364.84 1998 2 3 6 54 46 365.07 1998 2 3 7 35 51 364.84 1998 2 3 8 16 50 364.81 1998 2 19 0 35 14 363.41 1998 2 19 1 17 31 362.93 1998 2 19 1 58 32 363.86 1998 2 19 2 39 33 364.87 1998 2 19 3 20 34 366.05 1998 2 19 4 1 36 364.84 1998 2 19 4 42 36 364.77 1998 2 19 5 23 37 365.01 1998 2 19 6 4 39 365 1998 2 19 6 45 39 365.9 1998 2 19 7 26 40 366.24 1998 2 19 8 7 41 366.55 1998 3 5 0 50 20 363.13 1998 3 5 1 32 37 363.6 1998 3 5 2 13 39 364.26 1998 3 5 2 54 40 364.26 1998 3 5 3 35 41 364.39 1998 3 5 4 16 42 365.24 1998 3 5 4 57 42 365.48 1998 3 5 5 38 43 365.01 1998 3 5 6 19 44 365.43 1998 3 5 7 0 45 365.11 1998 3 5 7 41 46 368.54 1998 3 5 8 22 48 364.96 1998 3 19 0 46 36 363.25 1998 3 19 1 28 51 363.8 1998 3 19 2 9 53 364.21 1998 3 19 2 50 55 364.46 1998 3 19 3 31 58 365.61 1998 3 19 4 12 59 365.57 1998 3 19 4 53 59 365.53 1998 3 19 5 35 0 365.38 1998 3 19 6 16 3 366.23 1998 3 19 6 57 4 364.28 1998 3 19 7 38 6 367.08 1998 3 19 8 19 5 369.19 1998 4 2 1 8 1 363.76 1998 4 2 1 50 17 365.14 1998 4 2 2 31 18 365.26 1998 4 2 3 12 21 364.7 1998 4 2 3 53 24 364.04 1998 4 2 4
Re: [R] question about non-linear least squares in R
Below is one possibility: If you knew MA you would get a regular linear least-squares for parameters A,B and constant which can be easily solved. So now you can define a function f(MA) which returns that value. Now you must minimize that f - a function of one argument. It can have several local minima and so you must be careful but I believe that minimizing (even bad) function of one argument should be easier than your original problem. Regards, Moshe. P.S. if you do this I would be interested to know whether this works. --- Yu (Warren) Wang [EMAIL PROTECTED] wrote: Hi, everyone, My question is: It's not every time that you can get a converged result from the nls function. Is there any solution for me to get a reasonable result? For example: x - c(-0.06,-0.04,-0.025,-0.015,-0.005,0.005,0.015,0.025,0.04,0.06) y - c(1866760,1457870,1314960,1250560,1184850,1144920,1158850,1199910,1263850,1452520) fitOup- nls(y ~ constant + A*(x-MA)^4 + B*(x-MA)^2, start=list(constant=1000, A=1, B=-100, MA=0), control=nls.control(maxiter=100, minFactor=1/4096), trace=TRUE) For this one, I cannot get the converged result, how can I reach it? To use another funtion or to modify some settings for nls? Thank you very much! Yours, Warren __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] SQL like function?
observation %in% ID --- Takatsugu Kobayashi [EMAIL PROTECTED] wrote: Hi RUsers, I am wonder if I can search observations whose IDs matches any of the values in another vector, such as in MySQL. While I am learing MySQL for future database management, I appreciate if anyone could give me a hint. Suppose I have one 5*1 vector containing observation IDs and frequencies, and one 3*1 vector containing observation IDs. observation-c(1,2,3,4,5) ID-c(1,3,4) Then, I would like to program a code that returns a results showing matched observations like result: TRUE FALSE TRUE TRUE FALSE I am reading S programming, but I cannot find a way to do this. Thank you very much. Taka __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Problems with strsplit
unlist(strsplit(bA531F16-rep,\\-))[1] [1] bA531F16 --- Carlos Morales [EMAIL PROTECTED] wrote: Hello, I would like to know what can I do if I use strplit with a string and I want to use the middle left,I mean I have this: strsplit(bA531F16-rep,\\-) [[1]] [1] bA531F16 rep I would like to work just with bA531F16 in another variable, what could I do?, Thank you - Sé un Mejor Amante del Cine ¿Quieres saber cómo? ¡Deja que otras personas te ayuden!. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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] Generating Replicate Datasets (using loops or other means)
Hi Jonathan, What exactly do you mean by replication? Do you want to keep a1,b1,c1,... unchanged but have 30 different sets of random numbers? Regards, Moshe. --- VTLT1999 [EMAIL PROTECTED] wrote: Hello All, I have searched many help forums, message boards, etc. and I just can't apply the comments to what I need my program to do. I am running R 2.5.1 on an XP system, and my desire is to produce replicate datasets for a simulation study I am running. Essentially, I have sets of parameters (a's, b's, and c's) that define a function which produces a decimal value. This value is compared to a random uniform value, and is coded a 1 if the function is greater than the uniform value, 0 if it is = to the uniform value. My code thus far works great, but I just need it to run several times.Here we go: library(mvtnorm) library(sm) library(ltm) library(irtoys) k- 5000 set.seed(271828) t - rmvnorm(n=k,mean=c(-1,0,1),sigma=matrix(c(1,.8,.5,.8,1,.8,.5,.8,1),3,3)) #Using mv here because of the likely association of ability (theta = t) across time. t1-as.matrix(t[,1]) t2-as.matrix(t[,2]) t3-as.matrix(t[,3]) set.seed(271828) # Population item parameters (n=54) from which we will select relevant items # These are the parameters that are used in the function a - c(1.18120, 0.92613, 0.96886, 0.80503, 1.12384, 0.84073, 0.85544, 0.86801, 1.01054, 0.82278, 1.10353, 0.78865, 0.98421, 1.76071, 0.89603, 0.84671, 0.89737, 0.74775, 0.32190, 0.69730, 0.72059, 1.16762, 1.29257, 1.32902, 0.59540, 0.51022, 0.59259, 0.93951, 0.68568, 0.55649, 0.88084, 0.52940, 0.45735, 0.57560, 1.11779, 0.96984, 1.19692, 0.99102, 1.25847, 1.62555, 0.63049, 1.07807, 1.04897, 1.23138, 1.14014, 1.25230, 1.14844, 0.59287, 0.83143, 0.81723, 0.52141, 0.61980, 0.49945, 1.02749) b - c(-2.51737, -1.95897, -1.72667, -0.82988, -0.36093, 0.72554, 0.91442, 0.78061, 0.06088, 0.75733, -0.76371, 0.24552, -0.42050, 0.88232, -0.81761, 0.06466, -0.43866, -0.46042, 0.21636, -0.73147, -1.44086, -1.03718, 0.07275, -0.17197, 1.53796, -0.45631, -1.69826, -0.66506, 0.98921, 0.30714, -0.62245, 0.97253, 1.95894, 0.21277, 1.96346, 1.18825, 1.59917, -0.28401, -1.23530, -0.09671, -0.31581, -0.66149, -0.81284, -0.35399, -0.07623, 1.06442, -0.68559, 1.07591, 0.97458, 0.06436, 1.25622, 1.73954, 1.75052, 2.34088) c - c(0.0, 0.0, 0.0, 0.0, 0.19648, 0.31302, 0.26454, 0.19714, 0.06813, 0.21344, 0.0, 0.03371, 0.0, 0.16581, 0.11054, 0.08756, 0.07115, 0.26892, 0.0, 0.06883, 0.0, 0.14815, 0.32389, 0.19616, 0.17597, 0.0, 0.0, 0.04337, 0.19949, 0.20377, 0.0, 0.06243, 0.13639, 0.0, 0.18166, 0.15996, 0.20184, 0.08331, 0.24453, 0.26114, 0.16434, 0.20750, 0.32658, 0.31870, 0.45227, 0.35039, 0.31178, 0.17999, 0.22774, 0.21675, 0.10153, 0.17764, 0.15205, 0.19858) # Item parameters for generating 3PL data for all five testing occasions: # This selects the relevant parameters for a particular data generation run # Only parameters for the first testing occasion are shown to save space a1 - as.matrix(a[c(1:5,15:20,22:24,38:44)]) b1 - as.matrix(b[c(1:5,15:20,22:24,38:44)]) c1 - as.matrix(c[c(1:5,15:20,22:24,38:44)]) # Here is where I would like to begin my replications, but don't know how to make R do it. # The code below produces a matrix of 0's and 1's (which will be used by another program) # I would like to nest this in a do loop such that, say, 30 replicate datasets are produced using the #same parameters. N - nrow(t1) # number of examinees n - nrow(a1) # number of items d - 1.7 theta - t1 response - matrix (0,N,n) uni - matrix (runif(N*n),nrow = N) for (i in 1:N) { for (j in 1:n) { if ( c1[j]+(1-c1[j])/(1+exp(-d*a1[j]*(theta[i]-b1[j]))) uni[i,j] ) response[i,j] = 1 else response[i,j] = 0 } } write.table(response, file=C:/responses.dat, sep= ,row.names=FALSE, col.names=FALSE) I tried earlier nesting this in another for loop, but that indexes elements of matrices and vectors, and doesn't seem to apply to a global loop methodology. I am attempting to use replicate as we speak, but documentation is sparse (help(replicate) is nested in lapply information). Any guidance is greatly appreciated. Thanks in advance, Jonathan Beard -- View this message in context: http://www.nabble.com/Generating-Replicate-Datasets-%28using-loops-or-other-means%29-tf4418768.html#a12603580 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch mailing list
Re: [R] finding the minimum positive value of some data
Either min(diff(sort(x))[diff(sort(x))0]) or min(diff(sort(unique(x --- dxc13 [EMAIL PROTECTED] wrote: useRs, I am looking to find the minimum positive value of some data I have. Currently, I am able to find the minimum of data after I apply some other functions to it: x [1] 1 0 1 2 3 3 4 5 5 5 6 7 8 8 9 9 10 10 sort(x) [1] 0 1 1 2 3 3 4 5 5 5 6 7 8 8 9 9 10 10 diff(sort(x)) [1] 1 0 1 1 0 1 1 0 0 1 1 1 0 1 0 1 0 min(diff(sort(x))) [1] 0 The minimum is given as zero, which is clearly true, but I am interested in only the positive minimum, which is 1. Can I find this by using only 1 line of code, like I have above? Thanks! dxc13 -- View this message in context: http://www.nabble.com/finding-the-minimum-positive-value-of-some-data-tf4417250.html#a12599319 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch 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@stat.math.ethz.ch 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.