Re: [R] lmom package - Resending the email

2014-12-03 Thread Simon Zehnder
Katherine,

for a deeper understanding of differing values it makes sense to provide the 
list at least with an online description of the corresponding functions used in 
Minitab and SPSS…

Best 
Simon
On 03 Dec 2014, at 10:45, Katherine Gobin via R-help r-help@r-project.org 
wrote:

 Dear R forum
 I sincerely apologize as my earlier mail with the captioned subject, since 
 all the values got mixed up and the email is not readable. I am trying to 
 write it again. 
 My problem is I have a set of data and I am trying to fit some distributions 
 to it. As a part of this exercise, I need to find out the parameter values of 
 various distributions e.g. Normal distribution, Log normal distribution etc. 
 I am using lmom package to do the same, however the parameter values obtained 
 using lmom pacakge differ to a large extent from the parameter values 
 obtained using say MINITAB and SPSS as given below -
 _
 
 amounts =  
 c(38572.5599129508,11426.6705314315,21974.1571641187,118530.32782443,3735.43055996748,66309.5211176106,72039.2934132668,21934.8841708626,78564.9136114375,1703.65825161293,2116.89180930203,11003.495671332,19486.3296339113,1871.35861218795,6887.53851253407,148900.978055447,7078.56497101651,79348.1239806592,20157.6241066905,1259.99802108593,3934.45912233674,3297.69946631591,56221.1154121067,13322.0705174134,45110.2498756567,31910.3686613912,3196.71168501252,32843.0140437202,14615.1499458453,13013.9915051561,116104.176753387,7229.03056392023,9833.37962177814,2882.63239493673,165457.372543821,41114.066453219,47188.1677766245,25708.5883755617,82703.7378298092,8845.04197017415,844.28834047836,35410.8486123933,19446.3808445684,17662.2398792892,11882.8497070776,4277181.17817307,30239.0371267968,45165.7512343364,22102.8513746687,5988.69296597127,51345.0146170238,1275658.35495898,15260.4892854214,8861.76578480635,37647.1638704867,4979.53544046949,7012.48134772332,3385.20612391205,1911.03114395959,66886.5036605189,2223.47536156462,814.947809578378,234.028589468841,5397.4347625133,13346.3226579065,28809.3901352898,6387.69226236731,5639.42730553242,2011100.92675507,4150.63707173462,34098.7514446498,3437.10672573502,289710.315303182,8664.66947305203,13813.3867161134,208817.521491857,169317.624400274,9966.78447705792,37811.1721605562,2263.19211279927,80434.5581206454,19057.8093104899,24664.5067589624,25136.5042354789,3582.85741610706,6683.13898432794,65423.9991390846,134848.302304064,3018.55371579808,546249.641168158,172926.689143006,3074.15064180208,1521.70624812788,59012.4248281661,21226.928522236,17572.5682970983,226.646947337851,56232.2982652019,14641.0043361533,6997.94414914865)
 
 library(lmom)
 lmom  =  samlmu(amounts)
 # __
 # Normal Distribution parameters
 parameters_of_NOR  - pelnor(lmom); parameters_of_NOR
 
   mu  sigma 115148.4175945.8
   Location   Scale Minitab 115148.4 
 485173SPSS   115148.4 485173
 # __
 # Log Normal (3 Parameter) Distribution parameters
zetamu   sigma 3225.7988909.114879 
  2.240841
   LocationScale   Shape
 MINITAB   9.73361 1.76298  75.51864SPSS   
  9.73361.763  75.519   # 
 __
 
 Besides Genaralized extreme Value distributions, all the other distributions 
 e.g. Gamma, Exponential (2 parameter) distributions etc give different 
 results than MINITAB and SPSS.
 Can some one guide me?
 
 Regards
 Katherine
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Building R for better performance

2014-03-06 Thread Simon Zehnder
 calculation (vector calc) 
  1.310.38
 Creation of a 6000x6000 Hilbert matrix (matrix calc)  
0.77 0.99
 Grand common divisors of 400,000 pairs (recursion)
 0.63 0.56
 Creation of a 1000x1000 Toeplitz matrix (loops)   
   2.24 2.34
 Escoufier's method on a 90x90 matrix (mixed)  
  9.55 6.02
 Total 
 274.93
21.01
 
 Regards,
 Jonathan Anspach
 Sr. Software Engineer
 Intel Corp.
 jonathan.p.ansp...@intel.com
 713-751-9460
 
 
 -Original Message-
 From: Simon Zehnder [mailto:szehn...@uni-bonn.de] 
 Sent: Wednesday, March 05, 2014 3:55 AM
 To: Anspach, Jonathan P
 Cc: r-help@r-project.org
 Subject: Re: [R] Building R for better performance
 
 Jonathan,
 
 I myself tried something like this - comparing gcc, clang and intel on a Mac. 
 From my experiences in HPC on the university cluster (where we also use the 
 Xeon Phi, Landeshochleistungscluster University RWTH Aachen), the Intel 
 compiler has better code optimization in regard to vectorisation, etc. (clang 
 is up to now suffering from a not yet implemented OpenMP library).
 
 Here is a revolutionanalytics article about this topic: 
 http://blog.revolutionanalytics.com/2010/06/performance-benefits-of-multithreaded-r.html
 
 As I usually use the Rcpp package for C++ extensions this could give me 
 further performance. Though, I already failed when trying to compile R with 
 the Intel compiler and linking against the MKL (see my topic in the Intel 
 developer zone: http://software.intel.com/en-us/comment/1767418 and my 
 threads on the R-User list: 
 https://stat.ethz.ch/pipermail/r-sig-mac/2013-November/010472.html). 
 
 So, to your questions:
 
 1) I think that most admins do not even use the Intel compiler to compile R - 
 this seems to me rare. There are some people I know they do and I think they 
 could be aware of it - but these are only a few. As R is growing in usage and 
 I do know from regional user meetings that very large companies start using 
 it in their BI units - this should be of interest.
 
 2) I would really welcome this step because compilation with intel 
 (especially on a Mac) and linking to the MKL seems to be delicate. 
 
 I am interested in the data - so if it is possible send it via the list or 
 directly to my account. Further, could you show some code that you used for 
 the computations? 
 
 
 Best
 
 Simon
 
 
 On 04 Mar 2014, at 22:44, Anspach, Jonathan P jonathan.p.ansp...@intel.com 
 wrote:
 
 Greetings,
 
 I'm a software engineer with Intel.  Recently I've been investigating R 
 performance on Intel Xeon and Xeon Phi processors and RH Linux.  I've also 
 compared the performance of R built with the Intel compilers and Intel Math 
 Kernel Library to a default build (no config options) that uses the GNU 
 compilers.  To my dismay, I've found that the GNU build always runs on a 
 single CPU core, even during matrix operations.  The Intel build runs matrix 
 operations on multiple cores, so it is much faster on those operations.  
 Running the benchmark-2.5 on a 24 core Xeon system, the Intel build is 13x 
 faster than the GNU build (21 seconds vs 275 seconds). Unfortunately, this 
 advantage is not documented anywhere that I can see.
 
 Building with the Intel tools is very easy.  Assuming the tools are 
 installed in /opt/intel/composerxe, the process is simply (in bash shell):
 
 $ . /opt/intel/composerxe/bin/compilervars.sh intel64 $ ./configure 
 --with-blas=-L/opt/intel/composerxe/mkl/lib/intel64 -lmkl_intel_lp64 
 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm --with-lapack 
 CC=icc CFLAGS=-O2 CXX=icpc CXXFLAGS=-O2 F77=ifort FFLAGS=-O2 FC=ifort 
 FCFLAGS=-O2 $ make $ make check
 
 My questions are:
 1) Do most system admins and/or R installers know about this performance 
 difference, and use the Intel tools to build R?
 2) Can we add information on the advantage of building with the Intel tools, 
 and how to do it, to the installation instructions and FAQ?
 
 I can post my data if anyone is interested.
 
 Thanks,
 Jonathan Anspach
 Sr. Software Engineer
 Intel Corp.
 jonathan.p.ansp...@intel.com
 713-751-9460
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 R-benchmark-25_large.R

Re: [R] Building R for better performance

2014-03-05 Thread Simon Zehnder
Jonathan,

I myself tried something like this - comparing gcc, clang and intel on a Mac. 
From my experiences in HPC on the university cluster (where we also use the 
Xeon Phi, Landeshochleistungscluster University RWTH Aachen), the Intel 
compiler has better code optimization in regard to vectorisation, etc. (clang 
is up to now suffering from a not yet implemented OpenMP library).

Here is a revolutionanalytics article about this topic: 
http://blog.revolutionanalytics.com/2010/06/performance-benefits-of-multithreaded-r.html

As I usually use the Rcpp package for C++ extensions this could give me further 
performance. Though, I already failed when trying to compile R with the Intel 
compiler and linking against the MKL (see my topic in the Intel developer zone: 
http://software.intel.com/en-us/comment/1767418 and my threads on the R-User 
list: https://stat.ethz.ch/pipermail/r-sig-mac/2013-November/010472.html). 

So, to your questions:

1) I think that most admins do not even use the Intel compiler to compile R - 
this seems to me rare. There are some people I know they do and I think they 
could be aware of it - but these are only a few. As R is growing in usage and I 
do know from regional user meetings that very large companies start using it in 
their BI units - this should be of interest.

2) I would really welcome this step because compilation with intel (especially 
on a Mac) and linking to the MKL seems to be delicate. 

I am interested in the data - so if it is possible send it via the list or 
directly to my account. Further, could you show some code that you used for the 
computations? 


Best

Simon


On 04 Mar 2014, at 22:44, Anspach, Jonathan P jonathan.p.ansp...@intel.com 
wrote:

 Greetings,
 
 I'm a software engineer with Intel.  Recently I've been investigating R 
 performance on Intel Xeon and Xeon Phi processors and RH Linux.  I've also 
 compared the performance of R built with the Intel compilers and Intel Math 
 Kernel Library to a default build (no config options) that uses the GNU 
 compilers.  To my dismay, I've found that the GNU build always runs on a 
 single CPU core, even during matrix operations.  The Intel build runs matrix 
 operations on multiple cores, so it is much faster on those operations.  
 Running the benchmark-2.5 on a 24 core Xeon system, the Intel build is 13x 
 faster than the GNU build (21 seconds vs 275 seconds).  Unfortunately, this 
 advantage is not documented anywhere that I can see.
 
 Building with the Intel tools is very easy.  Assuming the tools are installed 
 in /opt/intel/composerxe, the process is simply (in bash shell):
 
 $ . /opt/intel/composerxe/bin/compilervars.sh intel64
 $ ./configure --with-blas=-L/opt/intel/composerxe/mkl/lib/intel64 
 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm 
 --with-lapack CC=icc CFLAGS=-O2 CXX=icpc CXXFLAGS=-O2 F77=ifort FFLAGS=-O2 
 FC=ifort FCFLAGS=-O2
 $ make
 $ make check
 
 My questions are:
 1) Do most system admins and/or R installers know about this performance 
 difference, and use the Intel tools to build R?
 2) Can we add information on the advantage of building with the Intel tools, 
 and how to do it, to the installation instructions and FAQ?
 
 I can post my data if anyone is interested.
 
 Thanks,
 Jonathan Anspach
 Sr. Software Engineer
 Intel Corp.
 jonathan.p.ansp...@intel.com
 713-751-9460
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Where did lost variables go

2013-12-31 Thread Simon Zehnder
A reproducible example would do well here David

Best

Simon
On 31 Dec 2013, at 02:42, David Parkhurst parkh...@indiana.edu wrote:

 I have several variables in a data frame that aren't listed by ls() 
 after I attach that data frame.  Where did they go, and how can I stop 
 the hidden ones from masking the local ones?
 Thanks for any help.
 David
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Estimation of AR(1) model by QML.

2013-12-30 Thread Simon Zehnder
Why not using optim on the likelihood in a) with normally distributed standard 
errors and for b) optim with a likelihood with t(3)-distributed standard errors?

Best

Simon

On 30 Dec 2013, at 21:19, Xuse Chuse chus...@gmail.com wrote:

 Dear Users,
 
 I am trying to estimate a model Y(t)=alpha+rho*Y(t-1)+e(t) where i know
 e(t)~t(3).
 
 a) I want to estimate (alpha, rho) by QML estimation assuming (wrongly)
 that e(t)~N(0,sigma2) and calculate the standard errors.
 b) Estimate (alpha, rho) by ML estimation assuming (correctly) e(t)~t(3)
 and compute standard errors.
 
 Can anyone help me out to figure out how I could answer this question in R?
 Thank you beforehand.
 
 Cheers,
 Chuse
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Milliken and Johnson Unbalanced Machines

2013-12-30 Thread Simon Zehnder
I don’t know anything about this topic, but looking into this book: 
http://books.google.de/books?id=3TVDQBAJpg=PA25lpg=PA25dq=milliken+and+johnson+unbalanced+machinessource=blots=lxqStiQju5sig=d5dG_cHsTzilCIklBrW9SAMKYRMhl=desa=Xei=ifrBUuvoM8qPtAbPqoCQDQved=0CGUQ6AEwBQ#v=onepageq=milliken%20and%20johnson%20unbalanced%20machinesf=false

should give you a little hint I guess.

Best

Simon

On 30 Dec 2013, at 23:37, Troels Ring tr...@gvdnet.dk wrote:

 Dear friends - reading Milliken and Johnson on messy data I failed to find R 
 code to master the unbalanced Machines data in ch 23.
 I wonder if anyone can lend me a hand - again
 Happy new year
 Troels Ring
 Aalborg, Denmark
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] testing for bimodal and for dip in between modes in R

2013-11-25 Thread Simon Zehnder
Testing for bimodality is rather testing for unimodality. Hartigan and Hartigan 
(1985) presented the Dip-Test which is implemented in the R package DipTest 
with a much better approximation of the test distribution. If the test 
statistic is too high unimodality is rejected. To estimate the dip point you 
could choose among several possibilities: (1) A very easy method is to use the 
kmeans function for a kmeans cluster and use the point in the middle of the 
connecting line between the kmeans cluster centers. (2) You could estimate a 
finite mixture distribution and take the middle of the connecting line of the 
modes.

Best

Simon
 
On 24 Nov 2013, at 20:41, Felix Breden bre...@sfu.ca wrote:

 Hi 
 I have distributions that are typically bimodal (see attached .pdf), and I 
 would like to test for bimodality, and then estimate the point between the 
 two modes, the dip in the distributions. any help would be greatly 
 appreciated.
 thanks
 felix 
 m66.junction.aln.pairwise.histogram.pdf__
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] MM estimator

2013-11-14 Thread Simon Zehnder
Thanks Martin, for making this clear to me, I thought of Pearson’s 
Method-of-Moments.
On 14 Nov 2013, at 16:08, Martin Maechler maech...@stat.math.ethz.ch wrote:

 SZ == Simon Zehnder szehn...@uni-bonn.de
on Thu, 14 Nov 2013 08:52:16 +0100 writes:
 
SZ Check the gmm package with a weighting matrix equal to
SZ the identity.  Best
 
SZ Simon
 
 No.MM is ambigous: The  gmm package is about
   Generalized Method of Moments
 
 but Izhak is really asking about MM estimation in the field of
 robust statistics, where an MM estimator is a special kind of 
 M-Estimator (and these, introduced by Peter Huber (1964), Annals,
  fortunately are well known unambigously,
  as a generalization of ML estimators (= MLE)),
 namely an M-estimator with redescending psi (== non-convex
 problem == solution typically depends on starting value),
 started with a high-breakdown point initial estimate.
 
 Izhak mentioned  nlrob() from package robustbase, and he is
 right that this does not provide an MM estimator, but only an M
 estimate started by Least Squares.
 
 Fortunately, we have had contributions to the robustbase
 package from Eduardo Conceição since last summer.
 He provided alternative versions of  nlrob(),
 notably a  nlrob.MM()  
 currently *hidden* {and undocumented}.
 
 I.e. you need
 
 install.packages(robustbase, repos=http://R-Forge.R-project.org;)
 
 library(robustbase)
 
 and then use  robustbase:::nlrob.MM()
 
 The reason for all this is that I plan to have one nlrob()
 function  but with  a new argument  'method'
 and so eventually
 
nlrob.MM(, method=MM)
 
 will correspond to current to the R-forge version
 
robustbase:::nlrob.MM()
 
 ---
 Please for all more, use the dedicated mailing list
 R-SIG-robust only (i.e. do *not* just reply-all to this
 e-mail!).
 
 Martin Maechler,
 ETH Zurich
 
 
SZ On 14 Nov 2013, at 04:37, IZHAK shabsogh
SZ ishaqb...@yahoo.com wrote:
 
 hi
 
 I have a nonlinear regression model with two parameter, i
 have estimated using nls in R and i want to also find the
 estimate using MM, someone refer me to this function
 nlrob but this is main for only M estimate. can you
 please help me findout a function in R for MM nonlinear
 regression
 
 thanks
 
 [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do
 read the posting guide
 http://www.R-project.org/posting-guide.html and provide
 commented, minimal, self-contained, reproducible code.
 
SZ __
SZ R-help@r-project.org mailing list
SZ https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do
SZ read the posting guide
SZ http://www.R-project.org/posting-guide.html and provide
SZ commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] optimization: multiple assignment problem

2013-11-14 Thread Simon Zehnder
It would be more clear if you tell, what you want to do instead of what you do 
not want to do. 

If you start with a usual cost matrix (whatever cost function you have) and you 
have to assign N to N this reduces to the well-known Munkre’s algorithm (see 
for example: http://gallery.rcpp.org/articles/minimal-assignment/). If you have 
to assign M to N, it is an extended version of the Munkre’s algorithm which has 
been covered by Bourgeois and Lassalle (1971). All this is graph theory. 

Best 

Simon
 
On 14 Nov 2013, at 19:24, Jean-Francois Chevalier 
jean-francois.cheval...@bisnode.com wrote:

 Hello,
 I'm trying to solve a multiple assignment problem.
 I found a package Adagio and its function mknapsack which
 
 maximize vstar = p(1)*(x(1,1) + ... + x(m,1)) + ... ... + p(n)*(x(1,n) + ... 
 + x(m,n))
 
 subject to w(1)*x(i,1) + ... + w(n)*x(i,n) = k(i) for i=1,...,m
 
 x(1,j) + ... + x(m,j) = 1 for j=1,...,n x(i,j) = 0 or 1 for i=1,...,m , 
 j=1,...,n ,
 It's close to what I'm trying to do except that
 1)k(i) = k for any I (not an issue)
 2)p is dependent of the item AND the knapsack
 3)each item must be assigned
 
 maximize vstar = p(1,1)*x(1,1) + ... + p(m,1)*x(m,1) + ... ... + 
 p(1,n)*x(1,n) + ... + p(m,n)*x(m,n)
 
 with p(j,i) profit of assigning item i to knapsack j
 
 subject to w(1)*x(i,1) + ... + w(n)*x(i,n) = k for i=1,...,m
 
 x(1,j) + ... + x(m,j) = 1 for j=1,...,n x(i,j) = 0 or 1 for i=1,...,m , 
 j=1,...,n ,
 It would be really helpful if you could indicate me any package, function 
 that would solve my problem?
 Thanks in advance,
 Best regards,
 
 Jean-François
 
 
 
  DISCLAIMER \ This e-mail and any attachments t...{{dropped:13}}
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] issues with calling predict.coxph.penal (survival) inside a function - subset-vector not found. Because of NextMethod?

2013-11-13 Thread Simon Zehnder
Works for me:

predicting_function(fit2,test1)
 1  2  3  4  5  6  7
-1.0481141  0.1495946  0.4492597  0.4492597  0.9982492 -0.4991246 -0.4991246

Best

Simon

On 13 Nov 2013, at 15:46, julian.bo...@elitepartner.de wrote:

 Hello everyone, 
 
 
 
 I got an issue with calling predict.coxph.penal inside a function. 
 
 
 
 Regarding the context: My original problem is that I wrote a function that
 uses predict.coxph and survfit(model) to predict
 
 a lot of survival-curves using only the basis-curves for the strata (as
 delivered by survfit(model) ) and then adapts them with 
 
 the predicted risk-scores. Because there are cases where my new data has
 strata which didn't exist in the original model I exclude 
 
 them, using a Boolean vector inside the function.
 
 I end up with a call like this: predict (coxph_model,
 newdata[subscript_vector,] ) 
 
 
 
 This works fine for coxph.model, but when I fit a model with a spline
 (class coxph.penal), I get an error: 
 
 Error in `[.data.frame`(newdata, [subscript_vector, ) : object
 '[subscript_vector ' not found
 
 
 
 I suppose this is because of NextMethod, but I am not sure how to work
 around it. I also read a little bit about all those
 matching-and-frame-issues, 
 
 But must confess I am not really into it. 
 
 
 
 I attach a reproducible example. 
 
 Any help or suggestions of work-arounds will be appreciated. 
 
 
 
 Thanks 
 
 Julian
 
 
 
 version
 
   _   
 
 platform   x86_64-w64-mingw32  
 
 arch   x86_64  
 
 os mingw32 
 
 system x86_64, mingw32 
 
 status 
 
 major  3   
 
 minor  0.1 
 
 year   2013
 
 month  05  
 
 day16  
 
 svn rev62743   
 
 language   R   
 
 version.string R version 3.0.1 (2013-05-16)
 
 nickname   Good Sport
 
 
 
 
 
 ##TEST-DATA
 
 
 
 # Create the simplest test data set 
 
 test1 - data.frame(time=c(4,3,1,1,2,2,3), 
 
  status=c(1,1,1,0,1,1,0), 
 
  x=c(0,2,1,1,1,0,0), 
 
  sex=c(0,0,0,0,1,1,1)) 
 
 
 
 # Fit a stratified model 
 
 fit1 - coxph(Surv(time, status) ~ x + strata(sex), test1) 
 
 summary(fit1)
 
 
 
 #fit stratified wih spline
 
 fit2 - coxph(Surv(time, status) ~ pspline(x, df=2) + strata(sex), test1) 
 
 summary(fit2)
 
 
 
 #function to predict within
 
 
 
 predicting_function - function(model, newdata){
 
  subs -vector(mode='logical', length=nrow(newdata))
 
  subs[1:length(subs)]- TRUE
 
 
 
  ret - predict (model, newdata=newdata[subs,])
 
  return(ret)
 
 }
 
 
 
 predicting_function(fit1, test1) # works
 
 
 
 predicting_function(fit2,test1) #doesnt work - Error in
 `[.data.frame`(newdata, subs, ) : object 'subs' not found
 
# probably because of NextMethod
 
 
 
 #
 
 traceback()
 
 #12: `[.data.frame`(newdata, subs, )
 
 #11: newdata[subs, ]
 
 #10: is.data.frame(data)
 
 #9: model.frame.default(data = newdata[subs, ], formula = ~pspline(x, 
 
 #   df = 2) + strata(sex), na.action = function (object, ...) 
 
 #   object)
 
 #8: model.frame(data = newdata[subs, ], formula = ~pspline(x, df = 2) + 
 
 #   strata(sex), na.action = function (object, ...) 
 
 #   object)
 
 #7: eval(expr, envir, enclos)
 
 #6: eval(tcall, parent.frame())
 
 #5: predict.coxph(model, newdata = newdata[subs, ])
 
 #4: NextMethod(predict, object, ...)
 
 #3: predict.coxph.penal(model, newdata = newdata[subs, ])
 
 #2: predict(model, newdata = newdata[subs, ]) at #5
 
 #1: predicting_function(fit2, test1)
 
 
 
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] MM estmator

2013-11-13 Thread Simon Zehnder
Check the gmm package with a weighting matrix equal to the identity.

Best

Simon
 
On 14 Nov 2013, at 04:37, IZHAK shabsogh ishaqb...@yahoo.com wrote:

 hi
 
 I have a nonlinear regression model with two parameter, i have estimated 
 using nls in R
 and i want to also find the estimate using MM, someone refer me to this 
 function nlrob
 but this is main for only M estimate. can you please help me findout  a 
 function in R for MM nonlinear regression
 
 thanks
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] S4; Setter function is not chaning slot value as expected

2013-11-09 Thread Simon Zehnder
If you want to set a slot you have to refer to it:

ac@CustomerID - “54321” 

or you use your setter:

ac - CustomerID(ac, “54321”)

What you did was creating a new symbol CustomerID referring to the String 
“54321”
CustomerID - “54321”
CustomerID
[1] “54321”

Best

Simon

On 09 Nov 2013, at 15:31, daniel schnaider dschnai...@gmail.com wrote:

 It is my first time programming with S4 and I can't get the setter fuction
 to actually change the value of the slot created by the constructor.
 
 I guess it has to do with local copy, global copy, etc. of the variable -
 but, I could't find anything relevant in documentation.
 
 Tried to copy examples from the internet, but they had the same problem.
 
 # The code
setClass (Account ,
   representation (
   customer_id = character,
   transactions = matrix)
)
 
 
Account - function(id, t) {
new(Account, customer_id = id, transactions = t)
}
 
 
setGeneric (CustomerID-, function(obj,
 id){standardGeneric(CustomerID-)})
setReplaceMethod(CustomerID, Account, function(obj, id){
obj@customer_id - id
obj
})
 
ac - Account(12345, matrix(c(1,2,3,4,5,6), ncol=2))
ac
CustomerID - 54321
ac
 
 #Output
 ac
An object of class Account
Slot customer_id:
[1] 12345
 
Slot transactions:
 [,1] [,2]
[1,]14
[2,]25
[3,]36
 
 # CustomerID is value has changed to 54321, but as you can see it does't
 CustomerID - 54321
 ac
An object of class Account
Slot customer_id:
[1] 12345
 
Slot transactions:
 [,1] [,2]
[1,]14
[2,]25
[3,]36
 
 
 Help!
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] S4 vs S3. New Package

2013-11-09 Thread Simon Zehnder
This depends very often of on the developer and what he is comfortable with. I 
like S4 classes, as I come from C++ and S4 classes approximate C++ classes at 
least more than S3 classes do (Reference Classes would do so even more and I 
know very good R programmers liking these most).

1) I wrote a package for MCMC simulation with S4 classes carrying all simulated 
values - fast enough for me: in less than 1.5 secs I have my sample of 100.000 
values together with several other 100T values like log-likelihoods, posterior 
hyper parameters, etc. I watch out for not copying too often an object but 
sometimes it is not avoidable. 

2) That is not true: 

Books:
http://www.amazon.de/Software-Data-Analysis-Programming-Statistics/dp/0387759352/ref=sr_1_1?ie=UTF8qid=1384014486sr=8-1keywords=John+chambers+data
http://www.amazon.de/Programming-Data-Language-John-Chambers/dp/0387985034/ref=sr_1_4?ie=UTF8qid=1384014486sr=8-4keywords=John+chambers+data

Online:
https://www.rmetrics.org/files/Meielisalp2009/Presentations/Chalabi1.pdf
https://www.stat.auckland.ac.nz/S-Workshop/Gentleman/S4Objects.pdf

And for a bunch of packages look into the Bioconductor packages.

Best

Simon

On 09 Nov 2013, at 16:22, daniel schnaider dschnai...@gmail.com wrote:

 Hi,
 
 I am working on a new credit portfolio optimization package. My question is
 if it is more recommended to develop in S4 object oriented or S3.
 
 It would be more naturally to develop in object oriented paradigm, but
 there is many concerns regarding S4.
 
 1) Performance of S4 could be an issue as a setter function, actually
 changes the whole object behind the scenes.
 
 2) Documentation. It has been really hard to find examples in S4. Most
 books and articles consider straightforward S3 examples.
 
 Thanks,
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] select .txt from .txt in a directory

2013-11-08 Thread Simon Zehnder
I do not understand the question. If you already know the names what is the 
problem to select the files by names? 

If you have the names but not inside of R you have to find a name pattern to 
avoid typing them in. Is there a pattern, e.g. da.txt, db.txt, dc.txt? 


On 08 Nov 2013, at 18:33, Zilefac Elvis zilefacel...@yahoo.com wrote:

 Hi,
 I have 300 .txt files in a directory. Out of this 300, I need just 100 of the 
 files.
 I have the names of the 100 .txt files which are also found in the 300 .txt 
 files.
 How can I extract only the 100 .txt files from the 300 ,txt files?
 
 e.g given d1.txt, ds.txt, dx.txt, df.txt...d300.txt, how can I select only 
 d1.txt and df.txt? Remember, I have 300 of such and want to extract 100 of 
 them with names known.
 
 Thanks for your great help.
 Atem.
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] select .txt from .txt in a directory

2013-11-08 Thread Simon Zehnder
If you want to type in the names by hand, you can simply use read.table to load 
them into R … I still don’t get the aim of your text file handling


On 08 Nov 2013, at 18:51, Zilefac Elvis zilefacel...@yahoo.com wrote:

 All files are text files. They are found in a folder on my computer. 
 Assume that I know the names of some of the files I want to select from the 
 300 txt files.
 How can I do this in R.
 Atem.
 
 
 On Friday, November 8, 2013 11:44 AM, Simon Zehnder szehn...@uni-bonn.de 
 wrote:
 I do not understand the question. If you already know the names what is the 
 problem to select the files by names? 
 
 If you have the names but not inside of R you have to find a name pattern to 
 avoid typing them in. Is there a pattern, e.g. da.txt, db.txt, dc.txt? 
 
 
 On 08 Nov 2013, at 18:33, Zilefac Elvis zilefacel...@yahoo.com wrote:
 
  Hi,
  I have 300 .txt files in a directory. Out of this 300, I need just 100 of 
  the files.
  I have the names of the 100 .txt files which are also found in the 300 .txt 
  files.
  How can I extract only the 100 .txt files from the 300 ,txt files?
  
  e.g given d1.txt, ds.txt, dx.txt, df.txt...d300.txt, how can I select only 
  d1.txt and df.txt? Remember, I have 300 of such and want to extract 100 of 
  them with names known.
  
  Thanks for your great help.
  Atem.
 
  [[alternative HTML version deleted]]
  
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] select .txt from .txt in a directory

2013-11-08 Thread Simon Zehnder
Elvis,

first, keep things on the list - so others can learn and comment. Second, as 
Sarah already commented: We do not like to open unsolicited binary attachments 
on the list. Sarah gives a good hint how to post data to the list.

What I would do so far is use the matching columns to get the names you need 
from BTemperature: 

temp_inv - read.table(Temperature Inventory, … ) (here I would change the 
.xlsx to a .csv and use read.csv instead of read.table)
btemp - read.table(“BTemperature_Stations.txt”, … ) (again think about 
converting via Excel to .csv - it makes things far more easy) 

Check ?read.table for options - you gonna need them.

Then match
mynames - btemp[(temp_inv[,3] %in% btemp[, 3]), 2]

Now you have the names of the stations and if your .txt files are named by the 
stations you can do something like:

for (name in mynames) {
tmp.table - read.table(paste(“path/to/your/Homog_daily_min_temp/“, name, 
“.txt”, sep = “”), … )
…. do things with the data
}



Best

Simon
 
On 08 Nov 2013, at 19:26, Zilefac Elvis zilefacel...@yahoo.com wrote:

 Hi Simon,
 Attached are my data files.
 Btemperature_Stations is my main file.
 Temperature inventory is my 'wanted' file and is a subset of 
 Btemperature_Stations.
 Using column 3 in both files, select the files in Temperature inventory from 
 Btemperature_Stations.
 The .zip file contains the .txt files which you will extract to a folder and 
 do the selection in R.
 
 Thanks,
 Atem.
  
 
 
 
 
 On Friday, November 8, 2013 11:54 AM, Simon Zehnder szehn...@uni-bonn.de 
 wrote:
 If you want to type in the names by hand, you can simply use read.table to 
 load them into R … I still don’t get the aim of your text file handling
 
 
 On 08 Nov 2013, at 18:51, Zilefac Elvis zilefacel...@yahoo.com wrote:
 
  All files are text files. They are found in a folder on my computer. 
  Assume that I know the names of some of the files I want to select from the 
  300 txt files.
  How can I do this in R.
  Atem.
  
  
  On Friday, November 8, 2013 11:44 AM, Simon Zehnder szehn...@uni-bonn.de 
  wrote:
  I do not understand the question. If you already know the names what is the 
  problem to select the files by names? 
  
  If you have the names but not inside of R you have to find a name pattern 
  to avoid typing them in. Is there a pattern, e.g. da.txt, db.txt, dc.txt? 
  
  
  On 08 Nov 2013, at 18:33, Zilefac Elvis zilefacel...@yahoo.com wrote:
  
   Hi,
   I have 300 .txt files in a directory. Out of this 300, I need just 100 of 
   the files.
   I have the names of the 100 .txt files which are also found in the 300 
   .txt files.
   How can I extract only the 100 .txt files from the 300 ,txt files?
   
   e.g given d1.txt, ds.txt, dx.txt, df.txt...d300.txt, how can I select 
   only d1.txt and df.txt? Remember, I have 300 of such and want to extract 
   100 of them with names known.
   
   Thanks for your great help.
   Atem.
  
  [[alternative HTML version deleted]]
   
   __
   R-help@r-project.org mailing list
   https://stat.ethz.ch/mailman/listinfo/r-help
   PLEASE do read the posting guide 
   http://www.R-project.org/posting-guide.html
   and provide commented, minimal, self-contained, reproducible code.
  
  
  
 
 
 BTemperature_Stations.txtTempearture 
 inventory.xlsxHomog_daily_min_temp.zip

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R with openblas and atlas

2013-11-01 Thread Simon Zehnder
There is no R code following


On 01 Nov 2013, at 05:34, Li Bowen bowenl...@gmail.com wrote:

 Hi,
 
 I have been trying to build R with optimized BLAS library.
 
 I am using a Ubuntu 13.10 x86_64 desktop, on which I am able to build R
 with openblas without any problem:
 
 #BEGIN_SRC sh
 ./configure --enable-BLAS-shlib --enable-R-shlib LIBnn=lib --disable-nls 
 --with-blas=-L/usr/lib/openblas-base/ -lopenblas --enable-memory-profiling
 make
 sudo make install
 #END_SRC
 
 However, on redhat 5.9, I am not able to install openblas.
 Firstly, there is no pre-built package, even for later version of
 redhat.
 Secondly, I am not able to build openblas from source, actually not even
 able to install newer gcc from source.
 
 I then tried to install ATLAS on redhat and build R with it. -t 2: 2 threads
 #BEGIN_SRC sh
 ../configure --shared -t 2 -b 64 -D c -DPentiumCPS=1600
 --with-netlib-lapack-tarfile=/path-to-lapack-3.4.2.tgz
 make build
 make check
 make ptcheck
 make time
 make install
 
 # R
 ../configure --enable-BLAS-shlib --enable-R-shlib --disable-nls
 --enable-memory-profiling --with-blas=-L/usr/local/atlas/lib
 -lptf77blas -lpthread -latlas
 make
 sudo make install
 #END_SRC
 
 Installation is successful. However when I run the following code, only
 one thread is used.
 
 I have looked through lots of manuals and forums and couldn't find an
 answer. Please advise. Thanks a lot.
 
 -- 
 Sincerely,
 Bowen
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Combinations of values in two columns

2013-11-01 Thread Simon Zehnder
You could use the data.table package
require(data.table)

DT - data.table(Friend1 = sample(LETTERS, 10, replace = TRUE), Friend2 = 
sample(LETTERS, 10, replace = TRUE), Indicator = 1)
ALL - data.table(unique(expand.grid(DT)))
setkey(ALL) 
OTHERS - ALL[!DT]
OTHERS[, Indicator := 0]

RESULT - rbind(DT, ALL)

Best

Simon



On 01 Nov 2013, at 10:32, Thomas thomas.ches...@nottingham.ac.uk wrote:

 I have data that looks like this:
 
 Friend1, Friend2
 A, B
 A, C
 B, A
 C, D
 
 And I'd like to generate some more rows and another column. In the new column 
 I'd like to add a 1 beside all the existing rows. That bit's easy enough.
 
 Then I'd like to add rows for all the possible directed combinations of rows 
 not included in the existing data. So for the above I think that would be:
 
 A, D
 D, A
 B, C
 C, B
 B, D
 C, A
 D, B
 D, C
 
 and then put a 0 in the column beside these.
 
 Can anyone suggest how to do this?
 
 I'm using R version 2.15.3.
 
 Thank you,
 
 Thomas Chesney
 This message and any attachment are intended solely for the addressee and may 
 contain confidential information. If you have received this message in error, 
 please send it back to me, and immediately delete it.   Please do not use, 
 copy or disclose the information contained in this message or in any 
 attachment.  Any views or opinions expressed by the author of this email do 
 not necessarily reflect the views of the University of Nottingham.
 
 This message has been checked for viruses but the contents of an attachment
 may still contain software viruses which could damage your computer system, 
 you are advised to perform your own checks. Email communications with the 
 University of Nottingham may be monitored as permitted by UK legislation.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] set difference between two data frames

2013-10-31 Thread Simon Zehnder
You could e.g. take the data.table package (every data.table is a data.frame) 
and make a join:

dt.x - data.table(x)
dt.y - data.table(y)
merge.xy - x[y, nomatch = 0]
diff.xy - x[!merge.xy]



On 31 Oct 2013, at 21:41, Yasin Gocgun entropy...@gmail.com wrote:

 Thanks. Actually, I forgot to add that both have the same number of columns.
 
 On Thu, Oct 31, 2013 at 4:07 PM, Bert Gunter gunter.ber...@gene.com wrote:
 lapply() setdiff() by columns.
 
 Unless you have failed to tell us something, you almost certainly will
 not get a data frame (same number of rows/column) as your answer.
 
 -- Bert
 
 On Thu, Oct 31, 2013 at 12:58 PM, Yasin Gocgun entropy...@gmail.com wrote:
 Hi,
 
 I have two data frames, say, x and y, where y is a subset of x. How
 can I find the set difference of these two data frames (i.e., x-y)?
 
 Thanks,
 
 --
 Yasin Gocgun
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 
 
 --
 
 Bert Gunter
 Genentech Nonclinical Biostatistics
 
 (650) 467-7374
 
 
 
 -- 
 Yasin Gocgun
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] speed of makeCluster (package parallel)

2013-10-28 Thread Simon Zehnder
See library(help = parallel”)


On 28 Oct 2013, at 17:19, Arnaud Mosnier a.mosn...@gmail.com wrote:

 Hi all,
 
 I am quite new in the world of parallelization and I wonder if there is a
 way to increase the speed of creation of a parallel socket cluster. The
 time spend to include threads increase exponentially with the number of
 thread considered and I use of computer with two 8 cores CPU and thus
 showing a total of 32 threads in windows 7.
 Currently, I use the default parameters (type = PSOCK), but is there any
 fine tuning parameters that I can use to take advantage of this system ?
 
 Thanks in advance for your help !
 
 Arnaud
 
 R version 3.0.1 (2013-05-16)
 Platform: x86_64-w64-mingw32/x64 (64-bit)
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] speed of makeCluster (package parallel)

2013-10-28 Thread Simon Zehnder
First,

use only the number of cores as a number of thread - i.e. I would not use hyper 
threading, etc.. Each core has its own caches and it is always fortunate if a 
process has enough memory; hyper threads use all the same cache on the core 
they are running on. detectCores() gives me for example 4 - I know I have 2. I 
would therefore call makeCluster() with nnode = 2. mcaffinity() lets you 
perform a technique called process-pinning (see process affinity) and is only 
possible if the OS supports it. It makes sometimes sense to assign certain 
processes to certain CPUs such that each process has enough memory in caches 
(e.g. for a 16 Core machine using 8 processes on CPUs 1, 3, 5, 7, 9, 11, 13 and 
15; so each process has the cache of two CPUs). 
A lot of functions though are not available for Windows. 

At first it comes always the problem you want to solve and then you look how 
much memory will be used in a process and how much you have (more often the 
memory bandwidth is the bottleneck and not the computing power). Look at the 
architecture of your chips (how much L1 Cache, how much L2 cache). Then you 
decide how many cores to use and if it makes sense to pin processes to certain 
cores. 

There are no general recipes for parallel computing - each problem is 
different. Some problems are even not scalable. 

Simon


On 28 Oct 2013, at 17:51, Arnaud Mosnier a.mosn...@gmail.com wrote:

 Thanks Simon,
 
 I already read the parallel vignette but I did not found what I wanted.
 May be you can be more specific on a part of the document that can provide me 
 hints !
 
 Arnaud
 
 
 2013/10/28 Simon Zehnder szehn...@uni-bonn.de
 See library(help = parallel”)
 
 
 On 28 Oct 2013, at 17:19, Arnaud Mosnier a.mosn...@gmail.com wrote:
 
  Hi all,
 
  I am quite new in the world of parallelization and I wonder if there is a
  way to increase the speed of creation of a parallel socket cluster. The
  time spend to include threads increase exponentially with the number of
  thread considered and I use of computer with two 8 cores CPU and thus
  showing a total of 32 threads in windows 7.
  Currently, I use the default parameters (type = PSOCK), but is there any
  fine tuning parameters that I can use to take advantage of this system ?
 
  Thanks in advance for your help !
 
  Arnaud
 
  R version 3.0.1 (2013-05-16)
  Platform: x86_64-w64-mingw32/x64 (64-bit)
 
[[alternative HTML version deleted]]
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] DEoptim inconsistent output

2013-09-29 Thread Simon Zehnder
It's always a good thing to trace the optimization. There can be a lot of 
reasons:

Is the function correctly implemented?

Is it defined for all values? If not, is it defined for all values in the 
constrained parameter space? Is it defined for the constraint? If it is not 
defined for the constrained (e.g. for zero, take the lower bound to be 1e-10; 
approximately zero)

How do the numerical gradients behave? They rely on difference methods, could 
be undefined for some values. 

So, as long as we have no error messages or more precise information, we cannot 
help you.


Best

Simon

On Sep 29, 2013, at 7:17 PM, Aya Anas aa...@feps.edu.eg wrote:

 Dear all,
 
 
 
 I wrote an R optimization code using the DEoptim package that is used to
 minimize a numerical integration subject to a single constraint. I noticed
 that the code takes a long time to give an output. The resulting output
 doesn't make sense at all. I got parameters that don't satisfy the
 constraint. In addition, when i substituted with the resulting parameters
 in the objective function, I didn't obtain the resulting optimal objective
 function value. What might be the reason behind this?
 
 .
 
 Thanks,
 
 Aya
 
 
 -- 
 Best Regards
 
 Faculty of Economics  Political Science
 Cairo University
 Tel:(202)35728055-(202)35728116-(202)35736608-(202)35736605
 Fax:(202)35711020
 Follow us on twitter:https://twitter.com/fepsnews
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Memory distribution using foreach

2013-09-27 Thread Simon Zehnder
First of all LSF is a batch scheduling software. It usually expects an .lsf 
script. Usually the compilers on a cluster are interchangeable via the 'module 
switch unload module load module' and MPI-2 is the message passing 
interface standard. This is also rather an topic for the high-performance R 
list. 

Next, doMC is a multicore package registering cores on one machine - AFAIK, 
i.e. you have to work on one machine with the 24 cores (inform yourself on the 
hardware on your cluster - there should also be introduction presentations 
online! To know what hardware you use and what architecture it has is the first 
step! Try 'bhosts' on your shell to see what hosts are available). If you want 
to use several machines, your backend for foreach should be doMPI and not doMC 
(see http://cran.r-project.org/web/packages/doMPI/vignettes/doMPI.pdf). 

If you found your host, you have to write an lsf-script like the one following 
(for OpenMP on ONE machine - using 24 cores, in most cases this suffices. 
Further, it is faster as you do not have to wait that long because you have to 
use just ONE machine. If you have BULL clusters - take these. A lot of cores 
32/64… and a lot of memory)

So in your case, write a script with extension .lsf containing:

### using the zsh shell
#!/usr/bin/env zsh 

### Job name
#BSUB -J OpenMP

### File/path where output will be written, the %J is the job ID
#BSUB -o OpenMP.%J

### (OFF) Different file for STDERR, if not to be merged with STDOUT
# #BSUB -e OpenM.e%J

### Request the time you need for execution in minutes 
### The format for the parameter is: [hour:]minute,
### that means for 80 minutes you could also use this: 1:20
#BSUB -W 3:00

### Request virtual memory you need for your job in MB (per process)
#BSUB -M 1024

### Request higher amount of stack site (per process)
#BSUB -S 1024
 
### Request the number of compute slots you want to use
#BSUB -n 24

### Specify your mail address 
#BSUB pko...@bgc-jena.mpg.de

### Send a mail when job is done
#BSUB -N

### Use esub for OpenMP 
#BSUB -a openmp

### (OFF) As R is usually compiled via gcc I would load the gcc module on your 
cluster 
# module switch pgi gcc/4.6

### (OFF) load another OpenMP (check which one is usually loaded!! should be 
now OpenMP 4.0) version than the default one
# module switch openmp openmp/3.0

### Set stack and address limits
ulimit -s unlimited
ulimit -v unlimited 

### Change to the work directory 
cd /home/your_username/

### Execute your application (make sure, that R can be loaded via 'R' on the 
shell!!!) 
R --no-restore --no-save --quiet --slave  your_R_script.R



In your R script file, load the packages

library(doMC)
library(foreach)

 
registerDoMC(24)  ## now, foreach knows the backend.

forach(...) %dopar% …..

## save your stuff to your work- or home directory (csv or database)

quit()

---

Then you send the script to LSF via 

bsub - my_LSF_script.lsf

Look via 'bjobs' if it is is send and what's its status (PEND or RUN). If the 
status is RUN you can look via 'bpeek your_job_ID' what the output looks like, 
while it runs.

Best

Simon



On Sep 27, 2013, at 10:48 PM, pakoun pko...@bgc-jena.mpg.de wrote:

 Dear R users,
 I am struggling with memory issues and try to understand a few things. I am
 using an LSF cluster with PGI compiler and parallel mpi2 computing (whatever
 does that means..) and i submit a job like:
 
 bsub -R rusage[mem=3] -q queue -n 24 R CMD BATCH arguments.. 
 myjob.r ..log
 
 According to that I am asking for 24 cores and 30GB RAM. 
 
 Then I am using 
 library(doMC)
 registerDoMC(24) 
 
 and a foreach loop either simple or nested with the %dopar% command. 
 
 1.  this 30 GB will be distributed among the 24 jobs or each will take 30?
 2. If i dont ask the -n 24 argument still the foreach loop will run in
 parallel as i check with TOP command. What is the purpose of using it? Just
 to reserve the nodes from other users?
 
 Thank you
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Memory-distribution-using-foreach-tp4677133.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Installing R, R packages

2013-09-24 Thread Simon Zehnder
Had yesterday something pretty similar on the 73 installing benchmark:

Error in untar2(tarfile, files, list, exdir, restore_times) :
  incomplete block on file

The downloaded source packages are in

‘/private/var/folders/n9/zxfxcd01557dc06bf3c1njy0gn/T/RtmpkwQyuQ/downloaded_packages’
Warning messages:
1: In download.file(url, destfile, method, mode = wb, ...) :
  downloaded length 247786 != reported length 474551
2: In download.file(url, destfile, method, mode = wb, ...) :
  downloaded length 280942 != reported length 403808
3: In install.packages(benchmark) :
  installation of package ‘relations’ had non-zero exit status
4: In install.packages(benchmark) :
  installation of package ‘benchmark’ had non-zero exit status

Today doing the same worked perfectly.


On Sep 24, 2013, at 6:10 PM, David Winsemius dwinsem...@comcast.net wrote:

 
 On Sep 24, 2013, at 12:47 AM, Jeff Newmiller wrote:
 
 I have stopped using the Berkeley mirror, and just automatically use UCLA 
 due to missing packages. However, I feel no compulsion to extrapolate and 
 say that there is some sort of corruption going on at CRAN mirrors because 
 it is only one data point.
 
 I have not been having this sort of difficulty with the Berkeley mirror. 
 Never have seen this particular error. Occasionally see the mirror down on a 
 weekend, but that happens with all mirrors. (I do use Berkeley as my standard 
 mirror and I do download a relatively large number of packages especially at 
 the point where I call `update.packages`.)
 
 
 -- 
 David.
 
 Sent from my phone. Please excuse my brevity.
 
 David Arnold dwarnol...@suddenlink.net wrote:
 All,
 
 Consider this attempt:
 install.packages(car)
 trying URL
 'http://cran.cnr.Berkeley.edu/bin/macosx/contrib/3.0/car_2.0-19.tgz'
 Content type 'application/x-gzip' length 1326903 bytes (1.3 Mb)
 opened URL
 ==
 downloaded 1.1 Mb
 
 car/data/Rdata.rdb: Truncated tar archive
 tar: Error exit delayed from previous errors.
 
 The downloaded binary packages are in
 
 /var/folders/qE/qEavkZWTFMmxjncuY+HnqE+++TI/-Tmp-//Rtmpnxi1hV/downloaded_packages
 Warning messages:
 1: In download.file(url, destfile, method, mode = wb, ...) :
 downloaded length 1126960 != reported length 1326903
 2: 'tar' returned non-zero exit code 1
 
 
 This seems to be happening frequently at the cran berkeley site. I've
 also
 had students try to install R from the Berkeley site and it just
 doesn't
 work for them.
 
 I had another student tell me today he tried all sorts of mirrors and
 could
 not get R and knitr installed until he tried the Washington site.
 
 Is there some sort of corruption going on at the CRAN mirrors?
 
 D.
 
 
 
 
 David Winsemius, MD
 Alameda, CA, USA
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Checking a large data set for normality

2013-09-24 Thread Simon Zehnder
Check the Jarque-Bera Test for univariate testing 
(http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/tseries/html/jarque.bera.test.html)
 and Mardia's test for multivariate testing 
(http://www.inside-r.org/packages/cran/MVN/docs/mardia.test).


On Sep 24, 2013, at 6:44 PM, steric rd023...@ohio.edu wrote:

 Hello,
 
 I have a large data set that includes many soil parameters (i.e. pH, calcium
 levels, enzyme activity, etc) Does anyone have any input as to the easiest
 way to check a large data set for normality? Is there an R function/package
 that can do this all at once?
 
 Thank you in advance,
 
 
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Checking-a-large-data-set-for-normality-tp4676858.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Hope u have some time for 2 more questions

2013-09-20 Thread Simon Zehnder
Hi,

as far as I know, there is no limitation on data size in regard to foreach. You 
should reserve though enough memory for your application on the cluster (via 
ulimit -s unlimited and ulimit -v unlimited). 

Furthermore I would check the following: 
Check if there are two versions of R on the cluster/your home directory on the 
frontend (LSF loads this frontend environment and uses the R version installed 
there). If you have two R executables (R and R64) make sure you use the 64bit 
version.

Run R and call memory.limit() to see what are the limits of memory in your 
system. 

If this is limited to sizes below your needed sizes, increase it by calling R 
in the LSF script with the options --max-mem-size=YourSize and if you get 
errors of kind  cannot allocate vector of size you should also use 
--max-vsize=YourVSize. 

Then, check if there is a memory leak in your application: If you compiled R 
with the --enable-memory-profiling you can use Rprof to do this otherwise you 
must rely on profiling instruments given by the cluster environment (I think 
you work there as well with modules, so type in the shell 'module avail' for 
listing available modules). 

If you detect a memory leak or if you see, that at certain points in your 
algorithms some objects are not used anymore call rm(ObjectName) and gc() for 
garbage collection. 


To your nested loop using foreach: That is a highly delicate issue in parallel 
computing and for the foreach syntax I refer to the must-read 
http://cran.r-project.org/web/packages/foreach/vignettes/nested.pdf. 

Using nested loops should be considered carefully in regard to organizing the 
nesting. In C++ you have the ability to determine how many cores should work on 
which loop. In the foreach environment using doMC this seems to me not 
possible. 


And, please keep the discussion to the r-help mailing list, so others can learn 
from it and researchers with more experience can also leave comments. 


Best

Simon


On Sep 19, 2013, at 9:24 PM, pko...@bgc-jena.mpg.de wrote:

 Hi again,
 
 if you have some time I would like to bother you again with 2 more questions. 
 After your response the parallel code is working perfect but when I implement 
 that to the real case (big matrices) I get an error for not numeric dimension 
 and i guess that again it returns NULL or something. Are you aware if foreach 
 loop can handle only a certain size objects? the equation that I am using 
 includes 3 objects with 2Gb size each. 
 
 The second question has to deal with the cores that foreach uses. Although I 
 am asking to our cluster (LSF) to give me certain number of cpus, and also i 
 am specifing that with
 library(doMC)
 registerDoMC(n) 
 
 it seems from the top command that I am using all the cores. I am using 2 
 foreach as nest  foreach(i in 1:16){
 foreach(j in 1:10)  etc etc..
 maybe i should do something with this kind of nest? I am not aware about that.
 
 I am sorry for the long text , and thank you for your nice solution
 
 _
 Sent from http://r.789695.n4.nabble.com
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Renaming variables

2013-09-20 Thread Simon Zehnder
You haven't said yet, what object your 'data file' is. If you mean a data.frame 
I would use colnames(dataName) - c(Col1Name, col2Name, ….)

Best
Simon

On Sep 20, 2013, at 4:10 PM, Preetam Pal lordpree...@gmail.com wrote:

 Hi,
 
 I guess this is pretty basic.
 
 I have 25 variables in the data file (name: score), i.e. X1,X2,.,X25.
 
 I dont want to use score$X1, score$X2 everytime I use these variables.
 
 Is there a way I can rename all these variables as simply X1,X2,.X25
 without writing 25 lines of code, one line for renaming each variable (eg:
 X1=score.X1  X2=score.X2  and so on) ?
 
 Thanks for your help.
 
 Regards,
 Preetam
 
 -- 
 Preetam Pal
 (+91)-9432212774
 M-Stat 2nd Year, Room No. N-114
 Statistics Division,   C.V.Raman
 Hall
 Indian Statistical Institute, B.H.O.S.
 Kolkata.
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] callNextMethod with dots argument

2013-09-19 Thread Simon Zehnder
Kien,

if you want to add variables in a function definition that is predefined by a 
Generic and calls CallNextMethod you have to add the '…' argument as well.

 setMethod(do,signature(a=Achild),
  function(a,msg,...)  {
print(do Achild)
callNextMethod()
  })

Best

Simon

On Sep 19, 2013, at 8:58 AM, Kiên Kiêu kien.k...@jouy.inra.fr wrote:

 Hi,
 
 I met a problem when invoking callNextMethod within a method associated with 
 a generic function taking ... as an argument.
 
 Here is the code
 
 setClass(Aparent,representation(x=numeric,y=numeric))
 setClass(Achild,contains=Aparent)
 
 setGeneric(do,def=function(a,...) standardGeneric(do))
 setMethod(do,signature(a=Aparent),
  function(a,msg) {
print(do Aparent)
  })
 setMethod(do,signature(a=Achild),
  function(a,msg)  {
print(do Achild)
callNextMethod()
  })
 
 myA - new(Achild)
 buf - do(a=myA)   # works
 buf - do(a=myA,msg=bonjour) # error
 
 The last call yields the following error message:
 
 Error in callNextMethod() :
  in processing 'callNextMethod', found a '...' in the matched call, but no 
 corresponding '...' argument
 
 which I do not understand. Replacing ... by msg in setGeneric makes it 
 work. But I don't like this limitation so much (unless I understand it).
 
 Regards.
 
 Kien
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R-3.0.1 g77 errors

2013-09-18 Thread Simon Zehnder
On my systems Linux Scientific and Mac OS X I use as well for the F77 the 
gfortran compiler and this works. You could give it a trial.

Best

Simon

On Sep 18, 2013, at 3:14 PM, Prigot, Jonathan jpri...@partners.org wrote:

 I am trying to build R-3.0.1 on our SPARC Solaris 10 system, but it
 fails part way through with g77 errors. Has anyone run into this? Any
 suggestions? For what it's worth, R-2.15.1 is the last one to build
 error free for us.
 ===
 Jon Prigot
 
 R is now configured for sparc-sun-solaris2.10
 
 Source directory:  .
 Installation directory:/usr/local
 
 C compiler:gcc -std=gnu99  -g -O2
 Fortran 77 compiler:   g77  -g -O2
 
 C++ compiler:  g++  -g -O2
 Fortran 90/95 compiler:gfortran 
 Obj-C compiler: 
 
 Interfaces supported:  X11, tcltk
 External libraries:readline, ICU
 Additional capabilities:   PNG, JPEG, TIFF, NLS
 Options enabled:   shared BLAS, R profiling
 
 Recommended packages:  yes
 
 make 
 ...
 
 g77 -fPIC  -g -O2 -ffloat-store -c dlamch.f -o dlamch.o
 dlamch.f: In function `dlamch':
 dlamch.f:89: warning:
   INTRINSIC  DIGITS, EPSILON, HUGE, MAXEXPONENT,
  ^
 Reference to unimplemented intrinsic `DIGITS' at (^) (assumed EXTERNAL)
 dlamch.f:89: 
   INTRINSIC  DIGITS, EPSILON, HUGE, MAXEXPONENT,
  ^
 Invalid declaration of or reference to symbol `digits' at (^) [initially
 seen at (^)]
 dlamch.f:89: warning:
   INTRINSIC  DIGITS, EPSILON, HUGE, MAXEXPONENT,
  ^
 Reference to unimplemented intrinsic `EPSILON' at (^) (assumed EXTERNAL)
 
 -- 
 Jonathan M. Prigot jpri...@partners.org
 Partners Healthcare Systems
 
 
 
 The information in this e-mail is intended only for th...{{dropped:18}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] foreach returns null first object in the list

2013-09-17 Thread Simon Zehnder
Use an extra call to return:

post.ls - foreach(i =1:2, .verbose = TRUE) %dopar% {

fun - func.list[[i]]
 if (i == 1) return(fun(Xa, Sa))
 if (i == 2) return(fun(Ta, Sa))
}

Best

Simon

On Sep 17, 2013, at 11:55 AM, pakoun pko...@bgc-jena.mpg.de wrote:

 Dear all,
 
 I am sending you a copy of a demo that i am trying to make that work and
 embedded in a more complicated structure. The problem is that i want to get
 a list of 2 objects (matrices in specific, and thats why i am not specifing
 a .combine argument if i am not wrong..), but the first element of the list
 is just null. the second one is perfectly fine. Any idea?
 Thank you in advance
 
 library(doMC)
 registerDoMC(2)
 
 Xa-matrix(1:100,ncol=10)
 Sa-matrix(101:200,ncol=10)
 Ta-matrix(1:100,ncol=10)
   get.Xp - function(Xa,Sa){
 result - Xa%*%Sa
 return(result)
   }
   get.Sa.post - function(Xa,Ta){
 
 result - Ta%*%Sa
 return(result)
   }
 
  func.list - list(get.Xp,get.Sa.post)   
 
  post.ls - foreach(i =1:2) %dopar% {
 
 fun - func.list[[i]]
 if(i==1){fun(Xa,Sa)}
 if(i==2){fun(Ta,Sa)}
 
 }
 
 post.ls
 [[1]]
 NULL
 
 [[2]]
[,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]
 [1,] 154855 169455 184055 198655 213255 227855 242455 257055 271655 286255
 [2,] 155910 170610 185310 200010 214710 229410 244110 258810 273510 288210
 [3,] 156965 171765 186565 201365 216165 230965 245765 260565 275365 290165
 [4,] 158020 172920 187820 202720 217620 232520 247420 262320 277220 292120
 [5,] 159075 174075 189075 204075 219075 234075 249075 264075 279075 294075
 [6,] 160130 175230 190330 205430 220530 235630 250730 265830 280930 296030
 [7,] 161185 176385 191585 206785 221985 237185 252385 267585 282785 297985
 [8,] 162240 177540 192840 208140 223440 238740 254040 269340 284640 299940
 [9,] 163295 178695 194095 209495 224895 240295 255695 271095 286495 301895
 [10,] 164350 179850 195350 210850 226350 241850 257350 272850 288350 303850
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/foreach-returns-null-first-object-in-the-list-tp4676303.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] foreach returns null first object in the list

2013-09-17 Thread Simon Zehnder
Use an extra call to 'return()' as posted below. 

On Sep 17, 2013, at 5:55 PM, pakoun pko...@bgc-jena.mpg.de wrote:

 This is the output with the debugging message. I don't really understand what
 is the problem.
 
  post.ls - foreach(i =1:2, .verbose=T) %dopar% {
 +
 +  fun - func.list[[i]]
 +  if(i==1){fun(Xa,Sa)}
 +  if(i==2){fun(Ta,Sa)}
 +
 +  }
 numValues: 2, numResults: 0, stopped: TRUE
 got results for task 1
 numValues: 2, numResults: 1, stopped: TRUE
 returning status FALSE
 got results for task 2
 numValues: 2, numResults: 2, stopped: TRUE
 calling combine function
 evaluating call object to combine results:
  fun(accum, result.1, result.2)
 returning status TRUE
 
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/foreach-returns-null-first-object-in-the-list-tp4676303p4676329.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Unrecognized token

2013-09-17 Thread Simon Zehnder
Maybe you should escape the single quotes something like  '  id  ' ?

On Sep 17, 2013, at 6:03 PM, srecko joksimovic sreckojoksimo...@gmail.com 
wrote:

 Hi,
 
 when I generate query using sqldf library, like this:
 query = paste(paste(select * from tbl_user where student_id = , id,
sep=),  order by date_time, sep=)
 
 student - sqldf(query)
 
 everything works fine in case the id is 21328, 82882, or something like
 that. But, when id is something like 78789D, there is an error:
 Error in sqliteExecStatement(con, statement, bind.data) :
  RS-DBI driver: (error in statement: unrecognized token: 78789D)
 
 I tried replacing single quotes with double, but it still doesn't work...
 
 thanks,
 Srecko
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] set breakpoint in debug

2013-09-16 Thread Simon Zehnder
You could just use debug(f) and then when the Browser opens and the loop begins 
type 'c', that jumps over the loop to next line after the loop. 

Best

Simon

On Sep 16, 2013, at 9:05 PM, Hui Du hui...@dataventures.com wrote:

 Hi All,
 
 I need some help regarding how to set up a breakpoint in debug. For example, 
 I have a very simple/naïve function (a useless function just for demo)
 
 f = function()
 {
x = 10;
len = 100;
 
a = 1;
for(i in 1:len)
{
a = a * i;
}
 
 y = x + a;
 y;
 }
 
 
 If I need to debug it, I can run debug(f). After I go into the debugger, if I 
 want to skip the loop and stop in the statement y = y + a, directly, how to 
 do that? I know R has a function named setBreakpoint but I have never used it 
 correctly.
 
 Your help is highly appreciated.
 
 HXD
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Package installation and path.package

2013-09-09 Thread Simon Zehnder
I am following your suggestion and move this discussion to the R-devel list. 

Best 

Simon
On Sep 9, 2013, at 7:58 AM, Prof Brian Ripley rip...@stats.ox.ac.uk wrote:

 On 09/09/2013 02:09, David Winsemius wrote:
 
 On Sep 8, 2013, at 8:00 AM, Simon Zehnder wrote:
 
 Dear R-Users and R-Devels,
 
 I am writing right now my own package that makes use of 'tempfile' and 
 there within with 'path.package'. When I install it, I get the error: Error 
 in path.package(mypackage) : none of the packages are loaded.
 
 I understand the error, but I would like to have a workaround. How can I 
 give the path to the package I am actually installing without getting this 
 error?
 
 (We do not have the code so this is speculation.) Your packages should be 
 assumed to be available in one of the directories in .libPaths()
 
 Not until the package is fully installed.  We do have no idea what is going 
 on here ... and R-devel seems the appropriate list.
 
 -- 
 Brian D. Ripley,  rip...@stats.ox.ac.uk
 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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Package installation and path.package

2013-09-08 Thread Simon Zehnder
Dear R-Users and R-Devels,

I am writing right now my own package that makes use of 'tempfile' and there 
within with 'path.package'. When I install it, I get the error: Error in 
path.package(mypackage) : none of the packages are loaded. 

I understand the error, but I would like to have a workaround. How can I give 
the path to the package I am actually installing without getting this error? 


Best

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] [dfoptim] 'Error in fn(ginv(par), ...) : object 'alpha' not found'

2013-09-03 Thread Simon Zehnder
Hi Carlos,

your problem is a wrong definition of your Likelihood function. You call 
symbols in the code (alpha, beta) which have no value assigned to. When L the 
long calculation in the last lines is assigned to L alpha and beta do not 
exist. The code below corrects it. But you have a problem with a divergent 
integral when calling integrate. A problem you can surely fix as you know what 
your function is doing.

Likelihood_cov - function(params, x, tx, T, IS) {
 r - params[1]
 alpha_zero - params[2]
 s - params[3]
 beta_zero - params[4]
 gamma_1 - params[5]
 gamma_2 - params[6]
 data$alpha - alpha_zero*exp(-gamma_1*IS)
 data$beta - beta_zero*exp(-gamma_2*IS)
 f - function(x, tx, T, alpha, beta)
 {
   g - function(y)
 (y + alpha)^(-( r + x))*(y + beta)^(-(s + 1))
   integrate(g, tx, T)$value
 }
 integral - mdply(data, f)
 L -
exp(lgamma(r+x)-lgamma(r)+r*(log(alpha_zero)-log(alpha_zero+T))-x*log(alpha_zero+T)+s*(log(beta_zero)-log(beta_zero+T)))+exp(lgamma(r+x)-lgamma(r)+r*log(alpha_zero)+log(s)+s*log(beta_zero)+log(integral$V1))
 f - -sum(log(L))
 return (f)
}


Best

Simon


On Sep 3, 2013, at 1:28 PM, Carlos Nasher carlos.nas...@googlemail.com wrote:

 Dear R helpers,
 
 I have problems to properly define a Likelihood function. Thanks to your
 help my basic model is running quite well, but I have problems to get the
 enhanced version (now incorporating covariates) running.
 
 Within my likelihood function I define a variable 'alpha'. When I want to
 optimize the function I get the error message:
 
 'Error in fn(ginv(par), ...) : object 'alpha' not found'
 
 I think it's actually not a problem with the optimization function (nmkb),
 but with the Likelihood function itself. I do not understand why 'alpha' is
 a missing object. 'alpha' should be part of the dataframe 'data' (as 'beta'
 should be too), like 'x', 'tx', ''T. But it obviously isn't.
 
 Here's a minimum example which reproduces my problem:
 
 ##
 
 library(plyr)
 library(dfoptim)
 
 ### Sample data ###
 x - c(3, 0, 2, 5, 1, 0, 0, 1, 0, 2)
 tx - c(24.57, 0.00, 26.86, 34.57, 2.14, 0.00, 0.00, 8.57, 0.00, 14.29)
 T - c(33.29, 30.71, 31.29, 34.57, 36.00, 35.43, 31.14, 33.86, 35.71, 35.86)
 IS - c(54.97, 13.97, 122.33, 110.84, 30.72, 14.96, 30.72, 20.74, 29.16,
 83.00)
 data - data.frame(x=x, tx=tx, T=T)
 rm(x, tx, T)
 
 ### Likelihood function ###
 Likelihood_cov - function(params, x, tx, T, IS) {
  r - params[1]
  alpha_zero - params[2]
  s - params[3]
  beta_zero - params[4]
  gamma_1 - params[5]
  gamma_2 - params[6]
  data$alpha - alpha_zero*exp(-gamma_1*IS)
  data$beta - beta_zero*exp(-gamma_2*IS)
  f - function(x, tx, T, alpha, beta)
  {
g - function(y)
  (y + alpha)^(-( r + x))*(y + beta)^(-(s + 1))
integrate(g, tx, T)$value
  }
  integral - mdply(data, f)
  L -
 exp(lgamma(r+x)-lgamma(r)+r*(log(alpha)-log(alpha+T))-x*log(alpha+T)+s*(log(beta)-log(beta+T)))+exp(lgamma(r+x)-lgamma(r)+r*log(alpha)+log(s)+s*log(beta)+log(integral$V1))
  f - -sum(log(L))
  return (f)
 }
 
 ### ML optimization ###
 params - c(0.2, 5, 0.2, 5, -0.02, -0.02)
 fit - nmkb(par=params, fn=Likelihood_cov, lower=c(0.0001, 0.0001, 0.0001,
 0.0001, -Inf, -Inf), upper=c(Inf, Inf, Inf, Inf, Inf, Inf), x=data$x,
 tx=data$tx, T=data$T, IS=IS)
 
 ##
 
 
 Maybe you could give me a hint were the flaw in my code is. Many thanks in
 advance.
 Carlos
 
 
 -
 Carlos Nasher
 Buchenstr. 12
 22299 Hamburg
 
 tel:+49 (0)40 67952962
 mobil:+49 (0)175 9386725
 mail:  carlos.nas...@gmail.com
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How to catch errors regarding the hessian in 'optim'

2013-09-02 Thread Simon Zehnder
Dear John,

thank you very much for your answer. I take a look at these packages (Rvmmin 
and Rcgmin). That sounds very interesting. 

For the example: The method relies on data which I always try to avoid to send 
on the r-help list - not that my data is confidential - but it becomes even 
more cumbersome when sending datasets over the list. I made the experience, 
that in such case answers are rare. Maybe you have a suggestion if and how 
users should send data with their code. I am interested in your opinion and 
thankful for sharing it.

Best

Simon


On Sep 2, 2013, at 4:42 PM, Prof J C Nash (U30A) nas...@uottawa.ca wrote:

 This may be one of the many mysteries of the internals of L-BFGS-B, which I 
 have found fails from time to time. That is one of the reasons for Rvmmin and 
 Rcgmin (and hopefully sooner rather than later Rtn - a truncated Newton 
 method, currently working for unconstrained problems, but still glitchy for 
 bounds constraints). These are all-R codes so that users and developers can 
 get inside to control special situations.
 
 If you have a test problem -- the infamous reproducible example -- there are 
 several of us who can likely help to sort out your troubles.
 
 JN
 
 
 On 13-09-02 06:00 AM, r-help-requ...@r-project.org wrote:
 Message: 10
 Date: Sun, 1 Sep 2013 17:09:35 +0200
 From: Simon Zehnderszehn...@uni-bonn.de
 To: R-help helpr-help@r-project.org
 Subject: [R] How to catch errors regarding the hessian in 'optim'
 Message-ID:eb37670e-8544-4c89-9172-245eb6cc5...@uni-bonn.de
 Content-Type: text/plain; charset=us-ascii
 
 Dear R-Users and R-Developers,
 
 in a comparison between two different estimation approaches I would like to 
 catch errors from optim regarding the hessian matrix.
 
 I use optim with method = L-BFGS-B thereby relying on numerical 
 differentiation for the hessian matrix. I do know, that the estimation 
 approach that uses numerical optimization has sometimes problems with 
 singular hessian matrices and I consider it as one of its disadvantages of 
 this method. To show the frequency of such problems in my simulation study I 
 have to set 'hessian = TRUE' and to collect the errors from optim regarding 
 the hessian.
 
 Now I am a little stucked how I could catch specifically errors from the 
 hessian matrix in 'optim'. I do know that such errors are thrown most 
 certainly from function 'La_solve' in Lapack.c. Does anyone has an idea how 
 I could solve this task (clearly with tryCatch but how to choose only errors 
 for the hessian)?
 
 
 Best
 
 Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] How to catch errors regarding the hessian in 'optim'

2013-09-01 Thread Simon Zehnder
Dear R-Users and R-Developers,

in a comparison between two different estimation approaches I would like to 
catch errors from optim regarding the hessian matrix.

I use optim with method = L-BFGS-B thereby relying on numerical 
differentiation for the hessian matrix. I do know, that the estimation approach 
that uses numerical optimization has sometimes problems with singular hessian 
matrices and I consider it as one of its disadvantages of this method. To show 
the frequency of such problems in my simulation study I have to set 'hessian = 
TRUE' and to collect the errors from optim regarding the hessian.

Now I am a little stucked how I could catch specifically errors from the 
hessian matrix in 'optim'. I do know that such errors are thrown most certainly 
from function 'La_solve' in Lapack.c. Does anyone has an idea how I could solve 
this task (clearly with tryCatch but how to choose only errors for the hessian)?


Best 

Simon
 
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Understanding S4 method dispatch

2013-08-14 Thread Simon Zehnder
Ambiguity is indeed detected by R and the user is informed on it. But in the 
case of Hadley's example, I still believe, that the specific multiple 
inheritance structure creates this behavior. If you call:

showMethods(f)
Function: f (package .GlobalEnv)
x=A, y=A
x=AB, y=AB
(inherited from: x=B, y=B)
x=B, y=B

you see, that AB inherited the method from B. As B inherits already from A and 
AB does as well, AB would have two A objects, the one from its inheritance from 
A and the one from B. In the case of direct inheritance, the only slot is the 
.xData slot, which is NULL. If one class contains another one, it contains all 
its slots. As the .xData slot is the same for each class A, B and AB, R must 
remove the ambiguity of slots, this is done by allowing only one slot .xData, 
which is the one from B. So in some way A is hidden by B in AB - but this is 
actually a common technique in OOP. In C++ for example we have the same problem 
and we can define which members should be inherited if multiple inheritance is 
aware, by using the keyword 'virtual'.

Towards Hervé's second point: 

If we build a class C:

setClass(C, contains = c(A))
setMethod(f, signature(C, C), function(x, y) C-C)

And then construct a multiple inheritance structure within a class BC:

setClass(BC, contains = c(B, C))
bc - new(BC)

we see, that indeed the first lexicographical signature is chosen: 

showMethods(f)
Function: f (package .GlobalEnv)
x=A, y=A
x=B, y=B
x=BC, y=BC
(inherited from: x=B, y=B)
x=C, y=C

f(bc, bc)
[1] B-B

In this case, we have a 'true' tie in the inheritance structure: B from A and C 
from A, any one of the two As have to be taken into the BC object and the one 
from B is chosen. But method dispatch is in this case independent of the second 
level inheritance. B and C are not related directly to each other, so method 
dispatch is chosen lexicographically.

In my opinion the reason for the behavior lies in the specific multiple 
inheritance structure between AB, B and A. 

Best

Simon




On Aug 14, 2013, at 2:11 AM, Hervé Pagès hpa...@fhcrc.org wrote:

 Hi Hadley,
 
 I suspect that the dispatch algorithm doesn't realize that selection
 is ambiguous in your example. For 2 reasons:
 
  (1) When it does realize it, it notifies the user:
 
setClass(A, NULL)
setGeneric(f, function(x, y) standardGeneric(f))
setMethod(f, signature(A, ANY), function(x, y) A-ANY)
setMethod(f, signature(ANY, A), function(x, y) ANY-A)
a - new(A)
 
  Then:
 
 f(a, a)
Note: method with signature ‘A#ANY’ chosen for function ‘f’,
 target signature ‘A#A’.
 ANY#A would also be valid
[1] A-ANY
 
  (2) When dispatch is ambiguous, the first method lexicographically in
  the ordering should be selected (according to ?Methods). So it
  should be A#A, not B#B.
 
 So it looks like a bug to me...
 
 Cheers,
 H.
 
 
 On 08/13/2013 06:08 AM, Hadley Wickham wrote:
 Hi all,
 
 Any insight into the code below would be appreciated - I don't
 understand why two methods which I think should have equal distance
 from the call don't.
 
 Thanks!
 
 Hadley
 
 # Create simple class hierarchy
 setClass(A, NULL)
 setClass(B, A)
 
 a - new(A)
 b - new(B)
 
 setGeneric(f, function(x, y) standardGeneric(f))
 setMethod(f, signature(A, A), function(x, y) A-A)
 setMethod(f, signature(B, B), function(x, y) B-B)
 
 # These work as I expect
 f(a, a)
 f(b, b)
 
 setClass(AB, contains = c(A, B))
 ab - new(AB)
 
 # Why does this return B-B? Shouldn't both methods be an equal distance?
 f(ab, ab)
 
 # These both return distance 1, as I expected
 extends(AB, A, fullInfo=TRUE)@distance
 extends(AB, B, fullInfo=TRUE)@distance
 # So why is signature(B, B) closer than signature(A, A)
 
 
 -- 
 Hervé Pagès
 
 Program in Computational Biology
 Division of Public Health Sciences
 Fred Hutchinson Cancer Research Center
 1100 Fairview Ave. N, M1-B514
 P.O. Box 19024
 Seattle, WA 98109-1024
 
 E-mail: hpa...@fhcrc.org
 Phone:  (206) 667-5791
 Fax:(206) 667-1319
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Grabbing numbers inside a character string

2013-08-14 Thread Simon Zehnder
Thomas,

se ?sub and regular expression. That should make it. Further, see the package 
gsubfn

Best

Simon


On Aug 14, 2013, at 2:14 PM, Thomas thomas.ches...@nottingham.ac.uk wrote:

 I have a string that contains something like:
 
 ...verified email at neu.edubrCited by 
 99853br/td/tr/table/div/div/div/divdiv...
 
 and I'd like to extract the number next to the text Cited by  - so it will 
 be whatever numbers are beside Cited by  until a non-numeric character is 
 reached.
 
 Can anyone suggest a good way of doing this please?
 
 Thank you,
 
 Thomas
 This message and any attachment are intended solely for the addressee and may 
 contain confidential information. If you have received this message in error, 
 please send it back to me, and immediately delete it.   Please do not use, 
 copy or disclose the information contained in this message or in any 
 attachment.  Any views or opinions expressed by the author of this email do 
 not necessarily reflect the views of the University of Nottingham.
 
 This message has been checked for viruses but the contents of an attachment
 may still contain software viruses which could damage your computer system, 
 you are advised to perform your own checks. Email communications with the 
 University of Nottingham may be monitored as permitted by UK legislation.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Understanding S4 method dispatch

2013-08-14 Thread Simon Zehnder
Because the signature is always (A,A) or (B,B). Then, as in AB we have A and B 
and no relationship between A and B, R chooses the method lexicographically. 
The result is as expected: f for A is chosen. 

If you would do something like: 

setClass(A, contains = list)
setClass(B, contains = list)
setClass(AB, contains = c(A, B))

setGeneric(f, function(x, y) standardGeneric(f))
setMethod(f, signature(A, ANY), function(x, y) A-ANY)
setMethod(f, signature(ANY, A), function(x, y) ANY-A)
setMethod(f, signature(B, B), function(x, y) B-B)

ab - new(AB)
f(ab, ab)

You get an ambiguity, as there is no function with signature pair that can be 
directly matched for A. But for B such a signature pair exists. So R chooses B. 
If you would specify again a function for signature (A,A) we are back:

setMethod(f, signature(A, A), function(x, y) A-A)
f(ab, ab)


Best 

Simon


On Aug 14, 2013, at 5:25 PM, Hadley Wickham h.wick...@gmail.com wrote:

 In my opinion the reason for the behavior lies in the specific multiple 
 inheritance structure between AB, B and A.
 
 So what if we don't make such a weird inheritance structure, and
 instead have A and B inherit from a common parent:
 
 setClass(A, contains = list)
 setClass(B, contains = list)
 setClass(AB, contains = c(A, B))
 
 setGeneric(f, function(x, y) standardGeneric(f))
 setMethod(f, signature(A, A), function(x, y) A-A)
 setMethod(f, signature(B, B), function(x, y) B-B)
 
 ab - new(AB)
 f(ab, ab)
 
 Why isn't there a warning about ambiguous dispatch?
 
 Hadley
 
 -- 
 Chief Scientist, RStudio
 http://had.co.nz/

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Understanding S4 method dispatch

2013-08-13 Thread Simon Zehnder
Hadley,

The class AB inherits from A and from B, but B already inherits from class A. 
So actually you only have an object of class B in your object of class AB. When 
you call the function f R looks for a method f for AB objects. It does not find 
such a method and looks for a method of the object inherited from, B. Such a 
method is present and is then executed. 

The inheritance structure has to be changed. The behavior is actually desired, 
as if this behavior weren't given a diamond class inheritance would be fatal. 


Best

Simon



On Aug 13, 2013, at 3:08 PM, Hadley Wickham h.wick...@gmail.com wrote:

 Hi all,
 
 Any insight into the code below would be appreciated - I don't
 understand why two methods which I think should have equal distance
 from the call don't.
 
 Thanks!
 
 Hadley
 
 # Create simple class hierarchy
 setClass(A, NULL)
 setClass(B, A)
 
 a - new(A)
 b - new(B)
 
 setGeneric(f, function(x, y) standardGeneric(f))
 setMethod(f, signature(A, A), function(x, y) A-A)
 setMethod(f, signature(B, B), function(x, y) B-B)
 
 # These work as I expect
 f(a, a)
 f(b, b)
 
 setClass(AB, contains = c(A, B))
 ab - new(AB)
 
 # Why does this return B-B? Shouldn't both methods be an equal distance?
 f(ab, ab)
 
 # These both return distance 1, as I expected
 extends(AB, A, fullInfo=TRUE)@distance
 extends(AB, B, fullInfo=TRUE)@distance
 # So why is signature(B, B) closer than signature(A, A)
 
 -- 
 Chief Scientist, RStudio
 http://had.co.nz/
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Understanding S4 method dispatch

2013-08-13 Thread Simon Zehnder
If you take an example which works with slots,

setClass(A, representation(a = numeric)
setClass(B, contains = c(A), representation(b = numeric))
a - new(A, a = 2)
b - new(B, a = 3, b = 2)

setClass(AB, contains = c(A, B))
new(AB, a = 2, b = 3)

You see, that there is only one @a slot, the one inherited from B, that B 
inherits from A. If this were not the case, which slot should be taken, if we 
would call @a? To avoid this kind of ambiguity, only one A class is inherited 
to AB: the one B already inherits from A. 

You could create a class, that contains another A object in a slot:

setClass(AandB, contains = c(B), representation(A = A))
new(AandB, a = 2, b = 3, A = new(A, a = 3))

Now back to your example: as there is only one A object inside the B object 
which is contained by the AB object, method dispatch works the way as it 
should: It looks for a method f for an AB object. It does not find one. Then it 
looks for a method f for the contained B object (as this is the only one 
contained in AB) and it finds a method. Then it calls this method on the B part 
of the object AB and the result is B-B

Best

Simon



On Aug 13, 2013, at 4:24 PM, Hadley Wickham h.wick...@gmail.com wrote:

 The class AB inherits from A and from B, but B already inherits from class 
 A. So actually you only have an object of class B in your object of class 
 AB. When you call the function f R looks for a method f for AB objects. It 
 does not find such a method and looks for a method of the object inherited 
 from, B. Such a method is present and is then executed.
 
 The inheritance structure has to be changed. The behavior is actually 
 desired, as if this behavior weren't given a diamond class inheritance would 
 be fatal.
 
 Are you sure? That behaviour doesn't agree with the description of
 method dispatch given in ?Methods, not with getClass(AB) which shows
 that AB inherits from both A and B. (I totally agree that this is a
 bad idea, and unlikely to be useful in real life, but I'm trying to
 understand the details of S4 dispatch)
 
 getClass(AB)
 Class AB [in .GlobalEnv]
 
 Slots:
 
 Name:  .xData
 Class:   NULL
 
 Extends:
 Class B, directly
 Class A, directly
 Class .NULL, by class A, distance 2
 Class NULL, by class A, distance 3, with explicit coerce
 Class OptionalFunction, by class A, distance 4, with explicit coerce
 Class optionalMethod, by class A, distance 4, with explicit coerce
 
 -- 
 Chief Scientist, RStudio
 http://had.co.nz/

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Outliers and overdispersion

2013-08-13 Thread Simon Zehnder
I do not know what you are exactly estimating, but if it is about count models 
and the model fit gets better when you drop the outliers, it does not say, that 
the model is now more correct. It just says, if the data were without the 
outliers, this model would fit good. 

Overdispersion in count data is sometimes a cue, that you have a mixture 
distribution as the generating process - for example instead of one, K 
different (sub)species of birds which were aggregated in the count data. In 
this case a mixture (negative binomial)- distribution with K components could 
fit the data better. 


Best

Simon
 
On Aug 13, 2013, at 5:28 PM, Marta Lomas lomasv...@hotmail.com wrote:

 
 
 
 Hi  again,  
 
 I have a question on some outliers that I have in my response variable (wich 
 are bird counts). At the beginning I did not drop them
 out because they are part of the normal counts and I considered them 
 ecologically correct. 
 
 However, I 
 tried some of the same models without ouliers and the AICs are thus better. I
 also have nice significances this way...
 
 
 So would you say that, even though the outliers are right 
 observations and taking into consideration that already the negative binomial 
 distribution that I am using is accounting for the some of the overdispersion 
 due to the outliers, it is
 better to drop them out as the models fit better this way? 
 
 
 Thanks for your patience!
 
 :)
 
 
 
 
 
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] decimal separator from comma to dot

2013-08-09 Thread Simon Zehnder
Hi,

I think this could help you: 
https://stat.ethz.ch/pipermail/r-help/2008-January/152378.html

Best

Simon

On Aug 9, 2013, at 12:19 PM, maxbre mbres...@arpa.veneto.it wrote:

 This is my reproducible example
 
 df-structure(list(IDANT = c(37837L, 37838L, 37839L, 37840L, 37841L, 
 37842L, 37843L, 40720L, 40721L, 40722L), N_TX = c(6L, 6L, 6L, 
 4L, 1L, 1L, 1L, 2L, 2L, 1L), TILT = c(0L, 0L, 0L, 0L, 6L, 6L, 
 6L, 0L, 0L, 0L), DIREZIONE = c(50L, 220L, 110L, 50L, 220L, 110L, 
 50L, 170L, 70L, 270L), DATA_INI = structure(c(2L, 2L, 2L, 2L, 
 2L, 2L, 2L, 1L, 1L, 1L), .Label = c(20/10/2004, 29/08/2002
 ), class = factor), POT_TX = structure(c(4L, 4L, 4L, 3L, 2L, 
 2L, 2L, 1L, 1L, 1L), .Label = c(10, 11,5, 4, 8), class = factor)),
 .Names = c(IDANT, 
 N_TX, TILT, DIREZIONE, DATA_INI, POT_TX), row.names = c(NA, 
 10L), class = data.frame)
 
 The data frame “df” it’s a “simplified snapshot” from an oracle table
 imported via odbc (by using the package RODBC) ;
 now, the problem with the data frame “df” is that one variable “POT_TX” have
 the decimal separator formatted with comma instead of dot; 
 so I’m in the need to reformat the variable from factor to numeric type to
 perform some useful calculations
 
 this is my code I worked so far
 
 #function to change comma to dot
 myfun - function(x) {sub(,,.,x)}
 
 #apply the function to all variables
 new - apply(df, 2, myfun )
 newdf - data.frame(apply(new, 2, as.numeric))
 str(newdf)
 
 #apply the function to one variables
 var-as.numeric(sapply(df[,6], FUN=myfun))
 df$POT_TX-var
 str(df)
 
 
 my questions:
 
 1.is it possible to use the more convenient function “apply” with 
 reference
 to just some variables and not all of them? I’ve been reading on the help
 that “Where X has named dimnames, it can be a character vector selecting
 dimension names.”, what does exactly means that?
 
 2.do you know some more convenient methods to perform comma to dot
 substitution in just some selected variables of a data frame? To note that
 the desired final result is again a dataframe
 
 3.do you know any workaround to set the decimal delimiter for a variable
 (field) “at the early steps” when connecting via odbc and then querying to
 the table?
 For these “early steps” I’ve been using something like:
 library(RODBC)
 con-odbcConnect(dsn, uid=user, pwd=password)
 df-sqlQuery(con, select * from table.name;)
 
 thank you all for the eventual support
 
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/decimal-separator-from-comma-to-dot-tp4673414.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Method dispatch in S4

2013-08-09 Thread Simon Zehnder
Hi Martin,

thank you very much for this profound answer! Your added design advice is very 
helpful, too! 

For the 'simple example': Sometimes I am still a little overwhelmed from a 
certain setting in the code and my ideas how I want to handle a process. But I 
learn from session to session. In future I will also span the lines more than 
80 columns. I am used to the indent in my vim editor. 

I have one further issue: I do know, that you are one of the leading developers 
of the bioconductor package which uses (as far as I have read) extensively OOP 
in R. Is there a package you could suggest to me to learn from by reading and 
understanding the code? Where can I find the source code? 

Best

Simon


On Aug 8, 2013, at 10:00 PM, Martin Morgan mtmor...@fhcrc.org wrote:

 On 08/04/2013 02:13 AM, Simon Zehnder wrote:
 So, I found a solution: First in the initialize method of class C coerce
 the C object into a B object. Then call the next method in the list with the
 B class object. Now, in the initialize method of class B the object is a B
 object and the respective generateSpec method is called. Then, in the
 initialize method of C the returned object from callNextMethod has to be
 written to the C class object in .Object. See the code below.
 
 setMethod(initialize, C, function(.Object, value) {.Object@c - value;
 object - as(.Object, B); object - callNextMethod(object, value);
 as(.Object, B) - object; .Object - generateSpec(.Object);
 return(.Object)})
 
 This setting works. I do not know though, if this setting is the usual way
 such things are done in R OOP. Maybe the whole class design is
 disadvantageous. If anyone detects a mistaken design, I am very thankful to
 learn.
 
 Hi Simon -- your 'simple' example is pretty complicated, and I didn't really 
 follow it in detail! The code is not formatted for easy reading (e.g., lines 
 spanning no more than 80 columns) and some of it (e.g., generateSpec) might 
 not be necessary to describe the problem you're having.
 
 A good strategy is to ensure that 'new' called with no arguments works (there 
 are other solutions, but following this rule has helped me to keep my classes 
 and methods simple). This is not the case for
 
  new(A)
  new(C)
 
 The reason for this strategy has to do with the way inheritance is 
 implemented, in particular the coercion from derived to super class. Usually 
 it is better to provide default values for arguments to initialize, and to 
 specify arguments after a '...'. This means that your initialize methods will 
 respects the contract set out in ?initialize, in particular the handling of 
 unnamed arguments:
 
 ...: data to include in the new object.  Named arguments
  correspond to slots in the class definition. Unnamed
  arguments must be objects from classes that this class
  extends.
 
 I might have written initialize,A-method as
 
  setMethod(initialize, A, function(.Object, ..., value=numeric()){
  .Object - callNextMethod(.Object, ..., a=value)
  generateSpec(.Object)
  })
 
 Likely in a subsequent iteration I would have ended up with (using the 
 convention that function names preceded by '.' are not exported)
 
  .A - setClass(A, representation(a = numeric, specA = numeric))
 
  .generateSpecA - function(a) {
  1 / a
   }
 
  A - function(a=numeric(), ...) {
  specA - .generateSpecA(a)
  .A(..., a=a, specA=specA)
  }
 
  setMethod(generateSpec, A, function(object) {
  .generateSpecA(object@a)
  })
 
 ensuring that A() returns a valid object and avoiding the definition of an 
 initialize method entirely.
 
 Martin
 
 
 Best
 
 Simon
 
 
 On Aug 3, 2013, at 9:43 PM, Simon Zehnder simon.zehn...@googlemail.com
 wrote:
 
 setMethod(initialize, C, function(.Object, value) {.Object@c - value;
 .Object - callNextMethod(.Object, value); .Object -
 generateSpec(.Object); return(.Object)})
 
 __ R-help@r-project.org mailing
 list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting
 guide http://www.R-project.org/posting-guide.html and provide commented,
 minimal, self-contained, reproducible code.
 
 
 
 -- 
 Computational Biology / Fred Hutchinson Cancer Research Center
 1100 Fairview Ave. N.
 PO Box 19024 Seattle, WA 98109
 
 Location: Arnold Building M1 B861
 Phone: (206) 667-2793

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Method dispatch in S4

2013-08-09 Thread Simon Zehnder
Hi Bert,

thank you very much for your suggestion! I will take a look at it soon!

Best

Simon

On Aug 9, 2013, at 4:45 PM, Bert Gunter gunter.ber...@gene.com wrote:

 Simon:
 
 Have a look at the proto package for which there is a vignette. You
 may find it suitable for your needs and less intimidating.
 
 Cheers,
 Bert
 
 On Fri, Aug 9, 2013 at 7:40 AM, Simon Zehnder szehn...@uni-bonn.de wrote:
 Hi Martin,
 
 thank you very much for this profound answer! Your added design advice is 
 very helpful, too!
 
 For the 'simple example': Sometimes I am still a little overwhelmed from a 
 certain setting in the code and my ideas how I want to handle a process. But 
 I learn from session to session. In future I will also span the lines more 
 than 80 columns. I am used to the indent in my vim editor.
 
 I have one further issue: I do know, that you are one of the leading 
 developers of the bioconductor package which uses (as far as I have read) 
 extensively OOP in R. Is there a package you could suggest to me to learn 
 from by reading and understanding the code? Where can I find the source code?
 
 Best
 
 Simon
 
 
 On Aug 8, 2013, at 10:00 PM, Martin Morgan mtmor...@fhcrc.org wrote:
 
 On 08/04/2013 02:13 AM, Simon Zehnder wrote:
 So, I found a solution: First in the initialize method of class C coerce
 the C object into a B object. Then call the next method in the list with 
 the
 B class object. Now, in the initialize method of class B the object is a 
 B
 object and the respective generateSpec method is called. Then, in the
 initialize method of C the returned object from callNextMethod has to 
 be
 written to the C class object in .Object. See the code below.
 
 setMethod(initialize, C, function(.Object, value) {.Object@c - value;
 object - as(.Object, B); object - callNextMethod(object, value);
 as(.Object, B) - object; .Object - generateSpec(.Object);
 return(.Object)})
 
 This setting works. I do not know though, if this setting is the usual 
 way
 such things are done in R OOP. Maybe the whole class design is
 disadvantageous. If anyone detects a mistaken design, I am very thankful to
 learn.
 
 Hi Simon -- your 'simple' example is pretty complicated, and I didn't 
 really follow it in detail! The code is not formatted for easy reading 
 (e.g., lines spanning no more than 80 columns) and some of it (e.g., 
 generateSpec) might not be necessary to describe the problem you're having.
 
 A good strategy is to ensure that 'new' called with no arguments works 
 (there are other solutions, but following this rule has helped me to keep 
 my classes and methods simple). This is not the case for
 
 new(A)
 new(C)
 
 The reason for this strategy has to do with the way inheritance is 
 implemented, in particular the coercion from derived to super class. 
 Usually it is better to provide default values for arguments to initialize, 
 and to specify arguments after a '...'. This means that your initialize 
 methods will respects the contract set out in ?initialize, in particular 
 the handling of unnamed arguments:
 
...: data to include in the new object.  Named arguments
 correspond to slots in the class definition. Unnamed
 arguments must be objects from classes that this class
 extends.
 
 I might have written initialize,A-method as
 
 setMethod(initialize, A, function(.Object, ..., value=numeric()){
 .Object - callNextMethod(.Object, ..., a=value)
 generateSpec(.Object)
 })
 
 Likely in a subsequent iteration I would have ended up with (using the 
 convention that function names preceded by '.' are not exported)
 
 .A - setClass(A, representation(a = numeric, specA = numeric))
 
 .generateSpecA - function(a) {
 1 / a
  }
 
 A - function(a=numeric(), ...) {
 specA - .generateSpecA(a)
 .A(..., a=a, specA=specA)
 }
 
 setMethod(generateSpec, A, function(object) {
 .generateSpecA(object@a)
 })
 
 ensuring that A() returns a valid object and avoiding the definition of an 
 initialize method entirely.
 
 Martin
 
 
 Best
 
 Simon
 
 
 On Aug 3, 2013, at 9:43 PM, Simon Zehnder simon.zehn...@googlemail.com
 wrote:
 
 setMethod(initialize, C, function(.Object, value) {.Object@c - value;
 .Object - callNextMethod(.Object, value); .Object -
 generateSpec(.Object); return(.Object)})
 
 __ R-help@r-project.org mailing
 list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the 
 posting
 guide http://www.R-project.org/posting-guide.html and provide commented,
 minimal, self-contained, reproducible code.
 
 
 
 --
 Computational Biology / Fred Hutchinson Cancer Research Center
 1100 Fairview Ave. N.
 PO Box 19024 Seattle, WA 98109
 
 Location: Arnold Building M1 B861
 Phone: (206) 667-2793
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented

Re: [R] Method dispatch in S4

2013-08-09 Thread Simon Zehnder
Hi Martin,

is proto in S3? 

I will take a look first at the simple package EBImage. 

Thank you very much for the suggestions!

Best

Simon


On Aug 9, 2013, at 5:01 PM, Martin Morgan mtmor...@fhcrc.org wrote:

 On 08/09/2013 07:45 AM, Bert Gunter wrote:
 Simon:
 
 Have a look at the proto package for which there is a vignette. You
 may find it suitable for your needs and less intimidating.
 
 Won't help much with S4, though! Some answers here
 
 http://stackoverflow.com/questions/5437238/which-packages-make-good-use-of-s4-objects
 
 including from Bioconductor simple class in EBImage, the advanced IRanges 
 package and the 'toy' StudentGWAS.
 
 Martin
 
 
 Cheers,
 Bert
 
 On Fri, Aug 9, 2013 at 7:40 AM, Simon Zehnder szehn...@uni-bonn.de wrote:
 Hi Martin,
 
 thank you very much for this profound answer! Your added design advice is 
 very helpful, too!
 
 For the 'simple example': Sometimes I am still a little overwhelmed from a 
 certain setting in the code and my ideas how I want to handle a process. 
 But I learn from session to session. In future I will also span the lines 
 more than 80 columns. I am used to the indent in my vim editor.
 
 I have one further issue: I do know, that you are one of the leading 
 developers of the bioconductor package which uses (as far as I have read) 
 extensively OOP in R. Is there a package you could suggest to me to learn 
 from by reading and understanding the code? Where can I find the source 
 code?
 
 Best
 
 Simon
 
 
 On Aug 8, 2013, at 10:00 PM, Martin Morgan mtmor...@fhcrc.org wrote:
 
 On 08/04/2013 02:13 AM, Simon Zehnder wrote:
 So, I found a solution: First in the initialize method of class C coerce
 the C object into a B object. Then call the next method in the list with 
 the
 B class object. Now, in the initialize method of class B the object is 
 a B
 object and the respective generateSpec method is called. Then, in the
 initialize method of C the returned object from callNextMethod has to 
 be
 written to the C class object in .Object. See the code below.
 
 setMethod(initialize, C, function(.Object, value) {.Object@c - value;
 object - as(.Object, B); object - callNextMethod(object, value);
 as(.Object, B) - object; .Object - generateSpec(.Object);
 return(.Object)})
 
 This setting works. I do not know though, if this setting is the usual 
 way
 such things are done in R OOP. Maybe the whole class design is
 disadvantageous. If anyone detects a mistaken design, I am very thankful 
 to
 learn.
 
 Hi Simon -- your 'simple' example is pretty complicated, and I didn't 
 really follow it in detail! The code is not formatted for easy reading 
 (e.g., lines spanning no more than 80 columns) and some of it (e.g., 
 generateSpec) might not be necessary to describe the problem you're having.
 
 A good strategy is to ensure that 'new' called with no arguments works 
 (there are other solutions, but following this rule has helped me to keep 
 my classes and methods simple). This is not the case for
 
  new(A)
  new(C)
 
 The reason for this strategy has to do with the way inheritance is 
 implemented, in particular the coercion from derived to super class. 
 Usually it is better to provide default values for arguments to 
 initialize, and to specify arguments after a '...'. This means that your 
 initialize methods will respects the contract set out in ?initialize, in 
 particular the handling of unnamed arguments:
 
 ...: data to include in the new object.  Named arguments
  correspond to slots in the class definition. Unnamed
  arguments must be objects from classes that this class
  extends.
 
 I might have written initialize,A-method as
 
  setMethod(initialize, A, function(.Object, ..., value=numeric()){
  .Object - callNextMethod(.Object, ..., a=value)
  generateSpec(.Object)
  })
 
 Likely in a subsequent iteration I would have ended up with (using the 
 convention that function names preceded by '.' are not exported)
 
  .A - setClass(A, representation(a = numeric, specA = numeric))
 
  .generateSpecA - function(a) {
  1 / a
   }
 
  A - function(a=numeric(), ...) {
  specA - .generateSpecA(a)
  .A(..., a=a, specA=specA)
  }
 
  setMethod(generateSpec, A, function(object) {
  .generateSpecA(object@a)
  })
 
 ensuring that A() returns a valid object and avoiding the definition of an 
 initialize method entirely.
 
 Martin
 
 
 Best
 
 Simon
 
 
 On Aug 3, 2013, at 9:43 PM, Simon Zehnder simon.zehn...@googlemail.com
 wrote:
 
 setMethod(initialize, C, function(.Object, value) {.Object@c - 
 value;
 .Object - callNextMethod(.Object, value); .Object -
 generateSpec(.Object); return(.Object)})
 
 __ R-help@r-project.org 
 mailing
 list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the 
 posting
 guide http://www.R-project.org/posting-guide.html and provide commented,
 minimal, self-contained, reproducible code.
 
 
 
 --
 Computational Biology

Re: [R] Method dispatch in S4

2013-08-04 Thread Simon Zehnder
So, I found a solution: First in the initialize method of class C coerce the 
C object into a B object. Then call the next method in the list with the B 
class object. Now, in the initialize method of class B the object is a B 
object and the respective generateSpec method is called. Then, in the 
initialize method of C the returned object from callNextMethod has to be 
written to the C class object in .Object. See the code below.

setMethod(initialize, C, function(.Object, value) {.Object@c - value; 
object - as(.Object, B); object - callNextMethod(object, value); 
as(.Object, B) - object; .Object - generateSpec(.Object); return(.Object)})

This setting works. I do not know though, if this setting is the usual way 
such things are done in R OOP. Maybe the whole class design is disadvantageous. 
If anyone detects a mistaken design, I am very thankful to learn.

Best

Simon

 
On Aug 3, 2013, at 9:43 PM, Simon Zehnder simon.zehn...@googlemail.com wrote:

 setMethod(initialize, C, function(.Object, value) {.Object@c - value; 
 .Object - callNextMethod(.Object, value); .Object - generateSpec(.Object); 
 return(.Object)})

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] fortran random number

2013-08-04 Thread Simon Zehnder
Hi Filippo,

I would advice you to use the RcppArmadillo package - maybe in combination with 
the inline package. It is extremely fast and very flexible. A simple RNGScope() 
lets you take the RNG from R. I myself programmed an MCMC sampler with this 
package and it is much faster, than the one I had before in matlab. 
Rcpp gives you a very nice API to communicate with R + you can reuse memory 
from R which makes it even faster. 

You find a lot of documentation online:

http://www.rcpp.org
http://cran.r-project.org/web/packages/RcppArmadillo/index.html
http://cran.r-project.org/web/packages/Rcpp/index.html


Best

Simon


On Aug 4, 2013, at 10:38 PM, Filippo ingf...@gmail.com wrote:

 Hi,
 I need an advice.
 I'm writing a Markov Chain Monte Carlo subroutine in F95 which should be 
 called from R. I noticed that the subroutine RANDOM_NUMBER takes the seed 
 from R (is that correct?).
 I was wondering if is better initialize the random seed inside the F95 
 subroutine or use the R seed.
 Regards,
 Filippo
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Method dispatch in S4

2013-08-03 Thread Simon Zehnder
Dear R-Users and R-Devels,

I am struggling with the method dispatch for S4 objects. First I will give a 
simple example to make clear what the general setting in my project is. 

Three classes will be build in the example: A is just a simple class. B 
contains A in a slot and C inherits from B. In addition all classes have a slot 
called specA, specB, specC that should be computed during initialization 
by calling the method generateSpec defined for each class A,B,C:

## Define class A##
setClass(A, representation(a = numeric, specA = numeric))
setMethod(initialize, A, function(.Object, value){.Object@a - value; 
.Object - generateSpec(.Object); return(.Object)})
setGeneric(generateSpec, function(object) standardGeneric(generateSpec))
setMethod(generateSpec, A, function(object) {object@specA - 1/object@a; 
return(object)})

## Define class B containing A in a slot ##
setClass(B, representation(b = numeric, specB = numeric, A = A))
setMethod(initialize, B, function(.Object, value = 0){.Object@b - value; 
.Object@A - new(A, value = value); .Object - generateSpec(.Object); 
return(.Object)})
setMethod(generateSpec, B, function(object) {aSpec - object@A@specA; 
object@specB - 1/aSpec; return(object)})
## all fine:
new(B, value = 2)

## Define class C inheriting from class B ##
setClass(C, representation(c = numeric, specC = numeric), contains=c(B))
setMethod(initialize, C, function(.Object, value) {.Object@c - value; 
.Object - callNextMethod(.Object, value); .Object - generateSpec(.Object); 
return(.Object)})
setMethod(generateSpec, C, function(object) {bSpec - object@specB; 
object@specC - 1/bSpec; return(object)})
## numeric(0) for specB and specC
new(C, value = 2)

Note, that the generateSpec method of B depends on the specA slot of A, the 
generateSpec method of class C depends on the specB slot of class B. Note 
also, that in the initialize method of B one has to give the value argument 
a default value. When calling the last line of code, specB and specC are 
numeric(0). As the slot b of the new C class is equal to 2, the 
initialize method of class B is called but the .Object in the initialize 
method of class B is an object of class C. And of course in this case method 
dispatching chooses the generateSpec method of class C - generateSpec for 
class B gets never called. In general this is exactly what we want, when 
inheriting classes: dispatching makes OOP programming not only simpler, but 
actually possible. 

My question: How can I make this setting work? How can I generate an 
initialize method for class B that calls ALWAYS the generateSpec method of 
class B (using as lets me end up in the initialize method of class C with 
an object of class B as it is only possible to convert a derived object to its 
base class but not the other way around)? 

Best

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] R number format with Hmisc and knitr

2013-07-31 Thread Simon Zehnder
Dear R-Users and R-Devels,

I have a problem when using knitr in combination with Hmisc. I generate a 
data.frame which has mixed scientific and non-scientific numbers inside. In my 
Latex Table I just want to have non-scientific format, so I call

latex(myDataFrame,
file = '',
cdec = c(0, rep(4, NROW(myDataFrame) - 1)),
)

Usually this works, but in this case it doesn't. I do not know why but suspect 
the mixed data format to be the culprit. What could I do?

Using format(, scientific = FALSE) before or options(scipen = 4) before has no 
influence. 


Best

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R number format with Hmisc and knitr

2013-07-31 Thread Simon Zehnder
Errata:

it must say:

latex(myDataFrame,
file = '',
cdec = c(0, rep(4, NCOL(myDataFrame) - 1))
)

But this does not work. Scientific notation is very robust :)

Apologize

Simon

On Jul 31, 2013, at 5:05 PM, Simon Zehnder szehn...@uni-bonn.de wrote:

 Dear R-Users and R-Devels,
 
 I have a problem when using knitr in combination with Hmisc. I generate a 
 data.frame which has mixed scientific and non-scientific numbers inside. In 
 my Latex Table I just want to have non-scientific format, so I call
 
 latex(myDataFrame,
 file = '',
 cdec = c(0, rep(4, NROW(myDataFrame) - 1)),
 )
 
 Usually this works, but in this case it doesn't. I do not know why but 
 suspect the mixed data format to be the culprit. What could I do?
 
 Using format(, scientific = FALSE) before or options(scipen = 4) before has 
 no influence. 
 
 
 Best
 
 Simon
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] MCMClogit: Cannot calculate marginal likelihood with improper prior

2013-07-29 Thread Simon Zehnder
Hi, 

what I see so far is that you have specified your user.prior.density correctly. 
The error comes from the prior precision matrix B0 in combination with the 
marginal.likelihood set to Laplace. B0, if not explicitly specified, defaults 
to zero, which results in eigenvalues of zero. If Laplace is indicated for 
the marginal.likelihood, the algorithm usually calls an optimization over 
logpost.logit in BayesianFactors.R where the matrix B0 is tried to be 
solved by solve(B0) ... as it is a zero matrix its linear equation system is 
exactly singular and cannot be solved. The Function MCMClogit knows about this 
fact and gives out a warning Cannot calculate marginal likelihood with 
improper prior while changing marginal.likelihood to none.

So concluding: Choose your user.prior.density with marginal.likelihood = none 
and all is fine (implicitly it is done so nevertheless).

Best

Simon

P.S. Using a name on a community help list will certainly improve the number of 
answers to your questions. 

 
On Jul 29, 2013, at 3:00 AM, ba0728 haleyb...@att.net wrote:

 I'm an undergrad who is new to MCMCpack and I haven't been able to find an
 answer to my problem online yet: I'm attempting to run MCMClogit with a
 Cauchy proper prior but I'm getting the warning Cannot calculate marginal
 likelihood with improper prior (my purposes require the marginal likelihood
 calculation so I understand that I need to use a proper prior).
 
 I'm trying to simulate the user-defined independent Cauchy prior with
 additional args as specified in the MCMCpack User Manual (p. 76, April 2013
 version). My input data has been standardized  (mean = 0, sd = 0.5 for
 non-binary variables, and binary variables with mean of 0 and difference of
 1 between upper and lower ends) according to the Gelman 2008 paper on
 logistic regression
 (www.stat.columbia.edu/~gelman/research/published/priors11.pdf‎). 
 
 When I run the example data set (birthwt) from the User Manual, the
 logpriorfun works correctly allowing the marginal likelihood to be
 generated. However, when I try running my data with the logprior fun, I get
 a warning that the prior is improper. Here is the code I am running:
 
 *logpriorfun = function(beta, location,scale){
  sum(dcauchy(beta, location, scale, log = TRUE))
 }*
 
 * MCMC.2= MCMClogit(DEAD ~ YEARS + MALE + x1 + x2 + x3+ x4 +x5 + x6 + x7 +
 x8 + x9, tune= 0.65,burnin =500, mcmc=5000, data = dat, marginal.likelihood
 = Laplace, user.prior.density=logpriorfun, logfun=TRUE, location = 0,
 scale=2.5)
 *
 
 *@
 The Metropolis acceptance rate was 0.27418
 @
 Warning message:
 In MCMClogit(DEAD ~ YEARS + MALE + x1 + x2 + x3 +  :
  Cannot calculate marginal likelihood with improper prior*
 
 Any advice on how to fix my arguments so it is a proper prior and will allow
 me to generate a marginal likelihood using the Laplace approximation? Or how
 should I be coding a Cauchy proper prior? I'm having problems defining the
 priors.
 
 Thanks, B.
 
 
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/MCMClogit-Cannot-calculate-marginal-likelihood-with-improper-prior-tp4672561.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Hmisc ctable rotate option obsolete?

2013-07-27 Thread Simon Zehnder
So, I downloaded the source files of Hmisc and changed in the file latex.s line 
688 'rotate' to 'sideways'. This does the work for landscape ctables in Latex. 

I also wrote an email to the package maintainer. I consider this thread as 
solved. 


Best

Simon

On Jul 26, 2013, at 5:34 PM, Simon Zehnder szehn...@uni-bonn.de wrote:

 Dear R-Users and R-Devels,
 
 I may have found a deprecated option for the 'latex' function in the Hmisc 
 package. I am working with Hmisc and knitr and tried the following code:
 
 \documentclass{article}
 \usepackage[utf8]{inputenc}
 \usepackage{amsmath, amssymb}
 \usepackage{ctable}
 %\usepackage{booktabs}
 \begin{document}
 results = 'asis'=
 library(Hmisc)
 iris.t - head(iris)
 iris.t[seq(2, NROW(iris.t), by = 2),] - format(iris.t[seq(2, NROW(iris.t), 
 by = 2),], scientific = TRUE)
 texMat - matrix(, ncol = 5, nrow = 6)
 texMat[seq(2,nrow(texMat), by = 2), ] - scriptsize
 latex(iris.t, 
 file = '', 
 landscape = TRUE,
 dcolumn = TRUE,
 col.just = c('r','c', 'r','c', 'l'),
 cdec = c(0, 0, 1, 1, 0),
 na.blank = TRUE,
 rowname = '',
 rowlabel = '', 
 cellTexCmd = texMat,
 ctable = TRUE, 
 cgroup = c('Observations', ''),
 n.cgroup = c(4, 1),
 rgroup = c('',''),
 n.rgroup = c(3, 3),
 caption = 'iris'
 )
 @
 \end{document}
 
 Everything runs fine but the 'landscape' option. It says in the help for 
 'latex' that if option 'ctable' is set to TRUE the 'rotate' option for ctable 
 is used if 'nadscape' is set TRUE. Looking at the ctable documentary 
 (http://texdoc.net/texmf-dist/doc/latex/ctable/ctable.pdf) in section Change 
 History, I get for version v1.07: General: Added option sideways, option 
 rotate now obsolete. Hasn't this been updated in the Hmisc package?
 
 Best
 
 Simon
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Hmisc ctable rotate option obsolete?

2013-07-26 Thread Simon Zehnder
Dear R-Users and R-Devels,

I may have found a deprecated option for the 'latex' function in the Hmisc 
package. I am working with Hmisc and knitr and tried the following code:

 \documentclass{article}
\usepackage[utf8]{inputenc}
\usepackage{amsmath, amssymb}
\usepackage{ctable}
%\usepackage{booktabs}
\begin{document}
results = 'asis'=
library(Hmisc)
iris.t - head(iris)
iris.t[seq(2, NROW(iris.t), by = 2),] - format(iris.t[seq(2, NROW(iris.t), by 
= 2),], scientific = TRUE)
texMat - matrix(, ncol = 5, nrow = 6)
texMat[seq(2,nrow(texMat), by = 2), ] - scriptsize
latex(iris.t, 
file = '', 
landscape = TRUE,
dcolumn = TRUE,
col.just = c('r','c', 'r','c', 'l'),
cdec = c(0, 0, 1, 1, 0),
na.blank = TRUE,
rowname = '',
rowlabel = '', 
cellTexCmd = texMat,
ctable = TRUE, 
cgroup = c('Observations', ''),
n.cgroup = c(4, 1),
rgroup = c('',''),
n.rgroup = c(3, 3),
caption = 'iris'
)
@
\end{document}

Everything runs fine but the 'landscape' option. It says in the help for 
'latex' that if option 'ctable' is set to TRUE the 'rotate' option for ctable 
is used if 'nadscape' is set TRUE. Looking at the ctable documentary 
(http://texdoc.net/texmf-dist/doc/latex/ctable/ctable.pdf) in section Change 
History, I get for version v1.07: General: Added option sideways, option rotate 
now obsolete. Hasn't this been updated in the Hmisc package?

Best

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Check the class of an object

2013-07-23 Thread Simon Zehnder
Dear R-Users and R-Devels,

I have large project based on S4 classes. While writing my unit tests I found 
out, that 'is' cannot test for a specific class, as also inherited classes can 
be treated as their super classes. I need to do checks for specific classes. 
What I do right now is sth. like

if (class(myClass) == firstClass) {

} else if (class(myClass) == secondClass) {

}

Is this the usual way how classes are checked in R? I was expecting some 
specific method (and 'inherits' or 'extends' is not what I look for)...


Best

Simon





[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Check the class of an object

2013-07-23 Thread Simon Zehnder
Hi David,

thanks for the reply. You are right. Using the %in% is more stable and I gonna 
change my code.

When testing for a specific class using 'is' one has to start at the lowest 
heir and walk up the inheritance structure. Starting at the checks at the root 
will always give TRUE. Having a structure which is quite complicated let me 
move to the check I suggested in my first mail. 

Best

Simon

On Jul 23, 2013, at 6:15 PM, David Winsemius dwinsem...@comcast.net wrote:

 
 On Jul 23, 2013, at 5:36 AM, Simon Zehnder wrote:
 
 Dear R-Users and R-Devels,
 
 I have large project based on S4 classes. While writing my unit tests I 
 found out, that 'is' cannot test for a specific class, as also inherited 
 classes can be treated as their super classes. I need to do checks for 
 specific classes. What I do right now is sth. like
 
 if (class(myClass) == firstClass) {
 
 I would think that you would need to use `%in%` instead. 
 
 if( firstClass %in% class(myObject) ){
 
 Objects can have more than one class, so testing with == would fail in 
 those instances.
 
 
 
 } else if (class(myClass) == secondClass) {
 
 }
 
 Is this the usual way how classes are checked in R?
 
 Well, `inherits` IS the usual way.
 
 I was expecting some specific method (and 'inherits' or 'extends' is not 
 what I look for)...
 
 
 Best
 
 Simon
 
  [[alternative HTML version deleted]]
 
 Plain-text format is the recommended format for Rhelp
 
 -- 
 David Winsemius
 Alameda, CA, USA
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Check the class of an object

2013-07-23 Thread Simon Zehnder
Hi Martin,

I didn't know that. But that is even more comfortable for checking. Thanks for 
the quick update!

Best


Simon

P.S. And thanks for the online documents about S4 - I could already learn a lot!


On Jul 23, 2013, at 7:11 PM, Martin Morgan mtmor...@fhcrc.org wrote:

 On 07/23/2013 09:59 AM, Simon Zehnder wrote:
 Hi David,
 
 thanks for the reply. You are right. Using the %in% is more stable and I 
 gonna change my code.
 
 you said you were you were using S4 classes. S4 classes do not report vectors 
 of length != 1, from ?class
 
 For objects which have a formal class, its name is returned by 'class'
 as a character vector of length one
 
 so a first unit test could be
 
  stopifnot(length(class(myObject)) != 1L)
 
 
 
 When testing for a specific class using 'is' one has to start at the lowest 
 heir and walk up the inheritance structure. Starting at the checks at the 
 root will always give TRUE. Having a structure which is quite complicated 
 let me move to the check I suggested in my first mail.
 
 Best
 
 Simon
 
 On Jul 23, 2013, at 6:15 PM, David Winsemius dwinsem...@comcast.net wrote:
 
 
 On Jul 23, 2013, at 5:36 AM, Simon Zehnder wrote:
 
 Dear R-Users and R-Devels,
 
 I have large project based on S4 classes. While writing my unit tests I 
 found out, that 'is' cannot test for a specific class, as also inherited 
 classes can be treated as their super classes. I need to do checks for 
 specific classes. What I do right now is sth. like
 
 if (class(myClass) == firstClass) {
 
 I would think that you would need to use `%in%` instead.
 
 if( firstClass %in% class(myObject) ){
 
 Objects can have more than one class, so testing with == would fail in 
 those instances.
 
 
 
 } else if (class(myClass) == secondClass) {
 
 }
 
 Is this the usual way how classes are checked in R?
 
 Well, `inherits` IS the usual way.
 
 I was expecting some specific method (and 'inherits' or 'extends' is not 
 what I look for)...
 
 
 Best
 
 Simon
 
[[alternative HTML version deleted]]
 
 Plain-text format is the recommended format for Rhelp
 
 --
 David Winsemius
 Alameda, CA, USA
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 
 
 -- 
 Computational Biology / Fred Hutchinson Cancer Research Center
 1100 Fairview Ave. N.
 PO Box 19024 Seattle, WA 98109
 
 Location: Arnold Building M1 B861
 Phone: (206) 667-2793

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Check the class of an object

2013-07-23 Thread Simon Zehnder
Hi Hervé,

thank you very much for your reply! This makes the different treatment of S3 
and S4 objects by 'class' clear.

Best

Simon


On Jul 23, 2013, at 7:20 PM, Hervé Pagès hpa...@fhcrc.org wrote:

 Hi,
 
 On 07/23/2013 09:59 AM, Simon Zehnder wrote:
 Hi David,
 
 thanks for the reply. You are right. Using the %in% is more stable and I 
 gonna change my code.
 
 Unlike with S3 objects, class() on an S4 object can only return 1 class.
 
 Also note that, on an S3 object, doing
 
  firstClass %in% class(myObject)
 
 is equivalent to doing inherits(myObject, firstClass), which is
 what you said you wanted to avoid. The most specific class should be
 the first so if that's what you wanted to check, you could do
 
  class(myObject)[1] == firstClass
 
 But that precaution is not needed if 'myObject' is guaranteed to be
 an S4 object (although when writing a unit test, one should probably
 discard any guarantee of that sort).
 
 Cheers,
 H.
 
 
 
 When testing for a specific class using 'is' one has to start at the lowest 
 heir and walk up the inheritance structure. Starting at the checks at the 
 root will always give TRUE. Having a structure which is quite complicated 
 let me move to the check I suggested in my first mail.
 
 Best
 
 Simon
 
 On Jul 23, 2013, at 6:15 PM, David Winsemius dwinsem...@comcast.net wrote:
 
 
 On Jul 23, 2013, at 5:36 AM, Simon Zehnder wrote:
 
 Dear R-Users and R-Devels,
 
 I have large project based on S4 classes. While writing my unit tests I 
 found out, that 'is' cannot test for a specific class, as also inherited 
 classes can be treated as their super classes. I need to do checks for 
 specific classes. What I do right now is sth. like
 
 if (class(myClass) == firstClass) {
 
 I would think that you would need to use `%in%` instead.
 
 if( firstClass %in% class(myObject) ){
 
 Objects can have more than one class, so testing with == would fail in 
 those instances.
 
 
 
 } else if (class(myClass) == secondClass) {
 
 }
 
 Is this the usual way how classes are checked in R?
 
 Well, `inherits` IS the usual way.
 
 I was expecting some specific method (and 'inherits' or 'extends' is not 
 what I look for)...
 
 
 Best
 
 Simon
 
[[alternative HTML version deleted]]
 
 Plain-text format is the recommended format for Rhelp
 
 --
 David Winsemius
 Alameda, CA, USA
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 
 -- 
 Hervé Pagès
 
 Program in Computational Biology
 Division of Public Health Sciences
 Fred Hutchinson Cancer Research Center
 1100 Fairview Ave. N, M1-B514
 P.O. Box 19024
 Seattle, WA 98109-1024
 
 E-mail: hpa...@fhcrc.org
 Phone:  (206) 667-5791
 Fax:(206) 667-1319

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] S4 method signature - integer matrix

2013-07-20 Thread Simon Zehnder
Hi David,

thanks for the reply! 

The prototype solution is I think not the appropriate solution, when the method 
adds the integer matrix next to the class object into the argument list. Also 
the .validObject method can only be applied to the object for which the method 
is called or I have to write an integer matrix class that has its own 
.validObject method.

From your suggestions I had a direction to search more precisely and I found 
this: 

Set the signature to matrix: checks for matrix when method is called.
Then check for typeof(mymatrix) == integer: Return an exception if FALSE.
In case of an index I can add then an all(apply(mymatrix, c(1, 2), , 0)): 
Return an exception if FALSE

Thank you very much for your advice. It helped me to find a good solution.


Best 

Simon


On Jul 20, 2013, at 12:03 AM, David Winsemius dwinsem...@comcast.net wrote:

 
 On Jul 19, 2013, at 9:54 AM, Simon Zehnder wrote:
 
 Dear R-Users and R-Devels,
 
 I am programming a package with S4 classes and I search for a solution of 
 the following problem: 
 
 If you want an S4 method to await an integer argument you set the signature 
 like this
 
 setMethod(myfunction, signature(object = myClass, y = integer), 
 function(object, y) {//Do sth})
 
 Now, if you want the method to await an integer matrix or array there is 
 only one way how to define it
 
 setMethod(myfunction, signature(object = myClass, y = matrix), 
 function(object, y) {//Do sth}) or
 setMethod(myfunction, signature(object = myClass, y = array), 
 function(object, y) {//Do sth}) 
 
 am I right?
 
 I don't think that is the only way. You could also test for mode being 
 'numeric' at the signature level and then within the validation check to see 
 whether it has dimensions. See ?is.numeric and ?validObject
 
 There is also the option of using a prototype.
 
 As I understand the S4 checks they will apply an is.mode or is.class test 
 as long as there is an correspond function that returns a logical. see ... ?is
 
 
 WIth this you can also pass a numeric matrix.
 
 How would you check for an integer matrix in an S4 method (in the index 
 function of R I think it is just coerced to integer)?
 
 It would need be both integer mode...   ?is.integer and have dimensions...  
 ?dim
 
 Furthermore: How would you check for an integer matrix with values 
 bigger one (so the typical R indices)?
 
 Just add a test in the 'validity' code. See ... ?setClass
 
 Is there a way how it is usually done in R (I think probably with apply())? 
 Or is it usually better to throw an exception?
 
 
 I don't see a reason to use apply.
 
 -- 
 David Winsemius
 Alameda, CA, USA
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] S4 method signature - integer matrix

2013-07-19 Thread Simon Zehnder
Dear R-Users and R-Devels,

I am programming a package with S4 classes and I search for a solution of the 
following problem: 

If you want an S4 method to await an integer argument you set the signature 
like this

setMethod(myfunction, signature(object = myClass, y = integer), 
function(object, y) {//Do sth})

Now, if you want the method to await an integer matrix or array there is only 
one way how to define it

setMethod(myfunction, signature(object = myClass, y = matrix), 
function(object, y) {//Do sth}) or
setMethod(myfunction, signature(object = myClass, y = array), 
function(object, y) {//Do sth}) 

am I right? WIth this you can also pass a numeric matrix.

How would you check for an integer matrix in an S4 method (in the index 
function of R I think it is just coerced to integer)? Furthermore: How would 
you check for an integer matrix with values 
bigger one (so the typical R indices)? Is there a way how it is usually done in 
R (I think probably with apply())? Or is it usually better to throw an 
exception?


Best

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] 'save' method for S4 class

2013-07-18 Thread Simon Zehnder
Hi Christopher,

I think, that save is no generic function like plot, show, etc. So at 
first you have to determine a generic.

setGeneric(save, function(x, file_Path) standardGeneric(save))

Now your definition via setMethod. 


Best

Simon



On Jul 18, 2013, at 12:09 PM, Christofer Bogaso bogaso.christo...@gmail.com 
wrote:

 Hello again,
 
 I am trying to define the 'save' method for my S4 class as below:
 
 setClass(MyClass, representation(
   Slot1 = data.frame
   ))  
   
 setMethod(save, MyClass, definition = function(x, file_Path) {
   
   write.table(x@Slot1, file = file_Path, append = FALSE, quote = 
 TRUE,
 sep = ,,
   eol = \n, na = NA, dec = 
 ., row.names = FALSE,
   col.names = TRUE, qmethod = 
 c(escape, double),
   fileEncoding = )
   })
 
 However while doing this I am getting following error:
 
 Error in conformMethod(signature, mnames, fnames, f, fdef, definition) :
  in method for ‘save’ with signature ‘list=MyClass’: formal
 arguments (list = MyClass, file = MyClass, ascii = MyClass,
 version = MyClass, envir = MyClass, compress = MyClass,
 compression_level = MyClass, eval.promises = MyClass, precheck =
 MyClass) omitted in the method definition cannot be in the signature
 
 
 Can somebody point me what will be the correct approach to define
 'save' method for S4 class?
 
 Thanks and regards,
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] 'save' method for S4 class

2013-07-18 Thread Simon Zehnder
Hi Christofer,

R's S4 class system builds on generics, that is one class defines the raw frame 
of a function and maybe an own distinct implementation and other classes, that 
(usually) inherit from this class, set their special implementations of it. 
Consider a generic function in R as being an abstract function, just telling 
inheriting classes how to build such a function, i.e. what to put in. 
Furthermore it enables the dispatches like callNextMethod, which calls the 
method from the parent (be aware if you have more than one direct parent!). 
Under the hood, R must know from the input arguments what function should be 
called. If you put in an B class that inherits from an A class it must know, 
what to do with a user-defined save function when the input argument is for 
example an object of class B. If you put in an A class object, it should 
know how to proceed this object and so on. So the generic method gives a kind 
of internal map for R how to call this function in which contexts.

The standardGeneric() function: We must provide a name and a definition of the 
generic function. In most cases you just want to tell, what arguments the 
function should expect. With the function standardGeneric(name) R creates a 
generic which can then be dispatched in form of distinct methods. So if you 
call save R sees, that it is a generic function and checks which of the 
distinct functions to call, evaluating the arguments. If it finds no specific 
function for the arguments at hand, it calls a default function. If, as in your 
case, the function name is already taken, R sets, when you use setGeneric, the 
usual save function as the default method, getting called when nothing else 
fits ... resulting in either an error because of mismatching function arguments 
or in saving in a way you do not want to. Solution is here to define inside the 
setGeneric also the default function via the argument useAsDefault = function() 
{}. 

So generic function (setGeneric): What should be done?
Distinct function (setMethod): How should it be done?


You find a very good introduction to the S4 class system here: 
http://cran.r-project.org/doc/contrib/Genolini-S4tutorialV0-5en.pdf

Deep down information can be found for example in the books: Programming with 
Data and Software for Data Analysis and Computing with R, both from John 
Chambers. 


Best

Simon



On Jul 18, 2013, at 8:23 PM, Christofer Bogaso bogaso.christo...@gmail.com 
wrote:

 Hi Simon,
 
 Thanks for your pointer.
 
 However could you please explain what 'function(x, file_Path) 
 standardGeneric(save)' means here? The underlying help files look quite 
 rocket science for me!
 
 Thanks for your time.
 
 Thanks and regards,
 
 
 On Thu, Jul 18, 2013 at 4:32 PM, Simon Zehnder szehn...@uni-bonn.de wrote:
 Hi Christopher,
 
 I think, that save is no generic function like plot, show, etc. So at 
 first you have to determine a generic.
 
 setGeneric(save, function(x, file_Path) standardGeneric(save))
 
 Now your definition via setMethod.
 
 
 Best
 
 Simon
 
 
 
 On Jul 18, 2013, at 12:09 PM, Christofer Bogaso bogaso.christo...@gmail.com 
 wrote:
 
  Hello again,
 
  I am trying to define the 'save' method for my S4 class as below:
 
  setClass(MyClass, representation(
Slot1 = data.frame
))
 
  setMethod(save, MyClass, definition = function(x, file_Path) {
 
write.table(x@Slot1, file = file_Path, append = FALSE, quote 
  = TRUE,
  sep = ,,
eol = \n, na = NA, dec = 
  ., row.names = FALSE,
col.names = TRUE, qmethod = 
  c(escape, double),
fileEncoding = )
})
 
  However while doing this I am getting following error:
 
  Error in conformMethod(signature, mnames, fnames, f, fdef, definition) :
   in method for ‘save’ with signature ‘list=MyClass’: formal
  arguments (list = MyClass, file = MyClass, ascii = MyClass,
  version = MyClass, envir = MyClass, compress = MyClass,
  compression_level = MyClass, eval.promises = MyClass, precheck =
  MyClass) omitted in the method definition cannot be in the signature
 
 
  Can somebody point me what will be the correct approach to define
  'save' method for S4 class?
 
  Thanks and regards,
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Looking for knitr example for beginner (NO RStudio)

2013-07-18 Thread Simon Zehnder
Hi Mike,

I found my way with this little blog: http://yihui.name/knitr/demo/editors/

The .Rnw files are created very well in a Latex editor. Everything else can be 
easily googled. The command via knitr::knit2pdf works very fine if you use the 
chunks. If you are trying to compile an Rtex file, this I do not know either (I 
like the symbols though in for example 
https://github.com/yihui/knitr-examples/blob/master/005-latex.Rtex). But the 
.Rnw files are compiled pretty nice in e.g. texmaker, as described in the blog. 
Use for example this source file: 
https://github.com/yihui/knitr/blob/master/inst/examples/knitr-minimal.Rnw


Hope this helps


Best

Simon





On Jul 18, 2013, at 8:52 PM, C W tmrs...@gmail.com wrote:

 How do you create a .Rnw file, in R or LaTex?  I don't think any
 tutorial mentions it.
 
 btw, I am very new to the terms like markdown, so I don't understand
 markdown to HTML.
 
 I am reading here http://biostat.mc.vanderbilt.edu/wiki/Main/KnitrHowto
 that you need to compile at terminal.  I do not know terminal, is
 there other ways?
 
 Could you do a video on just simple R?  I have seen 3 videos on R
 Studio including yours.
 
 Mike
 
 On Thu, Jul 18, 2013 at 2:43 PM, Yihui Xie x...@yihui.name wrote:
 I'm not sure what your question really is. You do not have to use
 RStudio, but it will be much easier to get started with RStudio,
 because it does a lot of automatic conversion behind the scenes (e.g.
 tex to PDF, markdown to HTML, ...). If you want a pure solution
 without any text editor support, the answer is
 
 library(knitr)
 knit('your_input_file')
 
 For example, knit('foo.Rnw') gives you foo.tex; if you are familiar
 with LaTeX, you can mess with this foo.tex now (outside of R).
 
 Minimal examples for different document formats are at
 http://yihui.name/knitr/demo/minimal/ (you must have read this page),
 and more examples at https://github.com/yihui/knitr-examples
 
 If you are asking about the internals of knitr, Luke, use the
 source: https://github.com/yihui/knitr Or for a more comprehensive
 introduction, see http://www.crcpress.com/product/isbn/9781482203530
 
 Regards,
 Yihui
 --
 Yihui Xie xieyi...@gmail.com
 Phone: 206-667-4385 Web: http://yihui.name
 Fred Hutchinson Cancer Research Center, Seattle
 
 
 On Thu, Jul 18, 2013 at 11:13 AM, C W tmrs...@gmail.com wrote:
 Hi everyone,
 
 I am using package knitr, FIRST TIME.  I don't have access to RStudio.
 
 Read through Yihui's page, didn't find it helpful.  Stuck on terms
 Rnw, GFM (GitHub Flavored Markdown).  Never used Sweave, so the
 reference is not helping.
 
 Is there a simple step-by-step example WITHOUT RStudio?
 
 My question:
 What is the procedure?  The documentation explains the functions, but
 does not say how to operate between R and LaTex.
 
 Mike
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Looking for knitr example for beginner (NO RStudio)

2013-07-18 Thread Simon Zehnder
The executable is in case of knitr always Rscript. 
On a mac it is simply Rscript on windows it is Rscript.exe. This should be on 
your PATH. If you are not sure, open the Mac Terminal and type Rscript 
--version. If it does not say Command not found all is fine. 

Best

Simon
On Jul 18, 2013, at 10:09 PM, C W tmrs...@gmail.com wrote:

 http://tex.stackexchange.com/questions/85154/knitr-with-texworks/85165#85165
 
 In step 3: add the executable file (step 3).
 
 What is the executable file?  Locate package knitr directory path in R?
 
 Mike
 
 On Thu, Jul 18, 2013 at 3:56 PM, C W tmrs...@gmail.com wrote:
 Actually, I see it at the bottom.  Sorry!
 
 On Thu, Jul 18, 2013 at 3:44 PM, C W tmrs...@gmail.com wrote:
 Hi Simon,
 I am on OS X Lion, I have TeXworks, I don't have knitr as an option.
 
 How do I install that into TeXworks?  Seems like I have to something
 in terminal?
 
 Mike
 
 On Thu, Jul 18, 2013 at 3:31 PM, Simon Zehnder szehn...@uni-bonn.de wrote:
 Hi Mike,
 
 I found my way with this little blog: http://yihui.name/knitr/demo/editors/
 
 The .Rnw files are created very well in a Latex editor. Everything else 
 can be easily googled. The command via knitr::knit2pdf works very fine if 
 you use the chunks. If you are trying to compile an Rtex file, this I do 
 not know either (I like the symbols though in for example 
 https://github.com/yihui/knitr-examples/blob/master/005-latex.Rtex). But 
 the .Rnw files are compiled pretty nice in e.g. texmaker, as described in 
 the blog. Use for example this source file: 
 https://github.com/yihui/knitr/blob/master/inst/examples/knitr-minimal.Rnw
 
 
 Hope this helps
 
 
 Best
 
 Simon
 
 
 
 
 
 On Jul 18, 2013, at 8:52 PM, C W tmrs...@gmail.com wrote:
 
 How do you create a .Rnw file, in R or LaTex?  I don't think any
 tutorial mentions it.
 
 btw, I am very new to the terms like markdown, so I don't understand
 markdown to HTML.
 
 I am reading here http://biostat.mc.vanderbilt.edu/wiki/Main/KnitrHowto
 that you need to compile at terminal.  I do not know terminal, is
 there other ways?
 
 Could you do a video on just simple R?  I have seen 3 videos on R
 Studio including yours.
 
 Mike
 
 On Thu, Jul 18, 2013 at 2:43 PM, Yihui Xie x...@yihui.name wrote:
 I'm not sure what your question really is. You do not have to use
 RStudio, but it will be much easier to get started with RStudio,
 because it does a lot of automatic conversion behind the scenes (e.g.
 tex to PDF, markdown to HTML, ...). If you want a pure solution
 without any text editor support, the answer is
 
 library(knitr)
 knit('your_input_file')
 
 For example, knit('foo.Rnw') gives you foo.tex; if you are familiar
 with LaTeX, you can mess with this foo.tex now (outside of R).
 
 Minimal examples for different document formats are at
 http://yihui.name/knitr/demo/minimal/ (you must have read this page),
 and more examples at https://github.com/yihui/knitr-examples
 
 If you are asking about the internals of knitr, Luke, use the
 source: https://github.com/yihui/knitr Or for a more comprehensive
 introduction, see http://www.crcpress.com/product/isbn/9781482203530
 
 Regards,
 Yihui
 --
 Yihui Xie xieyi...@gmail.com
 Phone: 206-667-4385 Web: http://yihui.name
 Fred Hutchinson Cancer Research Center, Seattle
 
 
 On Thu, Jul 18, 2013 at 11:13 AM, C W tmrs...@gmail.com wrote:
 Hi everyone,
 
 I am using package knitr, FIRST TIME.  I don't have access to RStudio.
 
 Read through Yihui's page, didn't find it helpful.  Stuck on terms
 Rnw, GFM (GitHub Flavored Markdown).  Never used Sweave, so the
 reference is not helping.
 
 Is there a simple step-by-step example WITHOUT RStudio?
 
 My question:
 What is the procedure?  The documentation explains the functions, but
 does not say how to operate between R and LaTex.
 
 Mike
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Looking for knitr example for beginner (NO RStudio)

2013-07-18 Thread Simon Zehnder
Hi Mike,

if you browse the folders, you find always the Rscript binary (the executable) 
under /Library/Frameworks/R.framework/Versions/.../Resources/Rscript. 

Do not forget to give your tex file the extension .Rnw! Then surround each 
Rcode with write here a nameadd later further options (important one; 
results = 'asis')= Here your r code as you do it in the R shell  at the 
end a @. Always inside the \begin{document} \end{document} tags.



Best

Simon

On Jul 18, 2013, at 10:22 PM, C W tmrs...@gmail.com wrote:

 Thanks, Simon.  I would never figured it out!
 
 I apologize if I sound frustrated, because I am.
 
 @package author: you have a great package, but I think a lot of the
 directions are hand waving.  For the newbies, this leads to more
 confusion.
 
 @Berend: I am using OS X.
 
 Mike
 
 
 On Thu, Jul 18, 2013 at 4:19 PM, Berend Hasselman b...@xs4all.nl wrote:
 
 On 18-07-2013, at 22:09, C W tmrs...@gmail.com wrote:
 
 http://tex.stackexchange.com/questions/85154/knitr-with-texworks/85165#85165
 
 In step 3: add the executable file (step 3).
 
 What is the executable file?  Locate package knitr directory path in R?
 
 
 From the window:  Executable == Program. So the executable is Rscript.exe.
 
 
 Berend
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Serialize data.frame to database

2013-07-16 Thread Simon Zehnder
Hi Jeff,

I think you are more right than me. I have not much experience with R-specific 
objects and Databases. I just know, that for languages like Java, C, etc. this 
process is a usual one for storing for example matrices (in case of course that 
I do not have to access specific elements inside of the matrix but only the 
matrix as a whole).

Thank you for the tip with the R-sig-DB list. I switch over to this list.

Best

Simon

On Jul 16, 2013, at 2:34 AM, Jeff Newmiller jdnew...@dcn.davis.ca.us wrote:

 I could be wrong, but I would guess that doing what you are describing is 
 very unusual. Most of the time the data frame is mapped to a table in the 
 database so the rows can be searched. Storing data frames as BLOBs really 
 seems odd.
 
 Note that there is an R-sig-db mailing list for questions of this type.
 ---
 Jeff NewmillerThe .   .  Go Live...
 DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
 Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
 /Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
 --- 
 Sent from my phone. Please excuse my brevity.
 
 Simon Zehnder szehn...@uni-bonn.de wrote:
 
 Dear R-Users,
 
 I need a very fast and reliable database solution so I try to serialize
 a data.frame (to binary data) and to store this data to an SQLite
 database. 
 
 This is what I tried to do:
 
 library(RSQLite)
 con - dbDriver(SQLite)
 db - dbConnect(con, test)
 dbSendQuery(db, 'CREATE TABLE frames(simID INT, data BLOB)')
 data.bin - serialize(iris, NULL, ascii = FALSE)
 dbSendQuery(db, paste(INSERT INTO frames VALUES(1, X', data.bin,
 '), sep = ))
 data.bin2 - dbGetQuery(db, SELECT DATA FROM frames WHERE simID = 1)
 data.bin2
 data
 1   58
 
 So, only the first entry of data.bin is saved to the database. I tried
 to first convert the binary data to raw data:
 data.raw - rawToChar(data.bin)
 Error in rawToChar(data.bin) :
 embedded nul in string:
 'X\n\0\0\0\002\0\003\0\001\0\002\003\0\0\0\003\023\0\0\0\005\0\0\0\016\0\0\0\x96@\024ff@\023\x99\x99\x99\x99\x99\x9a@\022\xcc\xcc\xcc\xcc\xcc\xcd@\022ff@\024\0\0\0\0\0\0@\025\x99\x99\x99\x99\x99\x9a@\022ff@\024\0\0\0\0\0\0@\021\x99\x99\x99\x99\x99\x9a@\023\x99\x99\x99\x99\x99\x9a@\025\x99\x99\x99\x99\x99\x9a@\02333@\02333@\02133@\02733@\026\xcc\xcc\xcc\xcc\xcc\xcd@\025\x99\x99\x99\x99\x99\x9a@\024ff@\026\xcc\xcc\xcc\xcc\xcc\xcd@\024ff@\025\x99\x99\x99\x99\x99\x9a@\024ff@\022ff@\024ff@\02333@\024\0\0\0\0\0\0@\024\0\0\0\0\0\0@\024\xcc\xcc\xcc\xcc\xcc\xcd@\024\xcc\xcc\xcc\xcc\xcc\xcd@\022\xcc\xcc\xcc\xcc\xcc\xcd@\02333@\025\x99\x99\x99\x99\x99\x9a@\024\xcc\xcc\xcc\xcc\xcc\xcd@\026\0\0\0\0\0\0@\023\x99\x99\x99\x99\x99\x9a@\024\0\0\0\0\0\0@\026\0\0\0\0\0\0@\023\x99\x99\x99\x99\x99\x9a@\021\x99\x99\x99\x99\x99\x9a@\024ff@\024\0\0\0\0\0\0@\022\0\0\0\0\0\0@\021\x99\x99\x99\x99\x99\x9a@\024\0\0\0\0\!
 
 
 0\0
 
 I don't know what this error should tell me. Then I tried to use the
 ASCII format
 
 data.ascii - serialize(iris, NULL, ascii = TRUE)
 data.raw - rawToChar(data.ascii)
 dbSendQuery(db, DELETE FROM frames)
 dbSendQuery(db, paste(INSERT INTO frames VALUES(1, X', data.raw,
 '), sep = ))
 Error in sqliteExecStatement(conn, statement, ...) :
 RS-DBI driver: (error in statement: unrecognized token: X'A
 
 This also does not work. It seems the driver does not deal that nicely
 with the regular INSERT query for BLOB objects in SQLite. Then I used a
 simpler way:
 
 dbSendQuery(db, DELETE FROM frames)
 dbSendQuery(db, DROP TABLE frames)
 dbSendQuery(db, 'CREATE TABLE frames(simID INT, data TEXT DEFAULT
 NULL)')
 dbSendQuery(db, paste(INSERT INTO frames VALUES(1, ', data.raw, '),
 sep = ))
 data.bin2 - dbGetQuery(db, SELECT data FROM frames WHERE simID = 1)
 
 Nice, that worked. Now I want to unserialize the data:
 
 unserialize(data.bin2)
 Error in unserialize(data.bin2) : 'connection' must be a connection
 
 unserialize(data.bin2[1, 'data'])
 Error in unserialize(data.bin2[1, data]) :
 character vectors are no longer accepted by unserialize()
 
 I feel a little stuck here, but I am very sure, that converting
 data.frames to binary data and storing them to a database is not that
 unusual. So I hope somebody has already done this and could give me the
 missing piece.
 
 
 Best
 
 Simon
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r

Re: [R] Serialize data.frame to database

2013-07-16 Thread Simon Zehnder
Hi Rainer,

dbWriteTable is a nice function but in my case I need something that can 
actually save a dataframe in one row of a table. That is why I want to 
serialize my data.frame. 


Best

Simon


On Jul 16, 2013, at 3:05 PM, Rainer Schuermann rainer.schuerm...@gmx.net 
wrote:

 Maybe a simple
 
 dbWriteTable( db, frames, iris )
 
 does what you want?
 
 
 
 On Monday 15 July 2013 23:43:18 Simon Zehnder wrote:
 Dear R-Users,
 
 I need a very fast and reliable database solution so I try to serialize a 
 data.frame (to binary data) and to store this data to an SQLite database. 
 
 This is what I tried to do:
 
 library(RSQLite)
 con - dbDriver(SQLite)
 db - dbConnect(con, test)
 dbSendQuery(db, 'CREATE TABLE frames(simID INT, data BLOB)')
 data.bin - serialize(iris, NULL, ascii = FALSE)
 dbSendQuery(db, paste(INSERT INTO frames VALUES(1, X', data.bin, '), sep 
 = ))
 data.bin2 - dbGetQuery(db, SELECT DATA FROM frames WHERE simID = 1)
 data.bin2
  data
 1   58
 
 So, only the first entry of data.bin is saved to the database. I tried to 
 first convert the binary data to raw data:
 data.raw - rawToChar(data.bin)
 Error in rawToChar(data.bin) :
  embedded nul in string: 
 'X\n\0\0\0\002\0\003\0\001\0\002\003\0\0\0\003\023\0\0\0\005\0\0\0\016\0\0\0\x96@\024ff@\023\x99\x99\x99\x99\x99\x9a@\022\xcc\xcc\xcc\xcc\xcc\xcd@\022ff@\024\0\0\0\0\0\0@\025\x99\x99\x99\x99\x99\x9a@\022ff@\024\0\0\0\0\0\0@\021\x99\x99\x99\x99\x99\x9a@\023\x99\x99\x99\x99\x99\x9a@\025\x99\x99\x99\x99\x99\x9a@\02333@\02333@\02133@\02733@\026\xcc\xcc\xcc\xcc\xcc\xcd@\025\x99\x99\x99\x99\x99\x9a@\024ff@\026\xcc\xcc\xcc\xcc\xcc\xcd@\024ff@\025\x99\x99\x99\x99\x99\x9a@\024ff@\022ff@\024ff@\02333@\024\0\0\0\0\0\0@\024\0\0\0\0\0\0@\024\xcc\xcc\xcc\xcc\xcc\xcd@\024\xcc\xcc\xcc\xcc\xcc\xcd@\022\xcc\xcc\xcc\xcc\xcc\xcd@\02333@\025\x99\x99\x99\x99\x99\x9a@\024\xcc\xcc\xcc\xcc\xcc\xcd@\026\0\0\0\0\0\0@\023\x99\x99\x99\x99\x99\x9a@\024\0\0\0\0\0\0@\026\0\0\0\0\0\0@\023\x99\x99\x99\x99\x99\x9a@\021\x99\x99\x99\x99\x99\x9a@\024ff@\024\0\0\0\0\0\0@\022\0\0\0\0\0\0@\021\x99\x99\x99\x99\x99\x9a@\024\0\0\0\!
 0\!
 0\0
 
 I don't know what this error should tell me. Then I tried to use the ASCII 
 format
 
 data.ascii - serialize(iris, NULL, ascii = TRUE)
 data.raw - rawToChar(data.ascii)
 dbSendQuery(db, DELETE FROM frames)
 dbSendQuery(db, paste(INSERT INTO frames VALUES(1, X', data.raw, '), sep 
 = ))
 Error in sqliteExecStatement(conn, statement, ...) :
  RS-DBI driver: (error in statement: unrecognized token: X'A
 
 This also does not work. It seems the driver does not deal that nicely with 
 the regular INSERT query for BLOB objects in SQLite. Then I used a simpler 
 way:
 
 dbSendQuery(db, DELETE FROM frames)
 dbSendQuery(db, DROP TABLE frames)
 dbSendQuery(db, 'CREATE TABLE frames(simID INT, data TEXT DEFAULT NULL)')
 dbSendQuery(db, paste(INSERT INTO frames VALUES(1, ', data.raw, '), sep 
 = ))
 data.bin2 - dbGetQuery(db, SELECT data FROM frames WHERE simID = 1)
 
 Nice, that worked. Now I want to unserialize the data:
 
 unserialize(data.bin2)
 Error in unserialize(data.bin2) : 'connection' must be a connection
 
 unserialize(data.bin2[1, 'data'])
 Error in unserialize(data.bin2[1, data]) :
  character vectors are no longer accepted by unserialize()
 
 I feel a little stuck here, but I am very sure, that converting data.frames 
 to binary data and storing them to a database is not that unusual. So I hope 
 somebody has already done this and could give me the missing piece.
 
 
 Best
 
 Simon
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 - - - - -
 
 Der NSA keine Chance: e-mail verschluesseln!
 http://www.gpg4win.org/
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Setting Derived Class Slots

2013-07-16 Thread Simon Zehnder
Hi Steve,

it seems to me, that you want to pass an object inside of getDelayProfile. In 
this case you must use setReplaceMethod(getDelayProfile-,) in your 
definition inside the virtual and of course also in the specification of 
OP_Appt. R does not know, that your function should give back an object, that 
replaces the one on which it was called.

Otherwise making any modifications on the slot of the object inside the 
function will never arrive at the object outside of the function. They have two 
different environments and as soon as you modify an object inside of a function 
R builds a copy of it. See more about this by googling R pass-by-value.  

By the way: I wouldn't call this function getDelayProfile as getters and 
setters have usually differing a different meaning (encapsulation) in OO 
programming. Use rather something like fetchDelayProfile or 
initializeDelayProfile. It will probably make it easier for Users and 
Programmers to work with your program. 


Hope this helps

Best 

Simon

On Jul 16, 2013, at 3:36 PM, Steve Creamer stephen.crea...@nhs.net wrote:

 


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Fitting Mixture distributions

2013-07-16 Thread Simon Zehnder
Hi Tjun Kiat Teo,

you try to fit a Normal mixture to some data. The Normal mixture is very 
delicate when it comes to parameter search: If the variance gets closer and 
closer to zero, the log Likelihood becomes larger and larger for any values of 
the remaining parameters. Furthermore for the EM algorithm it is known, that it 
takes sometimes very long until convergence is reached. 

Try the following: 

Use as starting values for the component parameters:

start.par - mean(your.data, na.rm = TRUE) + sd(your.data, na.rm = TRUE) * 
runif(K)

For the weights just use either 1/K or the R cluster function with K clusters

Here K is the number of components. Further enlarge the maximum number of 
iterations. What you could also try is to randomize start parameters and run an 
SEM (Stochastic EM). In my opinion the better method is in this case a Bayesian 
method: MCMC.


Best

Simon


On Jul 16, 2013, at 10:59 PM, Tjun Kiat Teo teotj...@gmail.com wrote:

 I was trying to use the normixEM in mixtools and I got this error message.
 
 And I got this error message
 
 One of the variances is going to zero;  trying new starting values.
 Error in normalmixEM(as.matrix(temp[[gc]][, -(f + 1)])) : Too many tries!
 
 Are there any other packages for fitting mixture distributions  ?
 
 
 Tjun Kiat Teo
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Serialize data.frame to database

2013-07-15 Thread Simon Zehnder
Dear R-Users,

I need a very fast and reliable database solution so I try to serialize a 
data.frame (to binary data) and to store this data to an SQLite database. 

This is what I tried to do:

library(RSQLite)
con - dbDriver(SQLite)
db - dbConnect(con, test)
dbSendQuery(db, 'CREATE TABLE frames(simID INT, data BLOB)')
data.bin - serialize(iris, NULL, ascii = FALSE)
dbSendQuery(db, paste(INSERT INTO frames VALUES(1, X', data.bin, '), sep = 
))
data.bin2 - dbGetQuery(db, SELECT DATA FROM frames WHERE simID = 1)
data.bin2
  data
1   58

So, only the first entry of data.bin is saved to the database. I tried to first 
convert the binary data to raw data:
data.raw - rawToChar(data.bin)
Error in rawToChar(data.bin) :
  embedded nul in string: 
'X\n\0\0\0\002\0\003\0\001\0\002\003\0\0\0\003\023\0\0\0\005\0\0\0\016\0\0\0\x96@\024ff@\023\x99\x99\x99\x99\x99\x9a@\022\xcc\xcc\xcc\xcc\xcc\xcd@\022ff@\024\0\0\0\0\0\0@\025\x99\x99\x99\x99\x99\x9a@\022ff@\024\0\0\0\0\0\0@\021\x99\x99\x99\x99\x99\x9a@\023\x99\x99\x99\x99\x99\x9a@\025\x99\x99\x99\x99\x99\x9a@\02333@\02333@\02133@\02733@\026\xcc\xcc\xcc\xcc\xcc\xcd@\025\x99\x99\x99\x99\x99\x9a@\024ff@\026\xcc\xcc\xcc\xcc\xcc\xcd@\024ff@\025\x99\x99\x99\x99\x99\x9a@\024ff@\022ff@\024ff@\02333@\024\0\0\0\0\0\0@\024\0\0\0\0\0\0@\024\xcc\xcc\xcc\xcc\xcc\xcd@\024\xcc\xcc\xcc\xcc\xcc\xcd@\022\xcc\xcc\xcc\xcc\xcc\xcd@\02333@\025\x99\x99\x99\x99\x99\x9a@\024\xcc\xcc\xcc\xcc\xcc\xcd@\026\0\0\0\0\0\0@\023\x99\x99\x99\x99\x99\x9a@\024\0\0\0\0\0\0@\026\0\0\0\0\0\0@\023\x99\x99\x99\x99\x99\x9a@\021\x99\x99\x99\x99\x99\x9a@\024ff@\024\0\0\0\0\0\0@\022\0\0\0\0\0\0@\021\x99\x99\x99\x99\x99\x9a@\024\0\0\0\0\!
 0\0

I don't know what this error should tell me. Then I tried to use the ASCII 
format

data.ascii - serialize(iris, NULL, ascii = TRUE)
data.raw - rawToChar(data.ascii)
dbSendQuery(db, DELETE FROM frames)
dbSendQuery(db, paste(INSERT INTO frames VALUES(1, X', data.raw, '), sep = 
))
Error in sqliteExecStatement(conn, statement, ...) :
  RS-DBI driver: (error in statement: unrecognized token: X'A

This also does not work. It seems the driver does not deal that nicely with the 
regular INSERT query for BLOB objects in SQLite. Then I used a simpler way:

dbSendQuery(db, DELETE FROM frames)
dbSendQuery(db, DROP TABLE frames)
dbSendQuery(db, 'CREATE TABLE frames(simID INT, data TEXT DEFAULT NULL)')
dbSendQuery(db, paste(INSERT INTO frames VALUES(1, ', data.raw, '), sep = 
))
data.bin2 - dbGetQuery(db, SELECT data FROM frames WHERE simID = 1)

Nice, that worked. Now I want to unserialize the data:

unserialize(data.bin2)
Error in unserialize(data.bin2) : 'connection' must be a connection

unserialize(data.bin2[1, 'data'])
Error in unserialize(data.bin2[1, data]) :
  character vectors are no longer accepted by unserialize()

I feel a little stuck here, but I am very sure, that converting data.frames to 
binary data and storing them to a database is not that unusual. So I hope 
somebody has already done this and could give me the missing piece.


Best

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Set window title for plot on any OS

2013-07-14 Thread Simon Zehnder
Haha, good point!

So, the correct code is:

any(names(formals(getOption(device))) == title) {
dev.new(title = title)
}

Thanks for correcting my approach Rolf!


Best 

Simon

On Jul 14, 2013, at 1:02 AM, Rolf Turner rolf.tur...@xtra.co.nz wrote:

 
 
 I think you ought to take a look at
 
fortune(Lewis Carroll)
 
cheers,
 
Rolf Turner
 
 On 14/07/13 01:45, Simon Zehnder wrote:
 Hi Duncan,
 
 thank you very much for your advice! That makes it all work.
 
 I check in addition for a title argument in the device via
 
 if (any(names(getOption(device)) == TRUE)) {
  dev.new(title = title)
 }
 
 Thanks again!
 
 Simon
 
 On Jul 13, 2013, at 2:40 PM, Duncan Murdoch murdoch.dun...@gmail.com wrote:
 
 On 13-07-13 1:33 PM, Simon Zehnder wrote:
 Dear R-Users,
 
 I am writing a package using S4 classes. In the generic method plot I 
 want to set the title for the plotting window as I will have several 
 windows and window titles help the users to distinguish the graphics 
 without putting a title into the plot itself (this can be done by users 
 whenever they want)
 
 So I created a helper function .setDeviceTitle which I called after the 
 plot has been done:
 
 .setDeviceTitle - function(title = title, dev = dev.cur()) {
dev - names(dev)[1]

## check for OS ##
if (dev == windows) {
windows(title = title)
} else if (dev == X11) {
X11(title = title)
} else {
quartz(title = title)
}
 }
 
 The result is a new device with the title in addition to the old. Is it 
 possible to give a window a title after the plot has been done? If not: 
 Before I plot the device I cannot know what device it will be, so I 
 thought about a check via capabilities():
 
 if (any(names(capabilities()) == X11)) {
X11(title = title)
 }
 else if (any(names(capabilities)) == windows) {
windows(title = title)
 } else {
quartz(title = title)
 }
 
 I want to have a safe method, which works on each OS R can run. How would 
 you solve the problem?
 Use dev.new() rather than picking a particular device.  If all the possible 
 devices support the title argument, then
 
 dev.new(title=title)
 
 will be fine.  If you might need more customization (or want to protect 
 against a user who chooses a device that doesn't have title as an 
 argument), use getOption(device) to examine the call that will be used.
 
 Duncan Murdoch
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Set window title for plot on any OS

2013-07-13 Thread Simon Zehnder
Dear R-Users,

I am writing a package using S4 classes. In the generic method plot I want to 
set the title for the plotting window as I will have several windows and window 
titles help the users to distinguish the graphics without putting a title into 
the plot itself (this can be done by users whenever they want) 

So I created a helper function .setDeviceTitle which I called after the plot 
has been done:

.setDeviceTitle - function(title = title, dev = dev.cur()) {
dev - names(dev)[1]

## check for OS ##
if (dev == windows) {
windows(title = title)
} else if (dev == X11) {
X11(title = title)
} else {
quartz(title = title)
}
}

The result is a new device with the title in addition to the old. Is it 
possible to give a window a title after the plot has been done? If not: Before 
I plot the device I cannot know what device it will be, so I thought about a 
check via capabilities():

if (any(names(capabilities()) == X11)) {
X11(title = title)
}
else if (any(names(capabilities)) == windows) {
windows(title = title)
} else {
quartz(title = title)
}

I want to have a safe method, which works on each OS R can run. How would you 
solve the problem? 


Best 

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Set window title for plot on any OS

2013-07-13 Thread Simon Zehnder
Hi Duncan,

thank you very much for your advice! That makes it all work. 

I check in addition for a title argument in the device via

if (any(names(getOption(device)) == TRUE)) {
dev.new(title = title)
} 

Thanks again!

Simon

On Jul 13, 2013, at 2:40 PM, Duncan Murdoch murdoch.dun...@gmail.com wrote:

 On 13-07-13 1:33 PM, Simon Zehnder wrote:
 Dear R-Users,
 
 I am writing a package using S4 classes. In the generic method plot I want 
 to set the title for the plotting window as I will have several windows and 
 window titles help the users to distinguish the graphics without putting a 
 title into the plot itself (this can be done by users whenever they want)
 
 So I created a helper function .setDeviceTitle which I called after the plot 
 has been done:
 
 .setDeviceTitle - function(title = title, dev = dev.cur()) {
  dev - names(dev)[1]
  
  ## check for OS ##
  if (dev == windows) {
  windows(title = title)
  } else if (dev == X11) {
  X11(title = title)
  } else {
  quartz(title = title)
  }
 }
 
 The result is a new device with the title in addition to the old. Is it 
 possible to give a window a title after the plot has been done? If not: 
 Before I plot the device I cannot know what device it will be, so I thought 
 about a check via capabilities():
 
 if (any(names(capabilities()) == X11)) {
  X11(title = title)
 }
 else if (any(names(capabilities)) == windows) {
  windows(title = title)
 } else {
  quartz(title = title)
 }
 
 I want to have a safe method, which works on each OS R can run. How would 
 you solve the problem?
 
 Use dev.new() rather than picking a particular device.  If all the possible 
 devices support the title argument, then
 
 dev.new(title=title)
 
 will be fine.  If you might need more customization (or want to protect 
 against a user who chooses a device that doesn't have title as an 
 argument), use getOption(device) to examine the call that will be used.
 
 Duncan Murdoch
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Optim seems not to work properly in foreach

2013-06-06 Thread Simon Zehnder
Hi Ilai,

after you sent this message I tried your code as well and it worked. As a 
result, I reconsidered the code written by me and of course also found the 
error in my function simulateEKOP. So for other people assuming errors or ill 
behavior in base functions: Forget it: these functions are tested this often, 
that the probabilitty is almost 1, the error is in your own code. 

Thanks for the help

Simon


On Jun 3, 2013, at 10:53 PM, ilai ke...@math.montana.edu wrote:

 On Mon, Jun 3, 2013 at 11:37 AM, Simon Zehnder szehn...@uni-bonn.de
 
 ... [Some not minimal, self contained, reproducible code]...

 Data simulation and thecreation of startpar works fine, but the parameters in 
 res$par are always the start parameters. If I run the same commands directly 
 on the shell I get in res$par the optimized parameters - only inside the 
 foreach loop optim seems not to work. What could that be?
 
 Don't know, but but this makes me doubt it has anything to do with optim 
 being inside foreach:   
 
 fr - function(x) { 
  x1 - x[1] ; x2 - x[2] 
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2
  }
 grr - function(x) {
  x1 - x[1] ; x2 - x[2]
  c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1) , 200 * (x2 - x1 * x1))
  }
 library(doMC) 
 registerDoMC(2)
 RNGkind(L'Ecuyer)
 set.seed(54321)
 foreach(i = 1:2) %do% {
   ret - foreach(j = 1:2) %do%{
strtpar - c(-2,2)+rnorm(2) 
optim(strtpar, fr, grr, method = 
 L-BFGS-B,control=list(trace=TRUE))$par
   } 
   ret 
 } 
 
 Also, wouldn't you want to register 4 cores by default if nesting 2 loops of 
 2 ? (to comment on the wisdom of doing so in terms of overhead is beyond my 
 expertise)  
  
 HTH 
 
 
 
 Best
 
 Simon
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] 3.0.1 update and compiler package

2013-06-04 Thread Simon Zehnder
Hi Christoph,

do you install from sources? 

Best 

Simon

On Jun 4, 2013, at 10:41 AM, Christoph Knapp christoph.knap...@gmail.com 
wrote:

 Hi,
 reinstalling R did not help. It still will not update the same packages. I 
 attached the terminal output. There were no errors while I was installing R 
 again. I might have to remove all libraries and reinstall them all. Would you 
 agree or do you think I should try something else first which is not as 
 extreme.
 
 Thanks for your help
 
 Christoph
 
 On 03/06/13 20:32, Pascal Oettli wrote:
 My mistake,
 
 Regards,
 Pascal
 
 On 06/03/2013 04:37 PM, Uwe Ligges wrote:
 
 
 On 03.06.2013 07:19, Pascal Oettli wrote:
 Hi,
 
 How did you upgraded your version of R? From source or from a Linux
 package?
 
 Actually the new R installation is just broken. It simply has to be
 reinstalled carefully (watch for errors).
 
 Best,
 Uwe Ligges
 
 
 
 
 
 
 
 Regards,
 Pascal
 
 
 On 05/31/2013 11:33 AM, Christoph Knapp wrote:
 Hi,
 I recently updated to R 3.0.1. I'm running linux ubuntu 12.04. I
 realized that I have to update all the installed packages so I run 
 update.packages(checkBuilt=TRUE)
 as described here
 http://www.r-bloggers.com/r-3-0-0-is-released-whats-new-and-how-to-upgrade/
  
 
 
 . The first thing it did was telling me that it will not update several
 packages. When it was finished I used the warnings function to have a
 look at what did not work. See the list below.
 
  warnings()
 Warning messages:
 1: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘boot’ had non-zero exit status
 2: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘cluster’ had non-zero exit status
 3: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘foreign’ had non-zero exit status
 4: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘KernSmooth’ had non-zero exit status
 5: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘MASS’ had non-zero exit status
 6: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘Matrix’ had non-zero exit status
 7: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘nlme’ had non-zero exit status
 8: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘nnet’ had non-zero exit status
 9: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘rpart’ had non-zero exit status
 10: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘spatial’ had non-zero exit status
 11: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘survival’ had non-zero exit status
 12: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘class’ had non-zero exit status
 13: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘epiR’ had non-zero exit status
 14: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘gmodels’ had non-zero exit status
 15: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘gplots’ had non-zero exit status
 16: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘mgcv’ had non-zero exit status
 17: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘gregmisc’ had non-zero exit status
 
 I tried to reinstall them manually but this always failed because of a
 package dependency to the compiler package. Now, if I try to install
 the compiler package it tells me.
 
  install.packages(compiler)
 Installing package into
 ‘/home/christoph/R/x86_64-pc-linux-gnu-library/3.0’
 (as ‘lib’ is unspecified)
 Warning message:
 package ‘compiler’ is not available (for R version 3.0.1)
 
 The last line also came up all the time when the packages were updated
 
 Doing a bit of research does not deliver much only that the compiler
 package was included into R at version 2.13.0
 (http://dirk.eddelbuettel.com/blog/2011/04/12/).
 
 Most of those packages which do not work any more are pretty important
 for some of my scripts and I would not even know what packages replace
 the packages above.
 
 Would anyone know how to fix this?
 
 Regards
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 
 

Re: [R] 3.0.1 update and compiler package

2013-06-04 Thread Simon Zehnder
Hi Christoph,

so what you could try to make sure it can run easily on your system is to 
install it from sources:

First check if you have at least gcc version 4.7.2:
gcc --version

Open a terminal and type 
wget ftp://ftp.stat.math.ethz.ch/Software/R/R-patched.tar.gz
tar -zxvf R-patched.tar.gz 
mkdir build  cd build
../R-patched/configure --prefix=/usr/opt/R --with-blas='-framework libVec' 
--with-lapack

Check the output at the end of the configuration process. If everything looks 
fine start the build:
make -j4

This should run without errors. Then check if all works fine:
make check 

If this is the case you can install it (up to here no installation have been 
done: so the folder /usr/opt/R does not exist  yet):
sudo make install

Start R via the shell:
/usr/opt/R/bin/R 

Install your packages and see what happens. If all is fine, add /usr/opt/R to 
your PATH variable, if not just delete the folder /usr/opt/R. It all worked 
perfectly on my scientific Linux system. 


Best

Simon


On Jun 4, 2013, at 11:50 AM, Christoph Knapp christoph.knap...@gmail.com 
wrote:

 Usually sudo apt-get install ... Sometimes synaptic package manager. I'm 
 using the 
 
 http://cran.ma.imperial.ac.uk/bin/linux/ubuntu/precise/
 http://cran.stat.ucla.edu/bin/linux/ubuntu/precise
 
 repositories. Not sure what the differences are. I had some problems before 
 when I installed packages from within Rstudio. It would install them in a 
 different directory than apt-get would and it would tell me that those 
 packages were not installed when I used the terminal to start R. Drove me 
 crazy for a while because I could use them in Rstudio so I was sure they were 
 installed. Do you think this could also cause this? 
 
 Christoph
 
 
 On 04/06/13 21:13, Simon Zehnder wrote:
 Hi Christoph,
 
 do you install from sources? 
 
 Best 
 
 Simon
 
 On Jun 4, 2013, at 10:41 AM, Christoph Knapp christoph.knap...@gmail.com 
 wrote:
 
 Hi,
 reinstalling R did not help. It still will not update the same packages. I 
 attached the terminal output. There were no errors while I was installing R 
 again. I might have to remove all libraries and reinstall them all. Would 
 you agree or do you think I should try something else first which is not as 
 extreme.
 
 Thanks for your help
 
 Christoph
 
 On 03/06/13 20:32, Pascal Oettli wrote:
 My mistake,
 
 Regards,
 Pascal
 
 On 06/03/2013 04:37 PM, Uwe Ligges wrote:
 
 On 03.06.2013 07:19, Pascal Oettli wrote:
 Hi,
 
 How did you upgraded your version of R? From source or from a Linux
 package?
 Actually the new R installation is just broken. It simply has to be
 reinstalled carefully (watch for errors).
 
 Best,
 Uwe Ligges
 
 
 
 
 
 
 Regards,
 Pascal
 
 
 On 05/31/2013 11:33 AM, Christoph Knapp wrote:
 Hi,
 I recently updated to R 3.0.1. I'm running linux ubuntu 12.04. I
 realized that I have to update all the installed packages so I run 
 update.packages(checkBuilt=TRUE)
 as described here
 http://www.r-bloggers.com/r-3-0-0-is-released-whats-new-and-how-to-upgrade/
  
 
 
 . The first thing it did was telling me that it will not update several
 packages. When it was finished I used the warnings function to have a
 look at what did not work. See the list below.
 
 warnings()
 Warning messages:
 1: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘boot’ had non-zero exit status
 2: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘cluster’ had non-zero exit status
 3: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘foreign’ had non-zero exit status
 4: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘KernSmooth’ had non-zero exit status
 5: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘MASS’ had non-zero exit status
 6: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘Matrix’ had non-zero exit status
 7: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘nlme’ had non-zero exit status
 8: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘nnet’ had non-zero exit status
 9: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘rpart’ had non-zero exit status
 10: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘spatial’ had non-zero exit status
 11: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘survival’ had non-zero exit status
 12: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘class’ had non-zero exit status
 13: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘epiR’ had non-zero exit status
 14: In install.packages(update[instlib == l, Package], l, ... :
 installation of package ‘gmodels’ had non-zero exit status
 15

[R] Optim seems not to work properly in foreach

2013-06-03 Thread Simon Zehnder
Hi guys,

I am working now for several years in R and I would say I manage things pretty 
easily, but with the foreach loop I have my problems. I call for a simulation a 
double foreach loop and this works fine. Inside the second loop (which I plan 
to parallelize later on) I call Rs optim-function:

   
simulateML - function(sim.it = 250, input.path,
output.path, T = 390, ncore = 2) {

## load packages ##
library(mmstruct)
library(doMC)
library(foreach)

## read input parameters ##
input - read.csv(input.path, header = FALSE)

## create container to store results ##
results - data.frame(NA, nrow = NROW(input) * sim.it, ncol = 12)

## initialize the parallel backend ##
registerDoMC(ncore)

result.list - foreach(i = 1:2) %do% {
input - as.matrix(input)
input.row - input[i, ]
list - foreach(i = 1:2, .combine = rbind, .packages = 
c(mmstruct, stats)) %do% {
data - simulateEKOP(size = input.row[1], alpha = 
input.row[2],   
epsilon = input.row[1], delta = 
input.row[4],
mu = input.row[5], T = T)
data - data[,4]
## create start values ##
tmp - mean(data)/T
startpar - c(0, tmp * 0.75/2, tmp * 0.25/2)

## set options for optimization ##
optim_fn- computeKokotLik
optim_method- L-BFGS-B
optim_lower - c(-1e+7, 0, 0)
optim_upper - rep(1e+7, 3)
optim_fnscale   - -1
optim_maxit - 200
optim_ctrl  - list(fnscale = optim_fnscale, maxit 
= optim_maxit)

## start optimization ##
res - optim(par = startpar, fn = optim_fn, data = 
data, T = T,
methodLik = approx, method = 
optim_method,
lower = optim_lower, upper = 
optim_upper,
control = optim_ctrl, hessian = TRUE)
res$par

}
print(list)
}

}

Data simulation and thecreation of startpar works fine, but the parameters in 
res$par are always the start parameters. If I run the same commands directly on 
the shell I get in res$par the optimized parameters - only inside the foreach 
loop optim seems not to work. What could that be? 


Best 

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Count function calls

2013-03-07 Thread Simon Zehnder
Dear Giovanni,

apologize for this late reply! I was testing and reading a lot of stuff. I 
tried your suggestions and the problem of singularity in the regressor cross 
product vanishes when using the Group Mean function 'pgm' instead of 'pvcm'. 
Nevertheless, I found the collinearity in the regressors and 'pvcm' ran 
further, until I got another error. This time a little different but still 
directed towards the same area:

'Error in solve.default(crossprod(X[[i]])[!coefna[i, ], !coefna[i, ]])) :
system is computationally singular: reciprocal condition number = 
1.86131e-17'

Following the description of R's function 'rcond', the reciprocal condition 
number measures how close a matrix is to be rank deficient. So in this case one 
regressor crossproduct seems to be sufficiently close to singularity, such that 
inversion becomes impossible. 

I attached you a simple subsample. If you let run the following commands on it:

require(plm)
data - read.csv(path, sep = ,, header = TRUE)
pdata - plm.data(data, index = c(gvkey, fyear))
debug(plm:::pvcm)
model.pvcm - pvcm(LOGDXSGA~LOGDSALE + DECRLOGDSALE + DECRLOGDSALEWDAVG, data = 
pdata, model = random)
.. go until
Browse[2] 
debug: ml - split(data cond)

Then type:
A - as.matrix(ml$'25062'[, 2:4])
XX - t(A) %*% A
qr(XX)$rank
[1] 2Aha, the rank is only 2! 
rcond(XX) 
[1] 9.339817e-16Alright, the matrix XX is pretty near to singularity, even 
it is not really singular!

Let us have a look at it:
A[, 2]/A[, 3]
.  A[, 2] is almost a multiple of A[, 3], if 
the latter is multiplied by -6.231979 or -6.231331 . at least it suffices 
to let the rank shrink to 2. 

The thing is, the same regression goes through in the Stata package. As I would 
say, I am an R fanatist, I would like to know, why it runs in Stata and not in 
R.and that one can let it run in R as well. Different number 
formats/precision? Maybe it suffices in such a case to enforce full rank of the 
crossproduct by enforcing positive definitess, for instance via the function 
'nearPD' in the R package 'Matrix':

Try 
library(Matrix)
qr(nearPD(XX)$mat)$rank
[1] 3 Yey!. :)

Tell me your opinion Giovanni! Give me a little help here please! 


Best 

Simon
 

P.S. Thank your for this fantastic package making panel data estimation 
possible in the R environment! 




On Feb 27, 2013, at 1:40 PM, Millo Giovanni giovanni_mi...@generali.com wrote:

 Hello.
 Another thing you may want to do depends on whether you are using
 model=within (the default) or model=random.
 
 In the first case, pvcm() estimates separate regressions, so you just
 need to loop lm() on individual indices to spot where it fails.
 
 In the second case, what you may want to do is try a similar estimator:
 pmg(..., model=mg), which is an unweighted version of Swamy's
 estimator in pvcm() and a simpler function to read and modify, possibly
 using the global assignment operator '-' as already suggested by
 William to output diagnostics.
 
 Best,
 Giovanni
 
 -Original Message-
 From: r-help-boun...@r-project.org
 [mailto:r-help-boun...@r-project.org] On Behalf
 Of Simon Zehnder
 Sent: Tuesday, February 26, 2013 2:53 AM
 To: r-help@r-project.org help
 Subject: [R] Count function calls
 
 Dear R-users,
 
 I have the following problem: I am running the function 'pvcm' from
 the 'plm' Panel Data
 package. Inside this function 'solve' is called and gives for a
 certain individual data series
 an exception because of singularity. I would like to know which
 individual data series
 causes this error. I tried to debug it, but this is truly painful, as
 solve is called inside of
 'lapply' and there are over 5,000 individual data series in the panel.
 
 Now, what I would like to do is to count the calls to 'solve' inside
 the function, so I can
 see, where the function throws the exception. I tried to use 'trace'
 with a count variable,
 but I have no clue how to define a global variable to be used by trace
 and updated at
 every call.is there another approach?
 
 
 Best
 
 Simon
 
 __
 
  
 Ai sensi del D.Lgs. 196/2003 si precisa che le informazioni contenute in 
 questo messaggio sono riservate ed a uso esclusivo del destinatario. Qualora 
 il messaggio in parola Le fosse pervenuto per errore, La invitiamo ad 
 eliminarlo senza copiarlo e a non inoltrarlo a terzi, dandocene gentilmente 
 comunicazione. Grazie.
 
 Pursuant to Legislative Decree No. 196/2003, you are hereby informed that 
 this message contains confidential information intended only for the use of 
 the addressee. If you are not the addressee, and have received this message 
 by mistake, please delete it and immediately notify us. You may not copy or 
 disseminate this message to anyone. Thank you.

__
R-help@r-project.org mailing list
https://stat.ethz.ch

[R] S4-classes: Assignment of values to slots by reference

2013-03-01 Thread Simon Zehnder
Dear R-users,

I am working on a project that uses S4-classes. At some point I encountered the 
problem - well known to R - that I have to pass 3 different objects to a 
function, that should modify several slots of them and of course there is no 
passing by reference in R. 

Then I read this thread by Steve Lianoglou: 
https://stat.ethz.ch/pipermail/r-help/2010-August/250468.html, which offers 
from viewpoint of OOP an elegant solution. But I do not use an Object Method, 
but a regular function, pass in several objects and modify them.

Here is an example with 1 Object using Steve's approach: 

setClass(example, 
representation(
par = matrix,
.cache = environment)
)
setMethod(initialize, example, function(.Object, ..., .cache = new.env()) {

callNextMethod(.Object, .cache = .cache, ...)
}
)

ex - new(example, par = matrix())

fmodify - function(O1) {
pm - mean(rnorm(100))
pm - matrix(pm)
O1@.cache$par - pm
} 

fmodify(ex)


So far so good. Everything worked nicely. But of course the value computed in 
the function 'fmodify' is now in the slot ex@.cache$par:

ex@.cache$par

and not in ex@par. Is there a possibility to get it into ex@par?

Best

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Count function calls

2013-02-27 Thread Simon Zehnder
Dear WIlliam,

thank you very much for this very useful information! I learned something new! 

Enjoy your day!

Simon

On Feb 27, 2013, at 1:40 AM, William Dunlap wdun...@tibco.com wrote:

 This is where - is helpful:
 N - 0 ; trace(solve, quote(N - N + 1), print=FALSE)
   Tracing function solve in package base
   [1] solve
 lapply(3:0, function(i)solve(diag(i,3), 1:3))
   Error in solve.default(diag(i, 3), 1:3) : 
 Lapack routine dgesv: system is exactly singular: U[1,1] = 0
 N
   [1] 4
 
 You can also set
   options(error=recover)
 to look at the state of things when the error occurs.
 
 Bill Dunlap
 Spotfire, TIBCO Software
 wdunlap tibco.com
 
 
 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
 Behalf
 Of Simon Zehnder
 Sent: Tuesday, February 26, 2013 2:53 AM
 To: r-help@r-project.org help
 Subject: [R] Count function calls
 
 Dear R-users,
 
 I have the following problem: I am running the function 'pvcm' from the 
 'plm' Panel Data
 package. Inside this function 'solve' is called and gives for a certain 
 individual data series
 an exception because of singularity. I would like to know which individual 
 data series
 causes this error. I tried to debug it, but this is truly painful, as solve 
 is called inside of
 'lapply' and there are over 5,000 individual data series in the panel.
 
 Now, what I would like to do is to count the calls to 'solve' inside the 
 function, so I can
 see, where the function throws the exception. I tried to use 'trace' with a 
 count variable,
 but I have no clue how to define a global variable to be used by trace and 
 updated at
 every call.is there another approach?
 
 
 Best
 
 Simon
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Count function calls

2013-02-26 Thread Simon Zehnder
Dear R-users,

I have the following problem: I am running the function 'pvcm' from the 'plm' 
Panel Data package. Inside this function 'solve' is called and gives for a 
certain individual data series an exception because of singularity. I would 
like to know which individual data series causes this error. I tried to debug 
it, but this is truly painful, as solve is called inside of 'lapply' and there 
are over 5,000 individual data series in the panel. 

Now, what I would like to do is to count the calls to 'solve' inside the 
function, so I can see, where the function throws the exception. I tried to use 
'trace' with a count variable, but I have no clue how to define a global 
variable to be used by trace and updated at every call.is there another 
approach? 


Best

Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Armadillo error in R extension

2013-02-02 Thread Simon Zehnder
Dear Michael,

during the last days of programming I came across RcppArmadillo (Rcpp) and for 
my kind of work this is indeed pretty interesting. So, I will take a closer 
look today what it is about and how it works. 

As a C++ programmer though, I am still interested why compilation of my package 
fails at this point of code. I hope I get a little more insight on the 
Rcpp-devel list suggested by Dirk! 

Thanks for your comment!

Simon


On Feb 2, 2013, at 12:14 AM, Michael Weylandt michael.weyla...@gmail.com 
wrote:

 Look up the terribly wonderful RcppArmadillo package. 
 
 MW
 
 On Feb 1, 2013, at 8:38 PM, Simon Zehnder szehn...@uni-bonn.de wrote:
 
 Is there anyway with some experience in using armadillo in R C++ extensions?
 
 My problem is the following:
 
 I programmed a function in a header looking like
 
 #include armadillo
 
 inline arma::vec foo(input) {
 
   ... do something
 
   return an arma::vec object 
 }
 
 compiling this via R CMD INSTALL packagename (PKG_CXXFLAGS = 
 -I/folder/of/armadillo and armadillo_bits in my package)
 
 I get the following error 
 
 error: expected initializer before 'foo' 
 
 The exact line and word is where inline arma::vec ends.
 
 Does anyone know what this error could be? 
 
 
 Best
 
 Simon
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Armadillo error in R extension

2013-02-01 Thread Simon Zehnder
Is there anyway with some experience in using armadillo in R C++ extensions?

My problem is the following:

I programmed a function in a header looking like

#include armadillo

inline arma::vec foo(input) {

... do something

return an arma::vec object 
}

compiling this via R CMD INSTALL packagename (PKG_CXXFLAGS = 
-I/folder/of/armadillo and armadillo_bits in my package)

I get the following error 

error: expected initializer before 'foo' 

The exact line and word is where inline arma::vec ends.

Does anyone know what this error could be? 


Best

Simon
 
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Modify objects in function

2013-01-31 Thread Simon Zehnder
Dear R community,

I do know, that an R function is constructing a copy of any object passed as 
argument into a function. I program on a larger S4 project for a package, and I 
arrived at a point where I have to think a little harder on implementation 
style (especially to spare users complex object handling). 

I have a function foo(), taking as input arguments two S4 objects of different 
class type

foo - function(o1, o2) {
o1@att1 - producesomething()
o2@att2 - producesomethingelse()
 
}

Of course, this functions does not change the objects in the global 
environment. Now I have two choices 

1. Change the objects and return a list with both objects:

foo - function(o1, o2) {
o1@att1 - producesomething()
o2@att2 - producesomethingelse()

l - list(O1 = o1, O2 = o2)
return(l)
}

This is cumbersome for users, as they have then to assign the objects inside 
the returned list to the symbols used in the global environment. But it is the 
way intended by R.

2. Change the objects of the global environment inside the function:

foo - function(o1, o2) {
o1@att1 - producesomething()
o2@att2 - producesomethingelse()

assign(o1, o1, .GlobalEnv)
assign(o2, o2, .GlobalEnv)
}

This is for users very practical, as the symbols used before refer now to 
objects with changed attributes and can be used immediately. BUT: It is not 
intended to be done like this.

I spared any solution using one function for every single object type, as this 
becomes even more cumbersome than choice 1 above, in my opinion.

What is your opinion on the trade-off between user-friendly handling of objects 
and R-intended programming style? Do I miss another choice, combining both?


Best 
Simon
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Modify objects in function

2013-01-31 Thread Simon Zehnder
Dear Barry,

thank you very much for this information. This looks pretty interesting! 

Best

Simon
On Jan 31, 2013, at 4:09 PM, Barry Rowlingson b.rowling...@lancaster.ac.uk 
wrote:

 it lets you do:
 
 (a~b~c) = foo()
 
 Mistook. should be:
 
 (a~b~c) %=% foo()
 
 because it defines the %=% operator.
 
 Barry

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Modify objects in function

2013-01-31 Thread Simon Zehnder
Hi Barry,

this actually a good idea, to put them together! Probably even creating an 
object containing both of them. Haven't thought about it before. 

Best
Simon
On Jan 31, 2013, at 3:49 PM, Barry Rowlingson b.rowling...@lancaster.ac.uk 
wrote:

 On Thu, Jan 31, 2013 at 11:52 AM, Simon Zehnder szehn...@uni-bonn.de wrote:
 Dear R community,
 
 I do know, that an R function is constructing a copy of any object passed as 
 argument into a function. I program on a larger S4 project for a package, 
 and I arrived at a point where I have to think a little harder on 
 implementation style (especially to spare users complex object handling).
 
 I have a function foo(), taking as input arguments two S4 objects of 
 different class type
 
 foo - function(o1, o2) {
o1@att1 - producesomething()
o2@att2 - producesomethingelse()
 
 }
 
 Of course, this functions does not change the objects in the global 
 environment. Now I have two choices
 
 1. Change the objects and return a list with both objects:
 
foo - function(o1, o2) {
o1@att1 - producesomething()
o2@att2 - producesomethingelse()
 
l - list(O1 = o1, O2 = o2)
return(l)
}
 
 This is cumbersome for users, as they have then to assign the objects inside 
 the returned list to the symbols used in the global environment. But it is 
 the way intended by R.
 
 By cumbersome you mean the following (anti-?) pattern has to be used,
 for example in a loop:
 
 o1 = something()
 o2 = somethingelse()
 o12 = foo(o1,o2)
 o1 = o12$o1
 o2 = o12$o2
 
 so that at the end of it you have a modified o1 and o2?
 
 I suggest that if you have a function that modifies more than one of
 its arguments, then there is a strong case for making those arguments
 a single argument. So that you'd do:
 
 foo = function(o12){
o12$o1@att=bar()
o12$o2@att=baz()
return(o12)
  }
 
 o12 = list(something(), somethingelse())
 o12 = foo(o12)
 
 and there you are, a modified o12 object. no packing/unpacking needed.
 
 I suspect that either your o1 and o2 are so closely related that you
 should pack them in a list (or ideally, an object) or differently
 related such that you should treat them separately and not tweak both
 of them within the same function. That approach is binding behaviour
 very tightly to two structures, and probably won't give you very
 debuggable modular code.
 
 
 B

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Irreproducible exception in R extension

2012-12-28 Thread Simon Zehnder
Dear R-Users,

I am having some trouble running an R extension on our cluster (linux). I call 
C++ code in which I use a) the Scythe Statistical Library and b) OpenMP. Most 
of the jobs run without a problem, but some arbitrary jobs throw an exception 
of the kind printed below while running in a parallel loop. The behavior is to 
me not reproducible, although at every run of 729 Jobs it at least happens 
once. I already tried to enforce the error by using the Intel Compiler instead 
of the GCC 4.7.2 and using the Options '-checkinit' for possibly uninitialized 
values and '-ftrapuv' for initialization of stack variables to unusual values 
to aid error detection. But I cannot enforce the exception, which makes it 
impossible to debug. I already googled this error,  but there only 2-3 small 
forum discussions about it, and so I hope I find some suggestions here. 

1. What do you think this error comes from? 
2. What do you advice me to do to detect the error?

I am excited for your answers.

Simon




*** glibc detected *** /home/simo_109//opt/R/lib64/R/bin/exec/R: double free or 
corruption (!prev): 0x01355610 ***
=== Backtrace: =
/lib64/libc.so.6[0x3845a75916]
/lib64/libc.so.6[0x3845a78443]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe9DataBlockIdE10deallocateEv+0x7b)[0x7f4645d5d6db]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe9DataBlockIdED1Ev+0x47)[0x7f4645d5d5bb]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe18DataBlockReferenceIdE17withdrawReferenceEv+0xaa)[0x7f4645d83442]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe18DataBlockReferenceIdED1Ev+0x94)[0x7f4645d82cc8]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe18DataBlockReferenceIdED2Ev+0x47)[0x7f4645d82d9f]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe6MatrixIdLNS_12matrix_orderE0ELNS_12matrix_styleE1EED1Ev+0x5d)[0x7f4645d8b26d]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_Z11MCgmmS_implIN6scythe8mersenneEEvRNS0_3rngIT_EERP7SEXPRECS8_jjRNS0_6MatrixIdLNS0_12matrix_orderE0ELNS0_12matrix_styleE0EEERKSC_jS8_S8_+0x478f)[0x7f4645d6d6af]
/opt/intel/Compiler/12.1/6.361/rwthlnk/compiler/lib/intel64/libiomp5.so(__kmp_invoke_microtask+0x93)[0x7f4645aabee3]
=== Memory map: 
0040-00401000 r-xp  00:21 32550383   
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/bin/exec/R
0060-00601000 rw-p  00:21 32550383   
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/bin/exec/R
00601000-029cc000 rw-p  00:00 0  [heap]
369ac0-369ac3a000 r-xp  08:01 263771 
/lib64/libreadline.so.6.0
369ac3a000-369ae3a000 ---p 0003a000 08:01 263771 
/lib64/libreadline.so.6.0
369ae3a000-369ae42000 rw-p 0003a000 08:01 263771 
/lib64/libreadline.so.6.0
369ae42000-369ae43000 rw-p  00:00 0 
369bc0-369bc16000 r-xp  08:01 263757 
/lib64/libgcc_s-4.4.6-20120305.so.1
369bc16000-369be15000 ---p 00016000 08:01 263757 
/lib64/libgcc_s-4.4.6-20120305.so.1
369be15000-369be16000 rw-p 00015000 08:01 263757 
/lib64/libgcc_s-4.4.6-20120305.so.1
369c00-369c0e8000 r-xp  08:01 656025 
/usr/lib64/libstdc++.so.6.0.13
369c0e8000-369c2e8000 ---p 000e8000 08:01 656025 
/usr/lib64/libstdc++.so.6.0.13
369c2e8000-369c2ef000 r--p 000e8000 08:01 656025 
/usr/lib64/libstdc++.so.6.0.13
369c2ef000-369c2f1000 rw-p 000ef000 08:01 656025 
/usr/lib64/libstdc++.so.6.0.13
369c2f1000-369c306000 rw-p  00:00 0 
369d00-369d00d000 r-xp  08:01 671072 
/usr/lib64/libgomp.so.1.0.0
369d00d000-369d20c000 ---p d000 08:01 671072 
/usr/lib64/libgomp.so.1.0.0
369d20c000-369d20d000 rw-p c000 08:01 671072 
/usr/lib64/libgomp.so.1.0.0
369fa0-369faf r-xp  08:01 660772 
/usr/lib64/libgfortran.so.3.0.0
369faf-369fcef000 ---p 000f 08:01 660772 
/usr/lib64/libgfortran.so.3.0.0
369fcef000-369fcf1000 rw-p 000ef000 08:01 660772 
/usr/lib64/libgfortran.so.3.0.0
369fcf1000-369fcf2000 rw-p  00:00 0 
36a360-36a373f000 r-xp  08:01 675210 
/usr/lib64/libicuuc.so.42.1
36a373f000-36a393f000 ---p 0013f000 08:01 675210 
/usr/lib64/libicuuc.so.42.1
36a393f000-36a395 rw-p 0013f000 08:01 675210 
/usr/lib64/libicuuc.so.42.1
36a395-36a3952000 rw-p  00:00 0 
36a3a0-36a3b88000 r-xp 

[R] Irreproducible exception in R extension

2012-12-28 Thread Simon Zehnder
Dear R-Users,

I am having some trouble running an R extension on our cluster (linux). I call 
C++ code in which I use a) the Scythe Statistical Library and b) OpenMP. Most 
of the jobs run without a problem, but some arbitrary jobs throw an exception 
of the kind printed below while running in a parallel loop. The behavior is to 
me not reproducible, although at every run of 729 Jobs it at least happens 
once. I already tried to enforce the error by using the Intel Compiler instead 
of the GCC 4.7.2 and using the Options '-checkinit' for possibly uninitialized 
values and '-ftrapuv' for initialization of stack variables to unusual values 
to aid error detection. But I cannot enforce the exception, which makes it 
impossible to debug. I already googled this error,  but there only 2-3 small 
forum discussions about it, and so I hope I find some suggestions here. 

1. What do you think this error comes from? 
2. What do you advice me to do to detect the error?

I am excited for your answers.

Simon




*** glibc detected *** /home/simo_109//opt/R/lib64/R/bin/exec/R: double free or 
corruption (!prev): 0x01355610 ***
=== Backtrace: =
/lib64/libc.so.6[0x3845a75916]
/lib64/libc.so.6[0x3845a78443]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe9DataBlockIdE10deallocateEv+0x7b)[0x7f4645d5d6db]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe9DataBlockIdED1Ev+0x47)[0x7f4645d5d5bb]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe18DataBlockReferenceIdE17withdrawReferenceEv+0xaa)[0x7f4645d83442]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe18DataBlockReferenceIdED1Ev+0x94)[0x7f4645d82cc8]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe18DataBlockReferenceIdED2Ev+0x47)[0x7f4645d82d9f]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_ZN6scythe6MatrixIdLNS_12matrix_orderE0ELNS_12matrix_styleE1EED1Ev+0x5d)[0x7f4645d8b26d]
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkg/libs/MCpkg.so(_Z11MCgmmS_implIN6scythe8mersenneEEvRNS0_3rngIT_EERP7SEXPRECS8_jjRNS0_6MatrixIdLNS0_12matrix_orderE0ELNS0_12matrix_styleE0EEERKSC_jS8_S8_+0x478f)[0x7f4645d6d6af]
/opt/intel/Compiler/12.1/6.361/rwthlnk/compiler/lib/intel64/libiomp5.so(__kmp_invoke_microtask+0x93)[0x7f4645aabee3]
=== Memory map: 
0040-00401000 r-xp  00:21 32550383   
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/bin/exec/R
0060-00601000 rw-p  00:21 32550383   
/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/bin/exec/R
00601000-029cc000 rw-p  00:00 0  [heap]
369ac0-369ac3a000 r-xp  08:01 263771 
/lib64/libreadline.so.6.0
369ac3a000-369ae3a000 ---p 0003a000 08:01 263771 
/lib64/libreadline.so.6.0
369ae3a000-369ae42000 rw-p 0003a000 08:01 263771 
/lib64/libreadline.so.6.0
369ae42000-369ae43000 rw-p  00:00 0 
369bc0-369bc16000 r-xp  08:01 263757 
/lib64/libgcc_s-4.4.6-20120305.so.1
369bc16000-369be15000 ---p 00016000 08:01 263757 
/lib64/libgcc_s-4.4.6-20120305.so.1
369be15000-369be16000 rw-p 00015000 08:01 263757 
/lib64/libgcc_s-4.4.6-20120305.so.1
369c00-369c0e8000 r-xp  08:01 656025 
/usr/lib64/libstdc++.so.6.0.13
369c0e8000-369c2e8000 ---p 000e8000 08:01 656025 
/usr/lib64/libstdc++.so.6.0.13
369c2e8000-369c2ef000 r--p 000e8000 08:01 656025 
/usr/lib64/libstdc++.so.6.0.13
369c2ef000-369c2f1000 rw-p 000ef000 08:01 656025 
/usr/lib64/libstdc++.so.6.0.13
369c2f1000-369c306000 rw-p  00:00 0 
369d00-369d00d000 r-xp  08:01 671072 
/usr/lib64/libgomp.so.1.0.0
369d00d000-369d20c000 ---p d000 08:01 671072 
/usr/lib64/libgomp.so.1.0.0
369d20c000-369d20d000 rw-p c000 08:01 671072 
/usr/lib64/libgomp.so.1.0.0
369fa0-369faf r-xp  08:01 660772 
/usr/lib64/libgfortran.so.3.0.0
369faf-369fcef000 ---p 000f 08:01 660772 
/usr/lib64/libgfortran.so.3.0.0
369fcef000-369fcf1000 rw-p 000ef000 08:01 660772 
/usr/lib64/libgfortran.so.3.0.0
369fcf1000-369fcf2000 rw-p  00:00 0 
36a360-36a373f000 r-xp  08:01 675210 
/usr/lib64/libicuuc.so.42.1
36a373f000-36a393f000 ---p 0013f000 08:01 675210 
/usr/lib64/libicuuc.so.42.1
36a393f000-36a395 rw-p 0013f000 08:01 675210 
/usr/lib64/libicuuc.so.42.1
36a395-36a3952000 rw-p  00:00 0 
36a3a0-36a3b88000 r-xp 

Re: [R] Generating an autocorrelated binary variable

2012-10-02 Thread Simon Zehnder
Thank you very much for your answer Rolf. It helped. 

I try to simulate a trade indicator model from market microstructure, where the 
1 or -1 indicate a buyer or seller initiated trade respectively. 
I use a Gaussian copula for simulation, so I can put in some correlation if I 
want to. So I generate my multivariate Gaussian sample and transform it to a 
uniform sample so I can use whatever marginal distribution I want to. 
Then I generate the autocorrelated binary by beginning with a transformation on 
the first observation and then by using the uniform sample to compare it to 
0.6. Afterwards I have a 0.6 autocorrelated binary variable:

sampleCop - function(n = 1000, rho = 0.6) {

require(splus2R)
mvrs - rmvnorm(n + 1, mean = rep(0, 3), cov = diag(3))
pmvrs - pnorm(mvrs, 0, 1)
var1 - matrix(0, nrow = n + 1, ncol = 1)
var1[1] - qbinom(pmvrs[1, 1], 1, 0.5)
if(var1[1] == 0) var1[1] - -1
for(i in  1:(nrow(pmvrs) - 1)) {
if(pmvrs[i + 1, 1] = rho) var1[i + 1] - var1[i]
else var1[i + 1] - var1[i] * (-1)
}
sample - matrix(0, nrow = n, ncol = 4)
sample[, 2] - var1[1:nrow(var1) - 1]
sample[, 3] - var1[2:nrow(var1)]
sample[, 4] - qnorm(pmvrs[1:nrow(var1) - 1, 2], 0, 1, 1, 0) + 
qnorm(pmvrs[1:nrow(var1) - 1, 3], 0, 1, 1, 0) - qnorm(pmvrs[2:nrow(var1), 3], 
0, 1, 1, 0)
sample[, 1] - sample[,2] * (0.015 + 0.03) - sample[,3] * (0.015 + 
0.2*0.03) + sample[,4]

sample

}

It is interesting though, that if I run my GMM-model on it it needs very high 
numbers of observations to converge to the real parameters. Using the 
gmm-function

gmmfun - function(theta, data) {

res = data[,1] - (theta[1] + theta[2]) * data[,2] + (theta[1] + 
theta[3] * theta[2]) * data[,3]

moments = matrix(0, nrow=nrow(data), ncol=3)
moments[,1] = data[,2] * data[,3] - data[,2]^2 * theta[3]
moments[,2] = res * data[,2]
moments[,3] = res * data[,3]

return(moments)
}

And the objective function 

obj_fun - function(theta, data) {
moments - gmmfun(theta, data)
moments_sum - colSums(moments)
val - moments_sum %*% diag(3) %*% moments_sum * 1.0 / nrow(data)
val
}

Running then the command 

optim(fn = obj_fun, par = c(0.05, 0.05, .1), method = BFGS)

gives a result with the first parameter fairly high away from its true value 
(true values for the three parameters were (0.015, 0.03, 0.2)). I wonder if 
this is due to the autocorrelation in the binary variable (btw. I also tried 
here lm, as the third parameter can be estimated via the correlation, but it 
remains highly dependent on the number of observations).

Any suggestions here?

Best regards

Simon 

On Sep 28, 2012, at 5:02 AM, Rolf Turner rolf.tur...@xtra.co.nz wrote:

 
 I have no idea what your code is doing, nor why you want correlated binary
 variables.  Correlation makes little or no sense in the context of binary 
 random
 variables --- or more generally in the context of discrete random variables.
 
 Be that as it may, it is an easy calculation to show that if X and Y are 
 binary
 random variables both with success probability of 0.5 then cor(X,Y) = 0.2 if
 and only if Pr(X=1 | Y = 1) = 0.6.  So just generate X and Y using that fact:
 
 set.seed(42)
 X - numeric(1000)
 Y - numeric(1000)
 for(i in 1:1000) {
   Y[i] - rbinom(1,1,0.5)
   X[i] - if(Y[i]==1) rbinom(1,1,0.6) else rbinom(1,1,0.4)
 }
 
 # Check:
 cor(X,Y) # Get 0.2012336
 
 Looks about right.  Note that the sample proportions are 0.484 and
 0.485 for X and Y respectively.  These values do not differ significantly
 from 0.5.
 
cheers,
 
Rolf Turner
 
 On 28/09/12 08:26, Simon Zehnder wrote:
 Hi R-fellows,
 
 I am trying to simulate a multivariate correlated sample via the Gaussian 
 copula method. One variable is a binary variable, that should be 
 autocorrelated. The autocorrelation should be rho = 0.2. Furthermore, the 
 overall probability to get either outcome of the binary variable should be 
 0.5.
 Below you can see the R code (I use for simplicity a diagonal matrix in 
 rmvnorm even if it produces no correlated sample):
 
 sampleCop - function(n = 1000, rho = 0.2) {
  
  require(splus2R)
  mvrs - rmvnorm(n + 1, mean = rep(0, 3), cov = diag(3))
  pmvrs - pnorm(mvrs, 0, 1)
  var1 - matrix(0, nrow = n + 1, ncol = 1)
  var1[1] - qbinom(pmvrs[1, 1], 1, 0.5)
  if(var1[1] == 0) var1[nrow(mvrs)] - -1
  for(i in  1:(nrow(pmvrs) - 1)) {
  if(pmvrs[i + 1, 1] = rho) var1[i + 1] - var1[i]
  else var1[i + 1] - var1[i] * (-1)
  }
  sample - matrix(0, nrow = n, ncol = 4)
  sample[, 1] - var1[1:nrow(var1) - 1]
  sample[, 2] - var1[2:nrow(var1)]
  sample[, 3] - qnorm(pmvrs[1:nrow(var1) - 1, 2], 0, 1, 1, 0)
  sample[, 4] - qnorm(pmvrs[1:nrow(var1) - 1, 3], 0, 1, 1, 0)
  
  sample

[R] Autocorrelated binary variable

2012-09-27 Thread Simon Zehnder
Hi R-fellows,

I am trying to simulate a multivariate correlated sample via the Gaussian 
copula method. One variable is a binary variable, that should be 
autocorrelated. The autocorrelation should be rho = 0.2. Furthermore, the 
overall probability to get either outcome of the binary variable should be 0.5. 
Below you can see the R code (I use for simplicity a diagonal matrix in rmvnorm 
even if it produces no correlated sample): 

sampleCop - function(n = 1000, rho = 0.2) {

require(splus2R)
mvrs - rmvnorm(n + 1, mean = rep(0, 3), cov = diag(3))
pmvrs - pnorm(mvrs, 0, 1)
var1 - matrix(0, nrow = n + 1, ncol = 1)
var1[1] - qbinom(pmvrs[1, 1], 1, 0.5)
if(var1[1] == 0) var1[nrow(mvrs)] - -1
for(i in  1:(nrow(pmvrs) - 1)) {
if(pmvrs[i + 1, 1] = rho) var1[i + 1] - var1[i]
else var1[i + 1] - var1[i] * (-1)
}
sample - matrix(0, nrow = n, ncol = 4)
sample[, 1] - var1[1:nrow(var1) - 1]
sample[, 2] - var1[2:nrow(var1)]
sample[, 3] - qnorm(pmvrs[1:nrow(var1) - 1, 2], 0, 1, 1, 0)
sample[, 4] - qnorm(pmvrs[1:nrow(var1) - 1, 3], 0, 1, 1, 0)

sample

}

Now, the code is fine, everything compiles. But when I compute the 
autocorrelation of the binary variable, it is not 0.2, but 0.6. Does anyone 
know why this happens? 

Best Regards

Simon


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Generating an autocorrelated binary variable

2012-09-27 Thread Simon Zehnder
Hi R-fellows,

I am trying to simulate a multivariate correlated sample via the Gaussian 
copula method. One variable is a binary variable, that should be 
autocorrelated. The autocorrelation should be rho = 0.2. Furthermore, the 
overall probability to get either outcome of the binary variable should be 0.5. 
Below you can see the R code (I use for simplicity a diagonal matrix in rmvnorm 
even if it produces no correlated sample): 

sampleCop - function(n = 1000, rho = 0.2) {

require(splus2R)
mvrs - rmvnorm(n + 1, mean = rep(0, 3), cov = diag(3))
pmvrs - pnorm(mvrs, 0, 1)
var1 - matrix(0, nrow = n + 1, ncol = 1)
var1[1] - qbinom(pmvrs[1, 1], 1, 0.5)
if(var1[1] == 0) var1[nrow(mvrs)] - -1
for(i in  1:(nrow(pmvrs) - 1)) {
if(pmvrs[i + 1, 1] = rho) var1[i + 1] - var1[i]
else var1[i + 1] - var1[i] * (-1)
}
sample - matrix(0, nrow = n, ncol = 4)
sample[, 1] - var1[1:nrow(var1) - 1]
sample[, 2] - var1[2:nrow(var1)]
sample[, 3] - qnorm(pmvrs[1:nrow(var1) - 1, 2], 0, 1, 1, 0)
sample[, 4] - qnorm(pmvrs[1:nrow(var1) - 1, 3], 0, 1, 1, 0)

sample

}

Now, the code is fine, everything compiles. But when I compute the 
autocorrelation of the binary variable, it is not 0.2, but 0.6. Does anyone 
know why this happens? 

Best Regards

Simon
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] openmp on R?

2011-09-19 Thread Simon Zehnder
Dear R-Users,

in the configuring file of R on the supercomputer I am using I found a line  

--disable-openmp  (version R 2.13.1)

Is there a possibility to enable this when I let run a BATCH file via R CMD?

Thanks for comments!


Simon

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Problem with memory consuming algorithm

2011-09-16 Thread Simon Zehnder
Hi guyz,

I have serious problems with an algorithm I let run on a supercomputer:

You find the functions under the following URLs:

simuFunctionCaller: http://pastebin.com/6gw2fJFb

calls Function simuFunctionNBM (http://pastebin.com/QeJDUnqx) 

after reading a csv-file ordered like the following:

,name,theta,std_error_theta,theta_phi,std_error_theta_phi,phi,N,sigma_epsilon,sigma_eta
31,25%,0.003044,2.28E-09,0.00742575,4.71E-09,0.00452525,68164.75,0.008706072,0.00774016
32,50%,0.004536,3.22E-09,0.01045,7.94E-09,0.0052762,110150,0.012387386,0.011169953
33,75%,0.007356625,1.24E-08,0.014249075,3.16E-08,0.007518,178840.5,0.018471461,0.016058424

Reading and running without an exception is no problem. But the algorithm needs 
hours to run and seems to use more memory than available:

I have given a time limit of 48 h and a memory limit of 128 GB RAM on the 
supercomputer (it is a cluster, a BATCH System)…..nevertheless the computer 
shuts down the algorithm after a certain time because of memory problems. 

Does anyone of u guys see a problem in the algorithm? Especially with so much 
RAM and time available?

I am very thankful for your suggestions!

Simon
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Exception in NeweyWest - Pre-Whitening necessary?

2011-09-09 Thread Simon Zehnder
Hi guyz,

I have run my algorithm in R (see http://pastebin.com/q84Tujfg) and got the 
following error: 

Error in ar.ols(x, aic = aic, order.max = order.max, na.action = na.action,  :
  'order.max' must be  'n.used'

I am pretty sure, that the error comes from the NeweyWest function in line 45, 
as the NeweyWest function uses the ar.ols() function for pre whitening. Does 
anyone has an idea how to circumvent this? Can I just shut off pre whitening?

I am thankful for every suggestion. 

Simon


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Exception while using NeweyWest function with doMC

2011-08-30 Thread Simon Zehnder
Hi Jay,

first: thank u very much for your comments! U made some very important points 
clear. I tried immediately to write directly the sample function from 

trade-as.big.matrix(matrix(sample(c(1,-1), (N+1)*K, replace=TRUE),ncol=K), 
backingpath=backingpath, backingfile=trade.bin,descriptorfile=trade.desc)

into the big matrix:

trade-big.matrix(sample(c(1,-1), (10+1), replace=TRUE),nrow=(10+1), ncol=10, 
type=double,backingpath=/Users/simon/Documents/R/BigMTest/, 
backingfile=terminaltest.bin, descriptorfile=terminaltest.desc)

But I either get only 1s:

trade[,]  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]111111111 1
 [2,]111111111 1
 [3,]111111111 1
 [4,]111111111 1
 [5,]111111111 1
 [6,]111111111 1
 [7,]111111111 1
 [8,]111111111 1
 [9,]111111111 1
[10,]111111111 1
[11,]111111111 1

or only -1s:

trade[,]  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
 [2,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
 [3,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
 [4,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
 [5,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
 [6,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
 [7,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
 [8,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
 [9,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
[10,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1
[11,]   -1   -1   -1   -1   -1   -1   -1   -1   -1-1

Is there another possibility? In addition I found under ?as.big.matrix, for the 
second example the usage of matrix() inside of as.big.matrix(). But if I 
understood u correctly: usage is possible but does not save memory? I used the 
big.memory package because I always got  - using simply matrix() - exceptions, 
telling me, that memory has reached its limits. After using big.memory all 
worked fine, BUT running http://pastebin.com/UxSkzrae and  
http://pastebin.com/MErGQsQd , I got this: http://pastebin.com/KrEncrSz. It 
seems, as there is a problem with memory allocation inside the underlying 
C-code, maybe a result from my matrix generation inside of the big matrix?

Any suggestions?

 

On Aug 29, 2011, at 6:24 PM, Jay Emerson wrote:

 Simon,
 
 Though we're please to see another use of bigmemory, it really isn't
 clear that it is gaining you
 anything in your example; anything like as.big.matrix(matrix(...))
 still consumes full RAM for both
 the inner matrix() and the new big.matrix -- is the filebacking really
 necessary.  It also doesn't
 appear that you are making use of shared memory, so I'm unsure what
 the gains are.  However,
 I don't have any particular insight as to the subsequent problem with
 NeweyWest (which doesn't
 seem to be using the big.matrix objects).
 
 Jay
 
 --
 Message: 32
 Date: Sat, 27 Aug 2011 21:37:55 +0200
 From: Simon Zehnder simon.zehn...@googlemail.com
 To: r-help@r-project.org
 Subject: [R] Exception while using NeweyWest function with doMC
 Message-ID:
   cagqvrp_gk+t0owbv1ste-y0zafmi9s_zwqrxyxugsui18ms...@mail.gmail.com
 Content-Type: text/plain
 
 Dear R users,
 
 I am using R right now for a simulation of a model that needs a lot of
 memory. Therefore I use the *bigmemory* package and - to make it faster -
 the *doMC* package. See my code posted on http://pastebin.com/dFRGdNrG
 
  snip 
 -
 
 -- 
 John W. Emerson (Jay)
 Associate Professor of Statistics
 Department of Statistics
 Yale University
 http://www.stat.yale.edu/~jay


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Descriptive Stats from Data Frame

2011-08-30 Thread Simon Zehnder
Hi Rich,

I do not know what u really want, because it seems to me, u want to calculate 
the mean of all rows, where the chemical is Arsenic??

But try this to get a little more inside:

mean(chemdata$quant[chemdata$param==Arsenic])

The vector chemdata[chemdata$param==Arsenic,] is a logical vector, returning 
TRUE for every row in which the variable param takes the value Arsenic. Try 
it in your R editor to see it and understand the R concept!

If u now want to get all values of a certain column, given all values have 
Arsenic as param, u just write:

chemdata$COLUMNNAME[chemdata$param==Arsenic]

I do not know if this helps, as it seems to me, that Arsenic only occurs once 
in your frame…..

Good luck Simon



On Aug 30, 2011, at 11:00 PM, Rich Shepard wrote:

  I don't find how to do what I need to do in Dalgaard or 'R Cookbook', so
 I'm asking here.
 
  I have a data frame with water chemistry data and I want to start
 exploring these data. There are three factors (site, date, chemical)
 associated with each measurement. The data frame looks like this:
 
 summary(chemdata)
 site_id.sample_date.param.quant
 BC-0.5|1996-04-19|Arsenic|0.01  :1
 BC-0.5|1996-04-19|Calcium|76.56 :1
 BC-0.5|1996-04-19|Chloride|12   :1
 BC-0.5|1996-04-19|Magnesium|43.23   :1
 BC-0.5|1996-04-19|Sulfate|175   :1
 BC-0.5|1996-04-19|Total Dissolved Solids|460:1
 (Other) :14880
 
  I want first to calculate (and plot) descriptive stats by chemical,
 ignoring site and date and telling R to ignore missing data. (Incorporating
 those factors will occur later.) What I have not been able to figure out is
 how to specify the command to, for example, calculate mean and sd for
 Arsenic. My floundering and thrashing includes attempts like these:
 
 mean(chemdata.param=Arsenic)
 Error in is.numeric(x) : 'x' is missing
 mean(chemdata.quant, param=Arsenic)
 Error in mean(chemdata.quant, param = Arsenic) :
  object 'chemdata.quant' not found
 mean(chemdata$quant, param=Arsenic)
 [1] NA
 Warning message:
 In mean.default(chemdata$quant, param = Arsenic) :
  argument is not numeric or logical: returning NA
 
  As a newcomer to R I've done a lot of reading, yet all the examples use
 nicely structured data to illustrate the point being made. I need to work
 with my data and learn how to specify columns and write correct commands for
 the analyses I need. Please point me in the right direction.
 
 Rich
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Exception while using NeweyWest function with doMC

2011-08-27 Thread Simon Zehnder
Dear R users,

I am using R right now for a simulation of a model that needs a lot of
memory. Therefore I use the *bigmemory* package and - to make it faster -
the *doMC* package. See my code posted on http://pastebin.com/dFRGdNrG

Now, if I use the foreach loop with the addon %do% (for sequential run) I
have no problems at all - only here and there some singularities in
regressor matrices which should be ok.
BUT if I run the loop on multiple cores I get very often a bad exception. I
have posted the exception on http://pastebin.com/eMWF4cu0 The exception
comes from the NeweyWest function loaded within the sandwich library.

I have no clue, what it want to say me and why it is so weirdly printed to
the terminal. I am used to receive here and there errorsbut the messages
never look like this.

Does anyone have a useful answer for me, where to look for the cause of this
weird error?

Here some additional information:

Hardware: MacBook Pro 2.66 GHz Intel Core Duo, 4 GB Memory 1067 MHz DDR3
Software System: Mac Os X Lion 10.7.1 (11B26)
Software App: R64 version 2.11.1 run via Mac terminal

I hope someone has a good suggestion!

Thank u all!

Simon

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.