Re: [R] How to Calculate the Mean by Multiple Groups in R

2023-10-24 Thread Gabor Grothendieck
A variation is to remove Well and then we can use dot to refer to the
remaining columns.

  aggregate(cbind(OD, ODnorm)  ~ . , subset(df, select = - Well), mean)


On Tue, Oct 24, 2023 at 8:32 AM Luigi Marongiu  wrote:
>
> Hello,
> I have a data frame with different groups (Time, Target, Conc) and
> each entry has a triplicate value of the measurements OD and ODnorm.
> How can I merge the triplicates into a single mean value?
> I tried the following:
> ```
> df = data.frame(Time=rep(1, 9), Well=paste("A", 1:9, sep=""),
> OD=c(666, 815, 815, 702, 739, 795, 657, 705, 663),
> Target=rep("BACT", 9),
> Conc=c(1,1,1,2,2,2,3,3,3),
> ODnorm=c(9, 158, 158,  45,  82, 138,   0,  48,   6),
> stringsAsFactors = FALSE)
> aggregate(.~ODnorm, df, mean)
>
> > aggregate(.~ODnorm, df, mean)
>   ODnorm Time Well OD Target Conc
> 1  0   NA   NA NA NA   NA
> 2  6   NA   NA NA NA   NA
> 3  9   NA   NA NA NA   NA
> 4 45   NA   NA NA NA   NA
> 5 48   NA   NA NA NA   NA
> 6 82   NA   NA NA NA   NA
> 7138   NA   NA NA NA   NA
> 8158   NA   NA NA NA   NA
>
>  aggregate(cbind(Time, Target, Conc) ~ ODnorm, df, mean)
>   ODnorm Time Target Conc
> 1  0   NA NA   NA
> 2  6   NA NA   NA
> 3  9   NA NA   NA
> 4 45   NA NA   NA
> 5 48   NA NA   NA
> 6 82   NA NA   NA
> 7138   NA NA   NA
> 8158   NA NA   NA
> ```
>
> Thank you.
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] group consecutive dates in a row

2023-08-07 Thread Gabor Grothendieck
It is best to use Date, rather than POSIXct, class if there are no times.

Use the cumsum expression shown to group the dates and then summarize
each group.

We assume that the dates are already sorted in ascending order.

  library(dplyr)

  mydf <- data.frame(date = as.Date(c("2012-02-05", "2012-02-06",
"2012-02-07", "2012-02-13", "2012-02-21")))

  mydf %>%
group_by(grp = cumsum(c(0, diff(date)) > 1)) %>%
summarize(start = first(date), end = last(date)) %>%
ungroup %>%
select(-grp)
  ## # A tibble: 3 × 2
  ##   start  end
  ##
  ## 1 2012-02-05 2012-02-07
  ## 2 2012-02-13 2012-02-13
  ## 3 2012-02-21 2012-02-21

or with only base R:

  smrz <- function(x) with(x, data.frame(start = min(date), end = max(date)))
  do.call("rbind", by(mydf, cumsum(c(0, diff(mydf$date)) > 1), smrz))
  ##startend
  ## 0 2012-02-05 2012-02-07
  ## 1 2012-02-13 2012-02-13
  ## 2 2012-02-21 2012-02-21


On Mon, Aug 7, 2023 at 12:42 PM Stefano Sofia
 wrote:
>
> Dear R users,
>
> I have a data frame with a single column of POSIXct elements, like
>
>
> mydf <- data.frame(data_POSIX=as.POSIXct(c("2012-02-05", "2012-02-06", 
> "2012-02-07", "2012-02-13", "2012-02-21"), format = "%Y-%m-%d", 
> tz="Etc/GMT-1"))
>
>
> I need to transform it in a two-columns data frame where I can get rid of 
> consecutive dates. It should appear like
>
>
> data_POSIX_init data_POSIX_fin
>
> 2012-02-05 2012-02-07
>
> 2012-02-13 NA
>
> 2012-02-21 NA
>
>
> I started with two "while cycles" and so on, but this is not an efficient way 
> to do it.
>
> Could you please give me an hint on how to proceed?
>
>
> Thank you for your precious attention and help
>
> Stefano
>
>
>  (oo)
> --oOO--( )--OOo--
> Stefano Sofia PhD
> Civil Protection - Marche Region - Italy
> Meteo Section
> Snow Section
> Via del Colle Ameno 5
> 60126 Torrette di Ancona, Ancona (AN)
> Uff: +39 071 806 7743
> E-mail: stefano.so...@regione.marche.it
> ---Oo-oO
>
> 
>
> AVVISO IMPORTANTE: Questo messaggio di posta elettronica può contenere 
> informazioni confidenziali, pertanto è destinato solo a persone autorizzate 
> alla ricezione. I messaggi di posta elettronica per i client di Regione 
> Marche possono contenere informazioni confidenziali e con privilegi legali. 
> Se non si è il destinatario specificato, non leggere, copiare, inoltrare o 
> archiviare questo messaggio. Se si è ricevuto questo messaggio per errore, 
> inoltrarlo al mittente ed eliminarlo completamente dal sistema del proprio 
> computer. Ai sensi dell'art. 6 della DGR n. 1394/2008 si segnala che, in caso 
> di necessità ed urgenza, la risposta al presente messaggio di posta 
> elettronica può essere visionata da persone estranee al destinatario.
> IMPORTANT NOTICE: This e-mail message is intended to be received only by 
> persons entitled to receive the confidential information it may contain. 
> E-mail messages to clients of Regione Marche may contain information that is 
> confidential and legally privileged. Please do not read, copy, forward, or 
> store this message unless you are an intended recipient of it. If you have 
> received this message in error, please forward it to the sender and delete it 
> completely from your computer system.
>
> --
> Questo messaggio  stato analizzato da Libraesva ESG ed  risultato non infetto.
> This message was scanned by Libraesva ESG and is believed to be clean.
>
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Could not read time series data using read.zoo()

2023-08-03 Thread Gabor Grothendieck
The header has white space in it so skip over it, use header = FALSE
and specify the
column headers yourself.  Also use fill=TRUE since the first row does
not have 3 entries.

   # generate test file
  cat("Date Adj Close lret
  02-01-1997 737.01
  03-01-1997 748.03 1.48416235
  06-01-1997 747.65 -0.050813009
  07-01-1997 753.23 0.743567202
  08-01-1997 748.41 -0.64196699
  09-01-1997 754.85 0.856809786
  10-01-1997 759.5 0.614126802", file = "1.csv")

  # test
  library(zoo)
  read.zoo("1.csv", skip = 1, header = FALSE, format = "%m-%d-%Y", fill = TRUE,
 col.names = c(NA, "Adj_Close", "Close"))

On Thu, Aug 3, 2023 at 10:53 AM Christofer Bogaso
 wrote:
>
> Hi,
>
> I have a CSV which contains data like below (only first few rows),
>
> Date Adj Close lret
> 02-01-1997 737.01
> 03-01-1997 748.03 1.48416235
> 06-01-1997 747.65 -0.050813009
> 07-01-1997 753.23 0.743567202
> 08-01-1997 748.41 -0.64196699
> 09-01-1997 754.85 0.856809786
> 10-01-1997 759.5 0.614126802
>
> However when I try to read this data using below code I get error,
>
> read.zoo("1.csv", sep = ',', format = '%d-%m-%Y')
>
> Error reads as,
>
> index has 4500 bad entries at data rows: 1 2 3 4 5 6 7 8 9.
>
> Could you please help to understand why I am getting this error?
>
> > sessionInfo()
>
> R version 4.2.2 (2022-10-31)
>
> Platform: x86_64-apple-darwin17.0 (64-bit)
>
> Running under: macOS Big Sur ... 10.16
>
>
> Matrix products: default
>
> BLAS:   
> /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
>
> LAPACK: 
> /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
>
>
> locale:
>
> [1] C/UTF-8/C/C/C/C
>
>
> attached base packages:
>
> [1] stats graphics  grDevices utils datasets  methods   base
>
>
> other attached packages:
>
> [1] zoo_1.8-12
>
>
> loaded via a namespace (and not attached):
>
> [1] compiler_4.2.2  tools_4.2.2 grid_4.2.2  lattice_0.20-45
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] preserve class in apply function

2023-02-08 Thread Gabor Grothendieck
Also try

apply(Filter(is.numeric, mydf), 1, sum)

On Tue, Feb 7, 2023 at 8:42 AM PIKAL Petr  wrote:
>
> Hi Naresh
>
> If you wanted to automate the function a bit you can use sapply to find
> numeric columns
> ind <- sapply(mydf, is.numeric)
>
> and use it in apply construct
> apply(mydf[,ind], 1, function(row) sum(row))
>  [1]  2.13002569  0.63305300  1.48420429  0.13523859  1.17515873 -0.98531131
>  [7]  0.47044467  0.23914494  0.26504430  0.02037657
>
> Cheers
> Petr
>
> > -Original Message-
> > From: R-help  On Behalf Of Naresh Gurbuxani
> > Sent: Tuesday, February 7, 2023 1:52 PM
> > To: r-help@r-project.org
> > Subject: [R] preserve class in apply function
> >
> >
> > > Consider a data.frame whose different columns have numeric, character,
> > > and factor data.  In apply function, R seems to pass all elements of a
> > > row as character.  Is it possible to preserve numeric class?
> > >
> > >> mydf <- data.frame(x = rnorm(10), y = runif(10))
> > >> apply(mydf, 1, function(row) {row["x"] + row["y"]})
> > > [1]  0.60150197 -0.74201827  0.80476392 -0.59729280 -0.02980335
> > 0.31351909
> > > [7] -0.63575990  0.22670658  0.55696314  0.39587314
> > >> mydf[, "z"] <- sample(letters[1:3], 10, replace = TRUE)
> > >> apply(mydf, 1, function(row) {row["x"] + row["y"]})
> > > Error in row["x"] + row["y"] (from #1) : non-numeric argument to binary
> > operator
> > >> apply(mydf, 1, function(row) {as.numeric(row["x"]) +
> > as.numeric(row["y"])})
> > > [1]  0.60150194 -0.74201826  0.80476394 -0.59729282 -0.02980338
> > 0.31351912
> > > [7] -0.63575991  0.22670663  0.55696309  0.39587311
> > >> apply(mydf[,c("x", "y")], 1, function(row) {row["x"] + row["y"]})
> > > [1]  0.60150197 -0.74201827  0.80476392 -0.59729280 -0.02980335
> > 0.31351909
> > > [7] -0.63575990  0.22670658  0.55696314  0.39587314
> >
> > __
> > 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] R emulation of FindRoot in Mathematica

2023-01-19 Thread Gabor Grothendieck
If the equations are in the form shown in your post then take the log of
both sides, expand the logs and replace log(whatever) with new variables so
now the equations are in linear form and are easy to solve.

__
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] Function for Distribution Fitting

2022-10-26 Thread Gabor Grothendieck
A Cullen & Frey graph (fitdistrplus::descdist) can be used to compare certain
common distributions.

On Wed, Oct 26, 2022 at 9:42 AM Paul Bernal  wrote:
>
> Dear friends from the R community,
>
> Hope you are all doing great. So far, whenever I need to perform
> distribution fitting on a particular dataset, I make use of R package
> fitdistrplus.
>
> However, distribution fitting using the fitdist() function from
> fitdistrplus is rather manual (you need to specify which probability
> distribution you are trying to fit), and it would be more convenient if the
> function tested all possible probability distributions (or a bunch of them)
> and then estimates the parameters and suggests the most suitable
> distribution.
>
> Is there any package besides fitdistrplus that does allow automatic
> distribution fitting?
>
> Best regards,
> Paul
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Adding page numbers to existing PDFs

2022-10-21 Thread Gabor Grothendieck
Create a pdf using latex that has only page numbers and then
superimpose that with your pdf using the free utility pdftk.  The
animation R package has an interface to pdftk. Google to locate
pdftk and again to locate instructions.

There are also freeware GUI Windows programs  that are easy
to use.  This has nothing to do with R.  The one I have used is
A-PDF Number.

On Fri, Oct 21, 2022 at 1:49 PM Dennis Fisher  wrote:
>
> R 4.2.1
> OS X
>
> Colleagues
>
> I have multipage PDF files that were created in R — the files do NOT have 
> page numbers.  I would like to add page numbers after the fact, i.e., read a 
> file into R, add margin text, then output as a PDF.
> Can this be done in R (base R or a package)?
> It can be done in Python using reportlab.pdfgen — however, the file size 
> increases prohibitively.
>
> Dennis
>
> Dennis Fisher MD
> P < (The "P Less Than" Company)
> Phone / Fax: 1-866-PLessThan (1-866-753-7784)
> www.PLessThan.com
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Unintended behaviour of stats::time not returning integers for the first cycle

2022-10-18 Thread Gabor Grothendieck
The zoo package implements tolerances internally. Converting the ts
object to zoo:

  library(zoo)
  identical(as.integer(time(as.zoo(x))), true.year)
  ## [1] TRUE

On Sat, Oct 15, 2022 at 3:26 AM Andreï V. Kostyrka
 wrote:
>
> Dear all,
>
>
>
> I was using stats::time to obtain the year as a floor of it, and
> encountered a problem: due to a rounding error (most likely due to its
> reliance on the base::seq.int internally, but correct me if I am wrong),
> the actual number corresponding to the beginning of a year X can still be
> (X-1).999..., resulting in the following undesirable behaviour.
>
>
>
> One of the simplest ways of getting the year from a ts object is
> floor(time(...)). However, if the starting time cannot be represented
> nicely as a power of 2, then, supposedly integer time does not have a
> .00... mantissa:
>
>
>
> x <- ts(2:252, start = c(2002, 2), freq = 12)
>
> d <- seq.Date(as.Date("2002-02-01"), to = as.Date("2022-12-01"), by =
> "month")
>
> true.year <- rep(2002:2022, each = 12)[-1]
>
> wrong.year <- floor(as.numeric(time(x)))
>
> tail(cbind(as.character(d), true.year, wrong.year), 15) # Look at 2022-01-01
>
> print(as.numeric(time(x))[240], 20) # 2021.7726, the floor of
> which is 2021
>
>
>
> Yes, I have read the 'R inferno' book and know the famous '0.3 != 0.7 -
> 0.4' example, but I believe that the expected / intended behaviour would be
> actually returning round years for the first observation in a year. This
> could be achieved by rounding the near-integer time to integers.
>
>
>
> Since users working with dates are expecting to get exact integer years for
> the first cycle of a ts, this should be changed. Thank you in advance for
> considering a fix.
>
>
>
> Yours sincerely,
>
> Andreï V. Kostyrka
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] rmoving dates from an xts object...

2022-07-10 Thread Gabor Grothendieck
Look at the examples at the end of ?xts for more info.

  library(quantmod)
  getSymbols("AAPL")

  class(AAPL)
  ## [1] "xts" "zoo"

  range(time(AAPL))
  ## [1] "2007-01-03" "2022-07-08"

  # everything up to indicated date
  a1 <- AAPL["/2018-02-01"]

  # remove non consecutive dates
  d <- as.Date(c("2022-07-01", "2022-07-06")) # dates to remove
  a2 <- AAPL[ ! time(AAPL) %in% d ]

On Sun, Jul 10, 2022 at 11:42 AM akshay kulkarni  wrote:
>
> Dear members,
>  I have OHLC data of 500 stocks: OHLCData and dates. 
> These are of xts object. I want to do the following:
>
>
>   1.  I want to remove a contiguous set of dates from the xts object. for 
> example, I want to remove data from the OHLC data from "2022-20-7"
>   2.   to "2018-2-2", continuously.
>   3.
>
>2. I want to remove a set of dates, which are not contiguous.
>
> Any idea on how to accomplish it? I can write an intricate for loop, but any 
> other method? Does an xts object behave like an atomic vector : 
> OHLCData[[i]][-dates[[i]]] ?
>
> Many thanks in advance
>
> THanking you,
> Yours sincerely,
> AKSHAY M KULKARNI
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Restoration of "rite" package of R as R-script editor or the like

2022-05-24 Thread Gabor Grothendieck
Soren Hojsgaard had also written a tcltk based editor R package and although
it is deprecated too it may be simpler than rite in which case it is
more likely to work.  Contact him (he has other packages on CRAN)
and see if he would send it to you and if it fits your needs.

On Sun, May 22, 2022 at 12:54 AM Akhilesh Singh
 wrote:
>
> Dear Sir,
>
> I am a professor at Indira Gandhi Agricultural University, Raipur,
> Chhattisgarh, India. I have been teaching statistics in my department since
> 1986.
>
> Almost since the last 10 years, I have included R capabilities in doing
> practical-classes of various statistics courses in my department. I use the
> default R-script editor that comes with RGUI.exe for that purpose. *But,
> this text editor is too simple without syntax highlighting etc. And RStudio
> is too heavy for the UG/PG level students who are basically from
> agriculture stream.*
>
> Meanwhile, I came across a "rite" package which was earlier maintained
> by Thomas
> J. Leeper . Now, it has been put into archives and
> not being maintained by anyone.
>
> At present, I am using the "rite" package for students' practical classes.
> But, I am facing difficulties using it. Sometimes, it hangs unnecessarily,
> and some of its features do not work, probably because it is not being
> maintained anymore by CRAN.
>
> Therefore, to remove this difficulty, the "rite" package must be maintained.
>
>
> Sir, "rite" package is a very light weight editor for R codes suitable for
> teaching to beginner students, with syntax highlighting facility and with
> the facility of generating reports into HTML, markdown etc.
>
> Therefore, I request you ,*one and all who are in control, to kindly make
> sure to include the "rite" package into the regular packages with proper
> maintenance, which will immensely help the beginner students to learn R*.
>
> Or,
>
> Kindly make sure that the *default R-script-editor gets the syntax
> highlighting features and capability of creating basic reports*.
>
> Or,
>
> Kindly make a lighter version of RStudio *for beginner student-level only,
> like the "rite" package,*  for executing R-scripts and for generating some
> basic reports.
>
> *In fact, the lack of a proper R-script editor is a major handicap among
> the beginner students to learn R.*
>
> Therefore, I earnestly request your kind attention to this problem.
>
> Thaing you.
>
>
> Dr. A.K. Singh
> Professor and Ex-Head (Agricultural Statistics)
> Department of Agricultural Statistics and Social Science (L)
> Indira Gandhi Krishi Vishwavidyalaya, Raipur-492 012,
> Chhattisgarh, India
> Mobile: +918770625795
> Email: akhileshsingh.i...@gmail.com
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Find "undirected" duplicates in a tibble

2021-08-20 Thread Gabor Grothendieck
Since you are dealing with graphs you could consider using
the igraph package.  This is more involved than needed for
what you are
asking but it might be useful for other follow on calculations.
We first define a 2 column matrix of edges, then convert it to
an igraph and simplify it to remove duplicate edges giving g.
At the end we get an edgelist back.

  library(igraph)
  m <- matrix(c(1, 2, 6, 6, 4, 9, 1, 5, 2, 1, 8, 7, 5, 10, 6, 10), 8, 2)
  g <- m |>
graph_from_edgelist(directed = FALSE) |>
simplify()

  plot(g)

  g |>
get.edgelist() |>
as.data.frame()



On Fri, Aug 20, 2021 at 5:00 AM Kimmo Elo  wrote:
>
> Hi!
>
> I am working with a large network data consisting of source-target
> pairs stored in a tibble. Now I need to transform the directed dataset
> to an undirected network data. This means, I need to keep only one
> instance for pairs with the same "nodes". In other words, if my data
> has one row with A (source) and B (target) and one with B (source) and
> A (target), only the pair A-B should be kept.
>
> Here an example how I have solved this problem so far:
>
> --- snip ---
>
> # Create some data
> x<-tibble(Source=rep(1:3,4), Target=c(rep(1,3),rep(2,3),rep(3,3),rep(4,3)))
> x   # print original data
>
> # Remove "undirected" duplicates
> x<-x %>% mutate(pair=mapply(function(x,y)
> paste0(sort(c(x,y)),collapse="-"), Source, Target)) %>% distinct(pair,
> .keep_all = T) %>% mutate(Source=sapply(pair, function(x)
> unlist(strsplit(x, split="-"))[1]), Target=sapply(pair, function(x)
> unlist(strsplit(x, split="-"))[2])) %>% select(-pair)
>
> x   # print cleaned data
>
> --- snip ---
>
> The good thing with my own solution is that it allows the creation of
> weighted pairs as well. One just needs to replace 'distinct(pair,
> .keep_all=T)' with 'count(pair)'.
>
> I have done a lot of searching but not found any function providing
> this functionality. Does someone know an alternative, maybe a more
> effective function/solution?
>
> Best,
>
> Kimmo Elo
>
>
> --
> Dr. Kimmo Elo
> Senior researcher in European Studies
> =
> University of Turku
> Centre for Parliamentary Studies
> Finland
> E-mail: kimmo@utu.fi
> =
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Testing optimization solvers with equality constraints

2021-05-27 Thread Gabor Grothendieck
In case it is of interest this problem can be solved with an
unconstrained optimizer,
here optim, like this:

proj <- function(x) x / sqrt(sum(x * x))
opt <- optim(c(0, 0, 1), function(x) f(proj(x)))
proj(opt$par)
## [1] 5.388907e-09 7.071068e-01 7.071068e-01

On Fri, May 21, 2021 at 11:01 AM Hans W  wrote:
>
> Just by chance I came across the following example of minimizing
> a simple function
>
> (x,y,z) --> 2 (x^2 - y z)
>
> on the unit sphere, the only constraint present.
> I tried it with two starting points, x1 = (1,0,0) and x2 = (0,0,1).
>
> #-- Problem definition in R
> f = function(x)  2 * (x[1]^2 - x[2]*x[3])   # (x,y,z) |-> 2(x^2 -yz)
> g = function(x)  c(4*x[1], 2*x[3], 2*x[2])  # its gradient
>
> x0 = c(1, 0, 0); x1 = c(0, 0, 1)# starting points
> xmin = c(0, 1/sqrt(2), 1/sqrt(2))   # true minimum -1
>
> heq = function(x)  1-x[1]^2-x[2]^2-x[3]^2   # staying on the sphere
> conf = function(x) {# constraint function
> fun = x[1]^2 + x[2]^2 + x[3]^2 - 1
> return(list(ceq = fun, c = NULL))
> }
>
> I tried all the nonlinear optimization solvers in R packages that
> allow for equality constraints: 'auglag()' in alabama, 'solnl()' in
> NlcOptim, 'auglag()' in nloptr, 'solnp()' in Rsolnp, or even 'donlp2()'
> from the Rdonlp2 package (on R-Forge).
>
> None of them worked from both starting points:
>
> # alabama
> alabama::auglag(x0, fn = f, gr = g, heq = heq)  # right (inaccurate)
> alabama::auglag(x1, fn = f, gr = g, heq = heq)  # wrong
>
> # NlcOptim
> NlcOptim::solnl(x0, objfun = f, confun = conf)  # wrong
> NlcOptim::solnl(x1, objfun = f, confun = conf)  # right
>
> # nloptr
> nloptr::auglag(x0, fn = f, heq = heq)   # wrong
> # nloptr::auglag(x1, fn = f, heq = heq) # not returning
>
> # Rsolnp
> Rsolnp::solnp(x0, fun = f, eqfun = heq) # wrong
> Rsolnp::solnp(x1, fun = f, eqfun = heq) # wrong
>
> # Rdonlp2
> Rdonlp2::donlp2(x0, fn = f, nlin = list(heq),   # wrong
>nlin.lower = 0, nlin.upper = 0)
> Rdonlp2::donlp2(x1, fn = f, nlin = list(heq),   # right
>nlin.lower = 0, nlin.upper = 0)  # (fast and exact)
>
> The problem with starting point x0 appears to be that the gradient at
> that point, projected onto the unit sphere, is zero. Only alabama is
> able to handle this somehow.
>
> I do not know what problem most solvers have with starting point x1.
> The fact that Rdonlp2 is the fastest and most accurate is no surprise.
>
> If anyone with more experience with one or more of these packages can
> give a hint of what I made wrong, or how to change calling the solver
> to make it run correctly, please let me know.
>
> Thanks  -- HW
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] xyplot.zoo trouble with the index when I use microseconds (type POSIXct)

2021-02-22 Thread Gabor Grothendieck
Also, it might be better to simply use seconds as the time. In that
case both plot and xyplot run without error.

 time(DF.w) <- as.POSIXlt(time(DF.w))$sec
 # now run plot or xyplot as before

On Mon, Feb 22, 2021 at 11:56 AM Gabor Grothendieck
 wrote:
>
> I assume that this is a lattice problem.  Replacing xyplot with plot
> and using all the same arguments there is no error.
>
>
> On Mon, Feb 22, 2021 at 11:26 AM Laurent Rhelp  wrote:
> >
> > Dear R-Help-List,
> >
> > I have to process time series with a sampling frequency of 1 MHz.
> > I use the POSIXct format for the date-times with microsecond in a zoo
> > object and the xyplot.zoo function to do the graphs.
> > As I show in the below example I had a trouble to plot the labels on the
> > x-axis with an error message. I found a solution but I would like to
> > know if I miss something.
> > The link
> > https://stackoverflow.com/questions/7726034/how-r-formats-posixct-with-fractional-seconds
> >
> > helped me to understand how to print the POSIXct value to see the
> > microseconds thanks to the function myformat.POSIXct:
> >
> > myformat.POSIXct <- function(x, digits=6)
> > {
> >x2 <- round(unclass(x), digits)
> >attributes(x2) <- attributes(x)
> >x <- as.POSIXlt(x2,origin="1970-01-01",tz="GMT")
> >x$sec <- round(x$sec, digits)
> >format.POSIXlt(x, paste("%Y-%m-%d %H:%M:%OS",digits,sep=""),tz="GMT")
> >
> > }
> >
> >
> > ## The example giving the error message:
> >
> > library(lattice)
> > library(zoo)
> >
> > ##
> > options(digits = 16) # to see all the digits on the screen
> > options(digits.secs = 6) # to see the microseconds
> > # mock data
> > # a sine with a frequency f0 and two others with a delay
> > Fs <- 1e+6 # sampling frequency 1 MHz
> > Ts <- 1/Fs
> > # frequency of the sinus
> > f0 <- 10
> > t0 <- 1/f0
> > time <- seq(0, length = 1000, by = Ts)
> > A1 <- 1
> > y1 <- A1 * sin(2*pi*f0*time)
> > y2 <- 2 * A1 * sin(2*pi*f0*(time+0.02))
> > y3 <- 3 * A1 * sin(2*pi*f0*(time+0.05))
> > ## creation of a dataframe:
> > ##
> > DF <- data.frame( time = time, y1 = y1, y2 = y2, y3 = y3)
> > # Since I want to allow for the datetime POSIXct format on the x-axis
> > # for the plot I transform my dataframe in a zoo object
> > #
> > # say that my acquisition began at "2021-02-08 09:15:50.00"
> > #
> > mystart <- as.POSIXct("2021-02-08 09:15:50.00", format = "%Y-%m-%d
> > %H:%M:%OS",tz="GMT")
> > mystart
> > # To see the correct datetime I use the myformat.POSIXct function
> > myformat.POSIXct(mystart)
> > ##
> > ## using the method seq.POSIXct as following doesn't work:
> > ## mydatetime <- seq( mystart , length = nrow(DF), by = "0.01 sec")
> > ## head( myformat.POSIXct(mydatetime) )
> > ## if I use the following command it works:
> > mydatetime <- seq( mystart , length = nrow(DF), by = 0.01)
> > head( myformat.POSIXct(mydatetime) )
> > ## I do the zoo object:
> > DF.z <- zoo(DF[,-1],order.by = mydatetime)
> > ## We don't see the correct value for the index:
> > head(DF.z)
> > # timey1 y2 y3
> > # 2021-02-08 09:15:50.00 0e+00 0.000e+00
> > 1.902113032590307e+00  3.673819061467132e-16
> > # 2021-02-08 09:15:50.00 1e-06 5.877852522924730e-01
> > 1.902113032590307e+00 -1.763355756877419e+00
> > # 2021-02-08 09:15:50.01 2e-06 9.510565162951535e-01
> > 1.175570504584947e+00 -2.853169548885460e+00
> > # 2021-02-08 09:15:50.03 3e-06 9.510565162951536e-01
> > 1.133099690464601e-15 -2.853169548885460e+00
> > # 2021-02-08 09:15:50.04 4e-06 5.877852522924736e-01
> > -1.175570504584946e+00 -1.763355756877420e+00
> > # 2021-02-08 09:15:50.05 5e-06 5.665498452323003e-16
> > -1.902113032590306e+00 -3.399299071393802e-15
> > # If I use myformat.POSIXct I see that the index is correct in the
> > object DF:
> > head(myformat.POSIXct(index(DF.z)))
> > ## and when I plot I have an error:
> > xyplot(   DF.z
> >, screens = c(1,1,1)
> >, type = "l"
> >, col = c("red","blue","black")
> >
> > )
> >
> > # Error in prettyDate_TMP(x, ...) : range too small for min.n
> >
> > # if I process by hand the plot of the labels on the x-axis it wo

Re: [R] xyplot.zoo trouble with the index when I use microseconds (type POSIXct)

2021-02-22 Thread Gabor Grothendieck
I assume that this is a lattice problem.  Replacing xyplot with plot
and using all the same arguments there is no error.


On Mon, Feb 22, 2021 at 11:26 AM Laurent Rhelp  wrote:
>
> Dear R-Help-List,
>
> I have to process time series with a sampling frequency of 1 MHz.
> I use the POSIXct format for the date-times with microsecond in a zoo
> object and the xyplot.zoo function to do the graphs.
> As I show in the below example I had a trouble to plot the labels on the
> x-axis with an error message. I found a solution but I would like to
> know if I miss something.
> The link
> https://stackoverflow.com/questions/7726034/how-r-formats-posixct-with-fractional-seconds
>
> helped me to understand how to print the POSIXct value to see the
> microseconds thanks to the function myformat.POSIXct:
>
> myformat.POSIXct <- function(x, digits=6)
> {
>x2 <- round(unclass(x), digits)
>attributes(x2) <- attributes(x)
>x <- as.POSIXlt(x2,origin="1970-01-01",tz="GMT")
>x$sec <- round(x$sec, digits)
>format.POSIXlt(x, paste("%Y-%m-%d %H:%M:%OS",digits,sep=""),tz="GMT")
>
> }
>
>
> ## The example giving the error message:
>
> library(lattice)
> library(zoo)
>
> ##
> options(digits = 16) # to see all the digits on the screen
> options(digits.secs = 6) # to see the microseconds
> # mock data
> # a sine with a frequency f0 and two others with a delay
> Fs <- 1e+6 # sampling frequency 1 MHz
> Ts <- 1/Fs
> # frequency of the sinus
> f0 <- 10
> t0 <- 1/f0
> time <- seq(0, length = 1000, by = Ts)
> A1 <- 1
> y1 <- A1 * sin(2*pi*f0*time)
> y2 <- 2 * A1 * sin(2*pi*f0*(time+0.02))
> y3 <- 3 * A1 * sin(2*pi*f0*(time+0.05))
> ## creation of a dataframe:
> ##
> DF <- data.frame( time = time, y1 = y1, y2 = y2, y3 = y3)
> # Since I want to allow for the datetime POSIXct format on the x-axis
> # for the plot I transform my dataframe in a zoo object
> #
> # say that my acquisition began at "2021-02-08 09:15:50.00"
> #
> mystart <- as.POSIXct("2021-02-08 09:15:50.00", format = "%Y-%m-%d
> %H:%M:%OS",tz="GMT")
> mystart
> # To see the correct datetime I use the myformat.POSIXct function
> myformat.POSIXct(mystart)
> ##
> ## using the method seq.POSIXct as following doesn't work:
> ## mydatetime <- seq( mystart , length = nrow(DF), by = "0.01 sec")
> ## head( myformat.POSIXct(mydatetime) )
> ## if I use the following command it works:
> mydatetime <- seq( mystart , length = nrow(DF), by = 0.01)
> head( myformat.POSIXct(mydatetime) )
> ## I do the zoo object:
> DF.z <- zoo(DF[,-1],order.by = mydatetime)
> ## We don't see the correct value for the index:
> head(DF.z)
> # timey1 y2 y3
> # 2021-02-08 09:15:50.00 0e+00 0.000e+00
> 1.902113032590307e+00  3.673819061467132e-16
> # 2021-02-08 09:15:50.00 1e-06 5.877852522924730e-01
> 1.902113032590307e+00 -1.763355756877419e+00
> # 2021-02-08 09:15:50.01 2e-06 9.510565162951535e-01
> 1.175570504584947e+00 -2.853169548885460e+00
> # 2021-02-08 09:15:50.03 3e-06 9.510565162951536e-01
> 1.133099690464601e-15 -2.853169548885460e+00
> # 2021-02-08 09:15:50.04 4e-06 5.877852522924736e-01
> -1.175570504584946e+00 -1.763355756877420e+00
> # 2021-02-08 09:15:50.05 5e-06 5.665498452323003e-16
> -1.902113032590306e+00 -3.399299071393802e-15
> # If I use myformat.POSIXct I see that the index is correct in the
> object DF:
> head(myformat.POSIXct(index(DF.z)))
> ## and when I plot I have an error:
> xyplot(   DF.z
>, screens = c(1,1,1)
>, type = "l"
>, col = c("red","blue","black")
>
> )
>
> # Error in prettyDate_TMP(x, ...) : range too small for min.n
>
> # if I process by hand the plot of the labels on the x-axis it works:
> #
> myend <- tail(mydatetime,1)
>
> myticks <- seq( mystart , to = myend , length = 5)
> mylabels <- format(myticks,"%H:%M:%OS")
>
> xyplot(  DF.z
>   , screens = c(1,1,1)
>   , type = "l"
>   , col = c("red","blue","black")
>   , scales = list(
>   y= list(relation = "free", abbreviate=TRUE),
>   x = list( at = myticks, labels = mylabels
> , rot = 45, cex = 0.5)
> )
> )
> # The microseconds are well taken into account with the window function
> # if I want to plot only 100 microseconds but there is of course the same
> # trouble for the plot
>
> myend <- as.POSIXct("2021-02-08 09:15:50.000100", format = "%Y-%m-%d
> %H:%M:%OS",tz="GMT")
> myformat.POSIXct(myend)
>
> DF.w <- window( DF.z ,start = mystart, end = myend)
>
> myticks <- seq( mystart , to = myend , length = 5)
> mylabels <- format(myticks,"%H:%M:%OS")
>
> xyplot(  DF.w
>   , screens = c(1,1,1)
>   , type = "l"
>   , col = c("red","blue","black")
> )
>
> ## Error in prettyDate_TMP(x, ...) : range too small for min.n
>
> xyplot(  DF.w
>   , screens = c(1,1,1)
>   , type = "l"
>   

Re: [R] Finance & R

2020-12-24 Thread Gabor Grothendieck
I personally use the packages you mention and
would describe them as well thought out and relatively
bug free through widespread use and years of improvement
and maintenance.   There are literally dozens of other packages
that depend on these packages so if you use them you will
also be able to seamlessly use a wide range of other
packages too.


On Wed, Dec 23, 2020 at 12:58 PM Ben van den Anker via R-help
 wrote:
>
>
>
> Hello everyone,
> Could anyonre recommend some good resources for finance applications in R? I 
> find the packages quantmod, TTR and PerformanceAnalytics a bit outdated. 
> There must be something more recent on the market. Any suggestions will be 
> much appreciated!
> Cheers,
> Ben van den Anker
>
>
>
>
>
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] nls() syntax

2020-12-11 Thread Gabor Grothendieck
The start= argument should be as follows:

 nls(y ~ x/(x - a[z]),start=list(a = strt),data=xxx)

On Fri, Dec 11, 2020 at 6:51 PM Rolf Turner  wrote:
>
>
>
> I want to fit a model y = x/(x-a) where the value of a depends
> on the level of a factor z.  I cannot figure out an appropriate
> syntax for nls().  The "parameter" a (to be estimated) should be a
> vector of length equal to the number of levels of z.
>
> I tried:
>
> strt <- rep(3,length(levels(z))
> names(strt=levels(z)
> fit <- nls(y ~ x/(x - a[z]),start=strt,data=xxx)
>
> but of course got an error:
>
> > Error in nls(y ~ x/(x - a[z]), start = strt, data = xxx) :
> >   parameters without starting value in 'data': a
>
> I keep thinking that there is something obvious that I should
> be doing, but I can't work out what it is.
>
> Does there *exist* an appropriate syntax for doing what I want
> to do?  Can anyone enlighten me?  The data set "xxx" is given
> in dput() form at the end of this message.
>
> cheers,
>
> Rolf Turner
>
> --
> Honorary Research Fellow
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>
> Data set "xxx":
>
> structure(list(x = c(30, 40, 50, 60, 70, 80, 90, 100, 110, 120,
> 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250,
> 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160,
> 170, 180, 190, 200, 210, 220, 230, 240, 250, 30, 40, 50, 60,
> 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190,
> 200, 210, 220, 230, 240, 250, 30, 40, 50, 60, 70, 80, 90, 100,
> 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230,
> 240, 250, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140,
> 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250), y = c(1.27,
> 1.16, 1.19, 1.15, 1.09, 1.07, 1.07, 1.05, 1.07, 1.03, 1.05, 1.07,
> 1.06, 1.03, 1.05, 1.04, 1.03, 1.03, 1.03, 1.02, 1.02, 1.01, 1.01,
> 1.21, 1.15, 1.1, 1.1, 1.06, 1.06, 1.05, 1.03, 1.07, 1.04, 1.04,
> 1.02, 1.04, 1.02, 1.04, 1.03, 1.01, 1.03, 1.01, 1, 1.02, 1.03,
> 1.02, 1.42, 1.27, 1.23, 1.14, 1.17, 1.08, 1.11, 1.06, 1.07, 1.08,
> 1.06, 1.07, 1.04, 1.03, 1.07, 1.04, 1.03, 1.03, 1.03, 1.04, 1.03,
> 1.03, 1.04, 1.85, 1.41, 1.35, 1.21, 1.22, 1.15, 1.14, 1.07, 1.1,
> 1.09, 1.1, 1.09, 1.08, 1.08, 1.09, 1.09, 1.07, 1.06, 1.03, 1.08,
> 1.05, 1.02, 1.05, 1.99, 1.6, 1.44, 1.4, 1.24, 1.3, 1.21, 1.23,
> 1.18, 1.18, 1.12, 1.15, 1.09, 1.07, 1.13, 1.1, 1.05, 1.13, 1.09,
> 1.03, 1.11, 1.07, 1.05), z = structure(c(1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
> 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
> 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
> 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
> 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
> 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("p1",
> "p2", "p3", "p4", "p5"), class = "factor")), class = "data.frame", row.names 
> = c(NA,
> -115L))
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] formula mungeing

2020-10-23 Thread Gabor Grothendieck
Recursively walk the formula performing the replacement:

g <- function(e, ...) {
if (length(e) > 1) {
if (identical(e[[2]], as.name(names(list(...) {
  e <- eval(e, list(...))
}
if (length(e) > 1) for (i in 1:length(e)) e[[i]] <- Recall(e[[i]], ...)
}
e
}

g(f, lambdas = 2:3)
## y ~ qss(x, lambda = 2L) + qss(z, 3L) + s

On Fri, Oct 23, 2020 at 9:33 AM Koenker, Roger W  wrote:
>
> Suppose I have a formula like this:
>
> f <- y ~ qss(x, lambda = lambdas[1]) + qss(z, lambdas[2]) + s
>
> I’d like a function, g(lambdas, f)  that would take g(c(2,3), f) and produce 
> the new
> formula:
>
> y ~ qss(x, lambda = 2) + qss(z, 3) + s
>
> For only two qss terms I have been using
>
> g <- function(lambdas, f){
> F <- deparse(f)
> F <- gsub("lambdas\\[1\\]",lambdas[1],F)
> F <- gsub("lambdas\\[2\\]",lambdas[2],F)
> formula(F)
> }
> but this is ugly and doesn’t extend nicely to more qss terms.  Isn’t there 
> some
> bquote() magic that can be invoked?  Or something else entirely?
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Help with read.csv.sql()

2020-07-30 Thread Gabor Grothendieck
Probably simplest to assign the names afterwards as others have
suggested but it could be done like this:

  library(sqldf)
  write.csv(BOD, "BOD.csv", quote = FALSE, row.names = FALSE)  # test data

  read.csv.sql("BOD.csv", "select Time as Time2, demand as demand2 from file")

giving the column names Time2 and demand2 rather than the original column names.

Time2 demand2
  1 1 8.3
  2 210.3
  3 319.0
  4 416.0
  5 515.6
  6 719.8

On Fri, Jul 17, 2020 at 9:28 PM H  wrote:
>
> I have created a dataframe with columns that are characters, integers and 
> numeric and with column names assigned by me. I am using read.csv.sql() to 
> read portions of a number of large csv files into this dataframe, each csv 
> file having a header row with columb names.
>
> The problem I am having is that the csv files have header rows with column 
> names that are slightly different from the column names I have assigned in 
> the dataframe and it seems that when I read the csv data into the dataframe, 
> the column names from the csv file replace the column names I chose when 
> creating the dataframe.
>
> I have been unable to figure out if it is possible to assign column names of 
> my choosing in the read.csv.sql() function? I have tried various variations 
> but none seem to work. I tried colClasses = c() but that did not work, I 
> tried field.types = c(...) but could not get that to work either.
>
> It seems that the above should be feasible but I am missing something? Does 
> anyone know?
>
> A secondary issue is that the csv files have a column with a date in 
> mm/dd/ format that I would like to make into a Date type column in my 
> dataframe. Again, I have been unable to find a way - if at all possible - to 
> force a conversion into a Date format when importing into the dataframe. The 
> best I have so far is to import is a character column and then use as.Date() 
> to later force the conversion of the dataframe column.
>
> Is it possible to do this when importing using read.csv.sql()?
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] sqldf and number of records affected

2020-06-11 Thread Gabor Grothendieck
There is no real difference between your example and the example I provided.
Both use a data.frame in R and both work for me under R 3.5 with RSQLite 2.2.0
See log below.

Note that there is a bug in R 4.0 related to tcltk that could possibly affect
sqldf as it uses tcltk.  A fix has been announced for R 4.0.2.

> library(sqldf)
Loading required package: gsubfn
Loading required package: proto
Loading required package: RSQLite
> con <- data.frame(V1 = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
> sqldf()

  Path: :memory:
  Extensions: TRUE
> sqldf(c("pragma count_changes = 1", "update con set V1 = 0 where V1 > 5 "))
  rows updated
15
Warning message:
In result_fetch(res@ptr, n = n) :
  SQL statements must be issued with dbExecute() or dbSendStatement()
instead of dbGetQuery() or dbSendQuery().
> ans <- sqldf("select * from main.con")
> sqldf()
NULL
> ans
   V1
1   1
2   2
3   3
4   4
5   5
6   0
7   0
8   0
9   0
10  0
> R.version.string
[1] "R version 3.5.3 (2019-03-11)"
> packageVersion("sqldf")
[1] ‘0.4.11’
> packageVersion("RSQLite")
[1] ‘2.2.0’
> packageVersion("DBI")
[1] ‘1.1.0’



On Thu, Jun 11, 2020 at 10:06 AM Ravi Jeyaraman  wrote:
>
> Thanks for the response Gabor.  Looks like the below example will work when 
> using SQLite, but in my case I'm just creating a dataframe in R and trying to 
> update it using sqldf as below and it doesn't seem to work ...
>
> con <- data.frame(V1 = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
> sqldf()
> sqldf(c("pragma count_changes = 1", "update con set V1 = 0 where V1 > 5 "))
> ans <- sqldf("select * from main.con")
> sqldf()
>
> -Original Message-
> From: Gabor Grothendieck [mailto:ggrothendi...@gmail.com]
> Sent: Thursday, June 11, 2020 9:12 AM
> To: Ravi Jeyaraman 
> Cc: r-help@r-project.org
> Subject: Re: [R] sqldf and number of records affected
>
> Here is an example.  Ignore the warning or use the workaround discussed here
> https://github.com/ggrothendieck/sqldf/issues/40
> to avoid the warning.
>
>   library(sqldf)
>   sqldf()  # use same connection until next sqldf()
>   sqldf(c("pragma count_changes = 1", "update BOD set demand = 99 where Time 
> > 4"))
>   sqldf("select * from main.BOD")
>   sqldf()
>
>
> On Thu, Jun 11, 2020 at 9:01 AM Ravi Jeyaraman  wrote:
> >
> > Hello all, When I execute a SQL using SQLDF, how do I get the number
> > of records affected?  I mean, if I run an UPDATE on a data frame, it
> > doesn't tell me if and how many records got updated.  I've read
> > through the documentation and there don't seem to be a way to get this
> > info unless it's done on a database.  Any ideas?
> >
> > Thanks
> > Ravi
> >
> >
> > --
> > This email has been checked for viruses by AVG.
> > https://www.avg.com
> >
> > __
> > 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.
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com
>
>
> --
> This email has been checked for viruses by AVG.
> https://www.avg.com
>


--
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] sqldf and number of records affected

2020-06-11 Thread Gabor Grothendieck
Here is an example.  Ignore the warning or use the workaround discussed here
https://github.com/ggrothendieck/sqldf/issues/40
to avoid the warning.

  library(sqldf)
  sqldf()  # use same connection until next sqldf()
  sqldf(c("pragma count_changes = 1", "update BOD set demand = 99
where Time > 4"))
  sqldf("select * from main.BOD")
  sqldf()


On Thu, Jun 11, 2020 at 9:01 AM Ravi Jeyaraman  wrote:
>
> Hello all, When I execute a SQL using SQLDF, how do I get the number of
> records affected?  I mean, if I run an UPDATE on a data frame, it doesn't
> tell me if and how many records got updated.  I've read through the
> documentation and there don't seem to be a way to get this info unless it's
> done on a database.  Any ideas?
>
> Thanks
> Ravi
>
>
> --
> This email has been checked for viruses by AVG.
> https://www.avg.com
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

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


Re: [R] How to find a split point in a curve?

2020-05-14 Thread Gabor Grothendieck
We can use nls2 to try each value in 10:100 as a possible split point
picking the one with lowest residual sum of squares:

library(nls2)
fm <- nls2(Y ~ cbind(1, pmin(X, X0)), start = data.frame(X0 = 10:100),
  algorithm = "plinear-brute")
plot(Y ~ X)
lines(fitted(fm) ~ X, col = "red")

> fm
Nonlinear regression model
  model: Y ~ cbind(1, pmin(X, X0))
   data: parent.frame()
 X0   .lin1   .lin2
18.  6.5570  0.2616
 residual sum-of-squares: 4.999

Number of iterations to convergence: 91
Achieved convergence tolerance: NA

On Thu, May 14, 2020 at 11:13 AM Luigi Marongiu
 wrote:
>
> Dear all,
> I am trying to find a turning point in some data. In the initial phase, the
> data increases more or less exponentially (thus it is linear in a nat log
> transform), then reaches a plateau. I would like to find the point that
> marks the end of the exponential phase.
> I understand that the function spline can build a curve; is it possible
> with it to find the turning points? I have no idea of how to use spline
> though.
> Here is a working example.
> Thank you
>
> ```
> Y = c(259, 716, 1404, 2173, 3944, 5403, 7140, 9121,
>   11220, 13809, 16634, 19869, 23753, 27447,
>   30590, 33975, 36627, 39600, 42067, 44082,
>   58190, 63280, 65921, 67929, 69977, 71865,
>   73614, 74005, 74894, 75717, 76365, 76579,
>   77087, 77493, 77926, 78253, 78680, 79253,
>   79455, 79580, 79699, 79838, 79981, 80080,
>   80124, 80164, 80183, 80207, 80222, 80230,
>   80241, 80261, 80261, 80277, 80290, 80303,
>   80337, 80376, 80422, 80461, 80539, 80586,
>   80653, 80708, 80762, 80807, 80807, 80886,
>   80922, 80957, 80988, 81007, 81037, 81076,
>   81108, 81108, 81171, 81213, 81259, 81358,
>   81466, 81555, 81601, 81647, 81673, 81998,
>   82025, 82041, 82053, 82064, 82094, 82104,
>   82110, 82122, 82133, 82136, 82142, 82164,
>   82168, 82180, 82181, 82184, 82187, 82188,
>   82190, 82192, 82193, 82194)
> Y = log(Y)
> X = 1:length(Y)
> plot(Y ~ X, ylab = "Ln(Y)", xlim=c(0,10, main="zoomed in"))
> abline(lm(Y[1:3] ~ X[1:3]))
> abline(lm(Y[1:5] ~ X[1:5]), lty=2)
> text(7, 6, "After third or fifth point, there is deviance", pos=3)
> text(2.5, 10, "Solid line: linear model points 1:3", pos =3)
> text(2.5, 9, "Dashed line: linear model points 1:5", pos =3)
> plot(Y ~ X, ylab = "Ln(Y)", xlim=c(0,10, main="overall"))
> abline(lm(Y[1:3] ~ X[1:3]))
> ```
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Relatively Simple Maximization Using Optim Doesnt Optimize

2020-03-14 Thread Gabor Grothendieck
It seems CG is having problems with the cube root.  This converges while
still using CG:

S1 <- optim(1001,function(x) (production1(x)^3), method = "CG",
  control = list(fnscale=-1))

On Thu, Mar 12, 2020 at 9:34 AM Skyler Saleebyan
 wrote:
>
> I am trying to familiarize myself with optim() with a relatively simple
> maximization.
>
> Description:
> L and K are two terms which are constrained to add up to a total 10
> (with respective weights to each). To map this constraint I plugged K into
> the function (to make this as simple as possible.)
>
> Together these two feed into one nonlinear function which is the product of
> two monotonic (on the positive interval) functions. Then that numbers is
> returned in a function fed to optim, which should maximize the output by
> adjusting L. The whole code is:
>
> production1 <- function(L){
>   budget=10
>   Lcost=12
>   Kcost=15
>   K=(budget-L*Lcost)/Kcost
>   machines=0.05*L^(2/3)*K^(1/3)
>   return(machines)
> }
>
> # production1(6000) #example of number with much higher output vs optim
> result
> S1=optim(1001,production1,method="CG",control=list(fnscale=-1))
> S1
>
> Output:
> $par
> [1] 1006.536
>
> $value
> [1] 90.54671
>
> $counts
> function gradient
>  201  101
>
> $convergence
> [1] 1
>
> $message
> NULL
>
>
> For some reason this never explores the problem space and just spits out
> some answer close to the initial condition. What am I doing wrong?
>
> Thanks,
> Skyler S.
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Converting irregular time series data into ts object

2020-02-19 Thread Gabor Grothendieck
Sorry there were some errors in my email. Use this instead.

Assuming that they both cover the same period of time then
if you are willing to throw away some points then
consider using only these 256 elements from the 305 series

  round(seq(1, 305, length = 256))
  ##  [1]   1   2   3   5

That is use the 1st 2nd, 3rd, 5th, etc. point in each year from the
305 series.  This aligns them by throwing away 305-256=49
points per year in the 305 series so that both series can be
set up with a frequency of 256 points per year.

On Wed, Feb 19, 2020 at 10:37 AM Gabor Grothendieck
 wrote:
>
> Assuming that they both cover the same period of time then
> if you are willing to throw away some points then
> consider using only these 256 elements from the 305 series
>
>   round(seq(1, 305, length = 50))
>   ## [1]   1   7  13  20  26  32  38  44 ...etc...
>
> That is use the 1st ,7th, 13th, etc. point in each year from the
> 305 series.  This aligns them by throwing away 305-256=49
> points per year in the 305 series so that both series can be
> set up with a frequency of 256 points per year.
>
>
>
> On Wed, Feb 19, 2020 at 8:36 AM Upananda Pani  wrote:
> >
> > Dear All,
> >
> > I want to convert irregular time series daily data in to ts objects. For
> > some years I have 305 days data and some years I have 256 days.
> >
> > I need your suggestion regarding the same.
> >
> > I have read tutorial on the same but not able to find solutions.
> >
> > With regards,
> > Upananda
> >
> > [[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.
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Converting irregular time series data into ts object

2020-02-19 Thread Gabor Grothendieck
Assuming that they both cover the same period of time then
if you are willing to throw away some points then
consider using only these 256 elements from the 305 series

  round(seq(1, 305, length = 50))
  ## [1]   1   7  13  20  26  32  38  44 ...etc...

That is use the 1st ,7th, 13th, etc. point in each year from the
305 series.  This aligns them by throwing away 305-256=49
points per year in the 305 series so that both series can be
set up with a frequency of 256 points per year.



On Wed, Feb 19, 2020 at 8:36 AM Upananda Pani  wrote:
>
> Dear All,
>
> I want to convert irregular time series daily data in to ts objects. For
> some years I have 305 days data and some years I have 256 days.
>
> I need your suggestion regarding the same.
>
> I have read tutorial on the same but not able to find solutions.
>
> With regards,
> Upananda
>
> [[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.



--
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Strange behaviour of R?

2020-01-20 Thread Gabor Grothendieck


Normally one uses match.fun to avoid such problems.
This will give the error shown even if FUN is defined in the global environment.

  test <- function(FUN, args) {
 FUN <- match.fun(FUN)
 print(FUN)
 FUN(args)
  }
  test(NULL, 1:10)
  ## Error in match.fun(FUN) : 'NULL' is not a function, character or symbol

On Fri, Jan 17, 2020 at 5:41 AM Duncan Murdoch  wrote:
>
> On 17/01/2020 2:33 a.m., Sigbert Klinke wrote:
> > Hi,
> >
> > I wrote a function like
> >
> > test <- function(FUN, args) {
> > print(FUN)
> > FUN(args)
> > }
> >
> > When I call it lieke this
> >
> > test(mean, 1:10)
> > test(NULL, 1:10)
> >
> > then the second call still uses mean, although I set FUN to NULL. Is
> > that ok?
>
> You probably have a function defined in your global environment that is
> named FUN and acts like mean.
>
> The general rule in R is that it only looks for objects of mode function
> when trying to find something used as a function.  So in your second
> case, when trying to evaluate FUN(args), R will look for a function
> named FUN in the local evaluation frame, and won't find one:  FUN is
> NULL there.  Then it will go to the environment of test, which is likely
> the global environment, and look there.  That's where it probably found
> the function.
>
> For example, try this:
>
> FUN <- function(...) print('FUN was called')
>
>
> test <- function(FUN, args) {
> print(FUN)
> FUN(args)
> }
>
> test(NULL, 1:10)
>
> Duncan Murdoch
>
> >
> > Actually, I used something like
> >
> > test(mean, list(x=1:10, na.rm=TRUE))
> >
> > which actually crashed R, but I can not reproduce it. Of course, when I
> > replaced FUN(args) with do.call(FUN, args) then everything works fine.
> >
> > Sigbert
> >
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] outer join of xts's

2020-01-02 Thread Gabor Grothendieck
join = "left" only applies with merge.xts if there are two objects.
If there are more it acts the same as join = TRUE..
See the Details section of ?merge.xts

On Thu, Jan 2, 2020 at 1:29 PM Duncan Murdoch  wrote:
>
> On 02/01/2020 9:31 a.m., Eric Berger wrote:
> > Hi Gabor,
> > This is great, thanks. It brought the time down to about 4 seconds.
> > The command
> > do.call("merge.xts",L)
> > also works in this case.
> > Suppose that instead of the default "outer" join I wanted to use, say, a
> > "left" join.
> > Is that possible? I tried a few ways of adding the
> > join="left"
> > parameter to the do.call() command but I could not get the syntax to work
> > (assuming it's even possible).
>
> This should work:
>
>do.call("merge", c(L, join = "left"))
>
> The second argument to do.call is a list which becomes the arguments to
> the function being called.  Your time series should be unnamed entries
> in the list, while other arguments to merge() should be named.
>
> Duncan Murdoch
>
> >
> > Thanks,
> > Eric
> >
> >
> > On Thu, Jan 2, 2020 at 3:23 PM Gabor Grothendieck 
> > wrote:
> >
> >> You don't need Reduce as xts already supports mutliway merges.  This
> >> perfroms one
> >> multiway merge rather than  k-1 two way merges.
> >>
> >>  do.call("merge", L)
> >>
> >> On Thu, Jan 2, 2020 at 6:13 AM Eric Berger  wrote:
> >>>
> >>> Hi,
> >>> I have a list L of about 2,600 xts's.
> >>> Each xts has a single numeric column. About 90% of the xts's have
> >>> approximately 500 rows, and the rest have fewer than 500 rows.
> >>> I create a single xts using the command
> >>>
> >>> myXts <- Reduce( merge.xts, L )
> >>>
> >>> By default, merge.xts() does an outer join (which is what I want).
> >>>
> >>> The command takes about 80 seconds to complete.
> >>> I have plenty of RAM on my computer.
> >>>
> >>> Are there faster ways to accomplish this task?
> >>>
> >>> Thanks,
> >>> Eric
> >>>
> >>>  [[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.
> >>
> >>
> >>
> >> --
> >> Statistics & Software Consulting
> >> GKX Group, GKX Associates Inc.
> >> tel: 1-877-GKX-GROUP
> >> email: ggrothendieck at gmail.com
> >>
> >
> >   [[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.
> >
>


-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] outer join of xts's

2020-01-02 Thread Gabor Grothendieck
> It is not clear what multiway left join means but merge.zoo (though
> not merge.xts) supports a generalized all= argument which is a logical
> vector having the same length as L that can be TRUE or FALSE for each
> object merged.  The objects corresponding to TRUE will have all their
> times included in the result but the ones with FALSE will only be
> included if they correspond to an already existing time. merge.zoo is
> R based whereas merge.xts is C based so I would not expect it to be as
> fast although it is more powerful.
>
> If All is the logical vector having the same length as L then:
>
> Lzoo <- lapply(L, as.zoo)
> do.call("merge",  c(Lzoo, list(all = All))
>
> On Thu, Jan 2, 2020 at 9:31 AM Eric Berger  wrote:
> >
> > Hi Gabor,
> > This is great, thanks. It brought the time down to about 4 seconds.
> > The command
> > do.call("merge.xts",L)
> > also works in this case.
> > Suppose that instead of the default "outer" join I wanted to use, say, a 
> > "left" join.
> > Is that possible? I tried a few ways of adding the
> > join="left"
> > parameter to the do.call() command but I could not get the syntax to work 
> > (assuming it's even possible).
> >
> > Thanks,
> > Eric
> >
> >
> > On Thu, Jan 2, 2020 at 3:23 PM Gabor Grothendieck  
> > wrote:
> >>
> >> You don't need Reduce as xts already supports mutliway merges.  This
> >> perfroms one
> >> multiway merge rather than  k-1 two way merges.
> >>
> >> do.call("merge", L)
> >>
> >> On Thu, Jan 2, 2020 at 6:13 AM Eric Berger  wrote:
> >> >
> >> > Hi,
> >> > I have a list L of about 2,600 xts's.
> >> > Each xts has a single numeric column. About 90% of the xts's have
> >> > approximately 500 rows, and the rest have fewer than 500 rows.
> >> > I create a single xts using the command
> >> >
> >> > myXts <- Reduce( merge.xts, L )
> >> >
> >> > By default, merge.xts() does an outer join (which is what I want).
> >> >
> >> > The command takes about 80 seconds to complete.
> >> > I have plenty of RAM on my computer.
> >> >
> >> > Are there faster ways to accomplish this task?
> >> >
> >> > Thanks,
> >> > Eric
> >> >
> >> > [[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.
> >>
> >>
> >>
> >> --
> >> Statistics & Software Consulting
> >> GKX Group, GKX Associates Inc.
> >> tel: 1-877-GKX-GROUP
> >> email: ggrothendieck at gmail.com
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] outer join of xts's

2020-01-02 Thread Gabor Grothendieck
It is not clear what multiway left join means but merge.zoo (though
not merge.xts) supports a generalized all= argument which is a logical
vector having the same length as L that can be TRUE or FALSE for each
object merged.  The objects corresponding to TRUE will have all their
times included in the result but the ones with FALSE will only be
included if they correspond to an already existing time. merge.zoo is
R based whereas merge.xts is C based so I would not expect it to be as
fast although it is more powerful.

If All is the logical vector having the same length as L then:

Lzoo <- lapply(L, as.zoo)
do.call("merge",  c(Lzoo, list(all = All))

On Thu, Jan 2, 2020 at 9:31 AM Eric Berger  wrote:
>
> Hi Gabor,
> This is great, thanks. It brought the time down to about 4 seconds.
> The command
> do.call("merge.xts",L)
> also works in this case.
> Suppose that instead of the default "outer" join I wanted to use, say, a 
> "left" join.
> Is that possible? I tried a few ways of adding the
> join="left"
> parameter to the do.call() command but I could not get the syntax to work 
> (assuming it's even possible).
>
> Thanks,
> Eric
>
>
> On Thu, Jan 2, 2020 at 3:23 PM Gabor Grothendieck  
> wrote:
>>
>> You don't need Reduce as xts already supports mutliway merges.  This
>> perfroms one
>> multiway merge rather than  k-1 two way merges.
>>
>> do.call("merge", L)
>>
>> On Thu, Jan 2, 2020 at 6:13 AM Eric Berger  wrote:
>> >
>> > Hi,
>> > I have a list L of about 2,600 xts's.
>> > Each xts has a single numeric column. About 90% of the xts's have
>> > approximately 500 rows, and the rest have fewer than 500 rows.
>> > I create a single xts using the command
>> >
>> > myXts <- Reduce( merge.xts, L )
>> >
>> > By default, merge.xts() does an outer join (which is what I want).
>> >
>> > The command takes about 80 seconds to complete.
>> > I have plenty of RAM on my computer.
>> >
>> > Are there faster ways to accomplish this task?
>> >
>> > Thanks,
>> > Eric
>> >
>> > [[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.
>>
>>
>>
>> --
>> Statistics & Software Consulting
>> GKX Group, GKX Associates Inc.
>> tel: 1-877-GKX-GROUP
>> email: ggrothendieck at gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] outer join of xts's

2020-01-02 Thread Gabor Grothendieck
You don't need Reduce as xts already supports mutliway merges.  This
perfroms one
multiway merge rather than  k-1 two way merges.

do.call("merge", L)

On Thu, Jan 2, 2020 at 6:13 AM Eric Berger  wrote:
>
> Hi,
> I have a list L of about 2,600 xts's.
> Each xts has a single numeric column. About 90% of the xts's have
> approximately 500 rows, and the rest have fewer than 500 rows.
> I create a single xts using the command
>
> myXts <- Reduce( merge.xts, L )
>
> By default, merge.xts() does an outer join (which is what I want).
>
> The command takes about 80 seconds to complete.
> I have plenty of RAM on my computer.
>
> Are there faster ways to accomplish this task?
>
> Thanks,
> Eric
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] giving priority to stats package

2019-11-25 Thread Gabor Grothendieck
Goold idea.  This seems to work.

  library(dplyr, pos = grep("package:stats", search()) + 1)

On Mon, Nov 25, 2019 at 1:27 PM Greg Snow <538...@gmail.com> wrote:
>
> You could use the `pos` arg to place the newly loaded package(s) on
> the search path after the stats package.  That would give priority for
> any functions in the stats package over the newly loaded package (but
> also give priority for any other packages earlier on the search path).
>
> On Sat, Nov 23, 2019 at 6:25 AM Gabor Grothendieck
>  wrote:
> >
> > library and require have new args in 3.6 giving additional control
> > over conflicts.  This seems very useful but I was wondering if there
> > were some, preferabley simple, way to give existing loaded packages
> > priority without knowing the actual conflicts in advance.  For example
> >
> > library(dplyr, exclude = c("filter", "lag"))
> >
> > works to avoid masking those names in stats that would otherwise be
> > masked by that package but I had to know in advance that filter and
> > lag were the conflicting names.
> >
> > --
> > Statistics & Software Consulting
> > GKX Group, GKX Associates Inc.
> > tel: 1-877-GKX-GROUP
> > email: ggrothendieck at gmail.com
> >
> > __
> > 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.
>
>
>
> --
> Gregory (Greg) L. Snow Ph.D.
> 538...@gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] giving priority to stats package

2019-11-23 Thread Gabor Grothendieck
library and require have new args in 3.6 giving additional control
over conflicts.  This seems very useful but I was wondering if there
were some, preferabley simple, way to give existing loaded packages
priority without knowing the actual conflicts in advance.  For example

library(dplyr, exclude = c("filter", "lag"))

works to avoid masking those names in stats that would otherwise be
masked by that package but I had to know in advance that filter and
lag were the conflicting names.

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] [r] Applying rollapply in VR test

2019-08-27 Thread Gabor Grothendieck
If I run this (which is identical to your code except I supplied some
random input since your
post did not include any) I get no error:

  library(vrtest)
  library(zoo)
  set.seed(123)
  x <- rnorm(100)
  y <- zoo(x)
  z <- rollapply(y,50, function(x) AutoBoot.test(x,nboot=30,
wild="Normal")$AutoBoot.test)


On Tue, Aug 27, 2019 at 1:20 AM Subhamitra Patra
 wrote:
>
> Dear R users,
>
> I want to use rollapply function from the zoo package with the VR test by
> using the following code.
>
> library(vrtest)
> library(zoo)
> x <- Data
> y <- zoo(x)
> z <- rollapply(y,50, function(x) AutoBoot.test(x,nboot=30,
> wild="Normal")$AutoBoot.test)
>
>
> But, I am getting some error message, and unable to do that. Is there any
> mistake in my code?
>
> Please help me with this.
>
> Thank you.
>
> --
> *Best Regards,*
> *Subhamitra Patra*
> *Phd. Research Scholar*
> *Department of Humanities and Social Sciences*
> *Indian Institute of Technology, Kharagpur*
> *INDIA*
>
> [image: Mailtrack]
> 
> Sender
> notified by
> Mailtrack
> 
> 08/27/19,
> 10:48:23 AM
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Multiple Lags with Dplyr

2019-04-23 Thread Gabor Grothendieck
lag.zoo supports vector-based lags on zoo objects.
A few caveats:

- dplyr's lag clobbers the base R lag (which you need to
invoke lag's methods) so if you have dplyr loaded be sure
to refer to stats::lag.

- dplyr's lag works backwards relative to the standard set
in base R so dplyr::lag(x, 1) corresponds to stat::lag(x, -1) in
base R

- zoo follows base R's standard

- you can use as.data.frame or fortify.zoo to convert a zoo
object to a data frame if you need that.  The first one
drops the time index and the second one includes it.

  library(zoo)
  stats::lag(zoo(d2$x1), 0:-2)

giving this zoo object:

   lag0 lag-1 lag-2
1 1NANA
2 2 1NA
3 3 2 1
4 4 3 2
5 5 4 3
6 6 5 4
7 7 6 5
8 8 7 6
9 9 8 7
10   10 9 8


On Tue, Apr 23, 2019 at 9:10 AM Lorenzo Isella  wrote:
>
> Dear All,
> I refer to the excellent post at
>
> https://purrple.cat/blog/2018/03/02/multiple-lags-with-tidy-evaluation/
>
> What I want to do is to create a function capable, ą la dplyr, to
> generate new columns which are a lagged version of existing columns in
> a data frame.
> For instance, you can do this manually as
>
>
> d2 <- tibble(x1 =1:10, x2=10:19,  x3=50:59)
>
>
> d3 <- d2%>%mutate(x1lag1=lag(x1, 1), x1lag2=lag(x1,2))
>
>
> but this becomes quickly tedious when you need to take several lags of
> different columns.
> One solution in the link above is the following
>
>
> lags <- function(var, n=10){
>   var <- enquo(var)
>
>   indices <- seq_len(n)
>   map( indices, ~quo(lag(!!var, !!.x)) ) %>%
> set_names(sprintf("lag_%s_%02d", quo_text(var), indices))
>
> }
>
>
> d4 <- d2 %>%
>   mutate( !!!lags(x1, 3), !!!lags(x2,3) )
>
>
> does anybody know how this could be made more general? I mean that I
> would like to take a fixed number of lags of a list of columns (x1 and
> x2, for instance), just by passing the list of columns and without
> repeating the commands for x1 and x2.
> Any suggestion is appreciated.
> Cheers
>
> Lorenzo
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Newbie Question on R versus Matlab/Octave versus C

2019-01-29 Thread Gabor Grothendieck
Two additional comments:

- depending on the nature of your problem you may be able to get an
analytic solution using branching processes. I found this approach
successful when I once had to model stem cell growth.

- in addition to NetLogo another alternative to R would be the Julia
language which is motivated to some degree by Octave but is actually
quite different and is particularly suitable in terms of performance
for iterative computations where one iteration depends on the prior
one.

On Mon, Jan 28, 2019 at 6:32 PM Gabor Grothendieck
 wrote:
>
> R has many similarities to Octave.  Have a look at:
>
> https://cran.r-project.org/doc/contrib/R-and-octave.txt
> https://CRAN.R-project.org/package=matconv
>
> On Mon, Jan 28, 2019 at 4:58 PM Alan Feuerbacher  wrote:
> >
> > Hi,
> >
> > I recently learned of the existence of R through a physicist friend who
> > uses it in his research. I've used Octave for a decade, and C for 35
> > years, but would like to learn R. These all have advantages and
> > disadvantages for certain tasks, but as I'm new to R I hardly know how
> > to evaluate them. Any suggestions?
> >
> > Thanks!
> >
> > ---
> > This email has been checked for viruses by Avast antivirus software.
> > https://www.avast.com/antivirus
> >
> > __
> > 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.
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] [FORGED] Newbie Question on R versus Matlab/Octave versus C

2019-01-28 Thread Gabor Grothendieck
This would be a suitable application for NetLogo.  The R package
RNetLogo provides an interface.  In a few lines of code you get a
simulation with graphics.

On Mon, Jan 28, 2019 at 7:00 PM Alan Feuerbacher  wrote:
>
> On 1/28/2019 4:20 PM, Rolf Turner wrote:
> >
> > On 1/29/19 10:05 AM, Alan Feuerbacher wrote:
> >
> >> Hi,
> >>
> >> I recently learned of the existence of R through a physicist friend
> >> who uses it in his research. I've used Octave for a decade, and C for
> >> 35 years, but would like to learn R. These all have advantages and
> >> disadvantages for certain tasks, but as I'm new to R I hardly know how
> >> to evaluate them. Any suggestions?
> >
> > * C is fast, but with a syntax that is (to my mind) virtually
> >incomprehensible.  (You probably think differently about this.)
>
> I've been doing it long enough that I have little problem with it,
> except for pointers. :-)
>
> > * In C, you essentially have to roll your own for all tasks; in R,
> >practically anything (well ...) that you want to do has already
> >been programmed up.  CRAN is a wonderful resource, and there's more
> >on github.
>  >
> > * The syntax of R meshes beautifully with *my* thought patterns; YMMV.
> >
> > * Why not just bog in and try R out?  It's free, it's readily available,
> >and there are a number of good online tutorials.
>
> I just installed R on my Linux Fedora system, so I'll do that.
>
> I wonder if you'd care to comment on my little project that prompted
> this? As part of another project, I wanted to model population growth
> starting from a handful of starting individuals. This is exponential in
> the long run, of course, but I wanted to see how a few basic parameters
> affected the outcome. Using Octave, I modeled a single person as a
> "cell", which in Octave has a good deal of overhead. The program
> basically looped over the entire population, and updated each person
> according to the parameters, which included random statistical
> variations. So when the total population reached, say 10,000, and an
> update time of 1 day, the program had to execute 10,000 x 365 update
> operations for each year of growth. For large populations, say 100,000,
> the program did not return even after 24 hours of run time.
>
> So I switched to C, and used its "struct" declaration and an array of
> structs to model the population. This allowed the program to complete in
> under a minute as opposed to 24 hours+. So in line with your comments, C
> is far more efficient than Octave.
>
> How do you think R would fare in this simulation?
>
> Alan
>
>
> ---
> This email has been checked for viruses by Avast antivirus software.
> https://www.avast.com/antivirus
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Newbie Question on R versus Matlab/Octave versus C

2019-01-28 Thread Gabor Grothendieck
R has many similarities to Octave.  Have a look at:

https://cran.r-project.org/doc/contrib/R-and-octave.txt
https://CRAN.R-project.org/package=matconv

On Mon, Jan 28, 2019 at 4:58 PM Alan Feuerbacher  wrote:
>
> Hi,
>
> I recently learned of the existence of R through a physicist friend who
> uses it in his research. I've used Octave for a decade, and C for 35
> years, but would like to learn R. These all have advantages and
> disadvantages for certain tasks, but as I'm new to R I hardly know how
> to evaluate them. Any suggestions?
>
> Thanks!
>
> ---
> This email has been checked for viruses by Avast antivirus software.
> https://www.avast.com/antivirus
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] R code for if-then-do code blocks

2018-12-17 Thread Gabor Grothendieck
There is some discussion of approaches to this here:

https://stackoverflow.com/questions/34096162/dplyr-mutate-replace-on-a-subset-of-rows/34096575#34096575


On Mon, Dec 17, 2018 at 10:30 AM Paul Miller via R-help
 wrote:
>
> Hello All,
>
> Season's greetings!
>
>  Am trying to replicate some SAS code in R. The SAS code uses if-then-do code 
> blocks. I've been trying to do likewise in R as that seems to be the most 
> reliable way to get the same result.
>
> Below is some toy data and some code that does work. There are some things I 
> don't necessarily like about the code though. So I was hoping some people 
> could help make it better. One thing I don't like is that the within function 
> reverses the order of the computed columns such that test1:test5 becomes 
> test5:test1. I've used a mutate to overcome that but would prefer not to have 
> to do so.
>
>  Another, perhaps very small thing, is the need to calculate an ID variable 
> that becomes the basis for a grouping.
>
> I did considerable Internet searching for R code that conditionally computes 
> blocks of code. I didn't find much though and so am wondering if my search 
> terms were not sufficient or if there is some other reason. It occurred to me 
> that maybe if-then-do code blocks like we often see in SAS as are frowned 
> upon and therefore not much implemented.
>
> I'd be interested in seeing more R-compatible approaches if this is the case. 
> I've learned that it's a mistake to try and make R be like SAS. It's better 
> to let R be R. Trouble is I'm not always sure how to do that.
>
> Thanks,
>
> Paul
>
>
> d1 <- data.frame(workshop=rep(1:2,4),
> gender=rep(c("f","m"),each=4))
>
> library(tibble)
> library(plyr)
>
> d2 <- d1 %>%
>   rownames_to_column("ID") %>%
>   mutate(test1 = NA, test2 = NA, test4 = NA, test5 = NA) %>%
>   ddply("ID",
> within,
> if (gender == "f" & workshop == 1) {
>   test1 <- 1
>   test1 <- 6 + test1
>   test2 <- 2 + test1
>   test4 <- 1
>   test5 <- 1
> } else {
>   test1 <- test2 <- test4 <- test5 <- 0
> })
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] year and week to date - before 1/1 and after 12/31

2018-10-18 Thread Gabor Grothendieck
Replace the week in the date with week 2, say -- a week in which
nothing will go wrong and then add or subtract the appropriate number
of weeks.

  d <- c('2016 00 Sun', '2017 53 Sun', '2017 53 Mon') # test data

  as.Date(sub(" .. ", "02", d), "%Y %U %a") +
 7 * (as.numeric(sub(" (..) ...", "\\1", d)) - 2)
  ## [1] "2015-12-27" "2017-12-31" "2018-01-01"

On Tue, Oct 16, 2018 at 10:23 AM peter salzman
 wrote:
>
> hi,
>
> to turn year and week into the date one can do the following:
>
> as.Date('2018 05 Sun', "%Y %W %a")
>
> however, when we want the Sunday (1st day of week) of the 1st week of 2018
> we get NA because 1/1/2018 was on Monday
>
> as.Date('2018 00 Mon',format="%Y %U %a")
> ## 2018-01-01
> as.Date('2018 00 Sun',format="%Y %U %a")
> ## NA
>
> btw the same goes for last week
> as.Date('2017 53 Sun',format="%Y %U %a")
> ## 2017-12-31
> as.Date('2017 53 Mon',format="%Y %U %a")
> ## NA
>
> So my question is :
> how do i get
> from "2018 00 Sun" to 2018-12-31
> and
> from "2017 53 Mon" to 2018-01-01
>
> i realize i can loop over days of week and do some if/then statements,
> but is there a built in function?
>
> thank you
> peter
>
>
>
>
>
> --
> Peter Salzman, PhD
> Department of Biostatistics and Computational Biology
> University of Rochester
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] looking for formula parser that allows coefficients

2018-08-24 Thread Gabor Grothendieck
The isChar function used in Parse is:

  isChar <- function(e, ch) identical(e, as.symbol(ch))
On Fri, Aug 24, 2018 at 10:06 PM Gabor Grothendieck
 wrote:
>
> Also here is a solution that uses formula processing rather than
> string processing.
> No packages are used.
>
> Parse <- function(e) {
>   if (length(e) == 1) {
> if (is.numeric(e)) return(e)
> else setNames(1, as.character(e))
>   } else {
> if (isChar(e[[1]], "*")) {
>x1 <- Recall(e[[2]])
>x2 <- Recall(e[[3]])
>setNames(unname(x1 * x2), paste0(names(x1), names(x2)))
> } else if (isChar(e[[1]], "+")) c(Recall(e[[2]]), Recall(e[[3]]))
> else if (isChar(e[[1]], "-")) {
>   if (length(e) == 2) -1 * Recall(e[[2]])
>   else c(Recall(e[[2]]), -Recall(e[[3]]))
> } else if (isChar(e[[1]], ":")) setNames(1, paste(e[-1], collapse = ":"))
>   }
> }
>
> # test
> fo <- y ~ 2 - 1.1 * x1 + x3 - x1:x3 + 0.2 * x2:x2
> Parse(fo[[3]])
>
> giving:
>
>  x1x3 x1:x3 x2:x2
>   2.0  -1.1   1.0  -1.0   0.2
> On Wed, Aug 22, 2018 at 11:50 AM Paul Johnson  wrote:
> >
> > Thanks as usual.  I owe you more KU decorations soon.
> > On Wed, Aug 22, 2018 at 2:34 AM Gabor Grothendieck
> >  wrote:
> > >
> > > Some string manipulation can convert the formula to a named vector such as
> > > the one shown at the end of your post.
> > >
> > > library(gsubfn)
> > >
> > > # input
> > > fo <- y ~ 2 - 1.1 * x1 + x3 - x1:x3 + 0.2 * x2:x2
> > >
> > > pat <- "([+-])? *(\\d\\S*)? *\\*? *([[:alpha:]]\\S*)?"
> > > ch <- format(fo[[3]])
> > > m <- matrix(strapplyc(ch, pat)[[1]], 3)
> > > m <- m[, colSums(m != "") > 0]
> > > m[2, m[2, ] == ""] <- 1
> > > m[3, m[3, ] == ""] <- "(Intercept)"
> > > co <- as.numeric(paste0(m[1, ], m[2, ]))
> > > v <- m[3, ]
> > > setNames(co, v)
> > > ## (Intercept)  x1  x3   x1:x3   x2:x2
> > > ## 2.0-1.1 1.0-1.0 0.2
> > > On Tue, Aug 21, 2018 at 6:46 PM Paul Johnson  wrote:
> > > >
> > > > Can you point me at any packages that allow users to write a
> > > > formula with coefficients?
> > > >
> > > > I want to write a data simulator that has a matrix X with lots
> > > > of columns, and then users can generate predictive models
> > > > by entering a formula that uses some of the variables, allowing
> > > > interactions, like
> > > >
> > > > y ~ 2 + 1.1 * x1 + 3 * x3 + 0.1 * x1:x3 + 0.2 * x2:x2
> > > >
> > > > Currently, in the rockchalk package, I have a function simulates
> > > > data (genCorrelatedData2), but my interface to enter the beta
> > > > coefficients is poor.  I assumed user would always enter 0's as
> > > > place holder for the unused coefficients, and the intercept is
> > > > always first. The unnamed vector is too confusing.  I have them specify:
> > > >
> > > > c(2, 1.1, 0, 3, 0, 0, 0.2, ...)
> > > >
> > > > I the documentation I say (ridiculously) it is easy to figure out from
> > > > the examples, but it really isnt.
> > > > It function prints out the equation it thinks you intended, thats
> > > > minimum protection against user error, but still not very good:
> > > >
> > > > dat <- genCorrelatedData2(N = 10, rho = 0.0,
> > > >   beta = c(1, 2, 1, 1, 0, 0.2, 0, 0, 0),
> > > >   means = c(0,0,0), sds = c(1,1,1), stde = 0)
> > > > [1] "The equation that was calculated was"
> > > > y = 1 + 2*x1 + 1*x2 + 1*x3
> > > >  + 0*x1*x1 + 0.2*x2*x1 + 0*x3*x1
> > > >  + 0*x1*x2 + 0*x2*x2 + 0*x3*x2
> > > >  + 0*x1*x3 + 0*x2*x3 + 0*x3*x3
> > > >  + N(0,0) random error
> > > >
> > > > But still, it is not very good.
> > > >
> > > > As I look at this now, I realize expect just the vech, not the whole 
> > > > vector
> > > > of all interaction terms, so it is even more difficult than I thought 
> > > > to get the
> > > > correct input.Hence, I'd like to let the user write a formula.
> > > >
> > > > The alternative for the user interface is to have named coefficients.
> > > > I can more or less easily allow a named vector for beta
> > > >

Re: [R] looking for formula parser that allows coefficients

2018-08-24 Thread Gabor Grothendieck
Also here is a solution that uses formula processing rather than
string processing.
No packages are used.

Parse <- function(e) {
  if (length(e) == 1) {
if (is.numeric(e)) return(e)
else setNames(1, as.character(e))
  } else {
if (isChar(e[[1]], "*")) {
   x1 <- Recall(e[[2]])
   x2 <- Recall(e[[3]])
   setNames(unname(x1 * x2), paste0(names(x1), names(x2)))
} else if (isChar(e[[1]], "+")) c(Recall(e[[2]]), Recall(e[[3]]))
else if (isChar(e[[1]], "-")) {
  if (length(e) == 2) -1 * Recall(e[[2]])
  else c(Recall(e[[2]]), -Recall(e[[3]]))
} else if (isChar(e[[1]], ":")) setNames(1, paste(e[-1], collapse = ":"))
  }
}

# test
fo <- y ~ 2 - 1.1 * x1 + x3 - x1:x3 + 0.2 * x2:x2
Parse(fo[[3]])

giving:

 x1x3 x1:x3 x2:x2
  2.0  -1.1   1.0  -1.0   0.2
On Wed, Aug 22, 2018 at 11:50 AM Paul Johnson  wrote:
>
> Thanks as usual.  I owe you more KU decorations soon.
> On Wed, Aug 22, 2018 at 2:34 AM Gabor Grothendieck
>  wrote:
> >
> > Some string manipulation can convert the formula to a named vector such as
> > the one shown at the end of your post.
> >
> > library(gsubfn)
> >
> > # input
> > fo <- y ~ 2 - 1.1 * x1 + x3 - x1:x3 + 0.2 * x2:x2
> >
> > pat <- "([+-])? *(\\d\\S*)? *\\*? *([[:alpha:]]\\S*)?"
> > ch <- format(fo[[3]])
> > m <- matrix(strapplyc(ch, pat)[[1]], 3)
> > m <- m[, colSums(m != "") > 0]
> > m[2, m[2, ] == ""] <- 1
> > m[3, m[3, ] == ""] <- "(Intercept)"
> > co <- as.numeric(paste0(m[1, ], m[2, ]))
> > v <- m[3, ]
> > setNames(co, v)
> > ## (Intercept)  x1  x3   x1:x3   x2:x2
> > ## 2.0-1.1 1.0-1.0 0.2
> > On Tue, Aug 21, 2018 at 6:46 PM Paul Johnson  wrote:
> > >
> > > Can you point me at any packages that allow users to write a
> > > formula with coefficients?
> > >
> > > I want to write a data simulator that has a matrix X with lots
> > > of columns, and then users can generate predictive models
> > > by entering a formula that uses some of the variables, allowing
> > > interactions, like
> > >
> > > y ~ 2 + 1.1 * x1 + 3 * x3 + 0.1 * x1:x3 + 0.2 * x2:x2
> > >
> > > Currently, in the rockchalk package, I have a function simulates
> > > data (genCorrelatedData2), but my interface to enter the beta
> > > coefficients is poor.  I assumed user would always enter 0's as
> > > place holder for the unused coefficients, and the intercept is
> > > always first. The unnamed vector is too confusing.  I have them specify:
> > >
> > > c(2, 1.1, 0, 3, 0, 0, 0.2, ...)
> > >
> > > I the documentation I say (ridiculously) it is easy to figure out from
> > > the examples, but it really isnt.
> > > It function prints out the equation it thinks you intended, thats
> > > minimum protection against user error, but still not very good:
> > >
> > > dat <- genCorrelatedData2(N = 10, rho = 0.0,
> > >   beta = c(1, 2, 1, 1, 0, 0.2, 0, 0, 0),
> > >   means = c(0,0,0), sds = c(1,1,1), stde = 0)
> > > [1] "The equation that was calculated was"
> > > y = 1 + 2*x1 + 1*x2 + 1*x3
> > >  + 0*x1*x1 + 0.2*x2*x1 + 0*x3*x1
> > >  + 0*x1*x2 + 0*x2*x2 + 0*x3*x2
> > >  + 0*x1*x3 + 0*x2*x3 + 0*x3*x3
> > >  + N(0,0) random error
> > >
> > > But still, it is not very good.
> > >
> > > As I look at this now, I realize expect just the vech, not the whole 
> > > vector
> > > of all interaction terms, so it is even more difficult than I thought to 
> > > get the
> > > correct input.Hence, I'd like to let the user write a formula.
> > >
> > > The alternative for the user interface is to have named coefficients.
> > > I can more or less easily allow a named vector for beta
> > >
> > > beta = c("(Intercept)" = 1, "x1" = 2, "x2" = 1, "x3" = 1, "x2:x1" = 0.1)
> > >
> > > I could build a formula from that.  That's not too bad. But I still think
> > > it would be cool to allow formula input.
> > >
> > > Have you ever seen it done?
> > > pj
> > > --
> > > Paul E. Johnson   http://pj.freefaculty.org
> > > Director, Center for Research Methods and Data Analysis 
> > > http://crmda.ku.edu
> > >
> > > To write to me directly, pleas

Re: [R] looking for formula parser that allows coefficients

2018-08-22 Thread Gabor Grothendieck
Some string manipulation can convert the formula to a named vector such as
the one shown at the end of your post.

library(gsubfn)

# input
fo <- y ~ 2 - 1.1 * x1 + x3 - x1:x3 + 0.2 * x2:x2

pat <- "([+-])? *(\\d\\S*)? *\\*? *([[:alpha:]]\\S*)?"
ch <- format(fo[[3]])
m <- matrix(strapplyc(ch, pat)[[1]], 3)
m <- m[, colSums(m != "") > 0]
m[2, m[2, ] == ""] <- 1
m[3, m[3, ] == ""] <- "(Intercept)"
co <- as.numeric(paste0(m[1, ], m[2, ]))
v <- m[3, ]
setNames(co, v)
## (Intercept)  x1  x3   x1:x3   x2:x2
## 2.0-1.1 1.0-1.0 0.2
On Tue, Aug 21, 2018 at 6:46 PM Paul Johnson  wrote:
>
> Can you point me at any packages that allow users to write a
> formula with coefficients?
>
> I want to write a data simulator that has a matrix X with lots
> of columns, and then users can generate predictive models
> by entering a formula that uses some of the variables, allowing
> interactions, like
>
> y ~ 2 + 1.1 * x1 + 3 * x3 + 0.1 * x1:x3 + 0.2 * x2:x2
>
> Currently, in the rockchalk package, I have a function simulates
> data (genCorrelatedData2), but my interface to enter the beta
> coefficients is poor.  I assumed user would always enter 0's as
> place holder for the unused coefficients, and the intercept is
> always first. The unnamed vector is too confusing.  I have them specify:
>
> c(2, 1.1, 0, 3, 0, 0, 0.2, ...)
>
> I the documentation I say (ridiculously) it is easy to figure out from
> the examples, but it really isnt.
> It function prints out the equation it thinks you intended, thats
> minimum protection against user error, but still not very good:
>
> dat <- genCorrelatedData2(N = 10, rho = 0.0,
>   beta = c(1, 2, 1, 1, 0, 0.2, 0, 0, 0),
>   means = c(0,0,0), sds = c(1,1,1), stde = 0)
> [1] "The equation that was calculated was"
> y = 1 + 2*x1 + 1*x2 + 1*x3
>  + 0*x1*x1 + 0.2*x2*x1 + 0*x3*x1
>  + 0*x1*x2 + 0*x2*x2 + 0*x3*x2
>  + 0*x1*x3 + 0*x2*x3 + 0*x3*x3
>  + N(0,0) random error
>
> But still, it is not very good.
>
> As I look at this now, I realize expect just the vech, not the whole vector
> of all interaction terms, so it is even more difficult than I thought to get 
> the
> correct input.Hence, I'd like to let the user write a formula.
>
> The alternative for the user interface is to have named coefficients.
> I can more or less easily allow a named vector for beta
>
> beta = c("(Intercept)" = 1, "x1" = 2, "x2" = 1, "x3" = 1, "x2:x1" = 0.1)
>
> I could build a formula from that.  That's not too bad. But I still think
> it would be cool to allow formula input.
>
> Have you ever seen it done?
> pj
> --
> Paul E. Johnson   http://pj.freefaculty.org
> Director, Center for Research Methods and Data Analysis http://crmda.ku.edu
>
> To write to me directly, please address me at pauljohn at ku.edu.
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] plotmath and logical operators?

2018-08-21 Thread Gabor Grothendieck
Try this:

plot(1)
tmp <- x >= 3 ~ "&" ~ y <= 3
mtext(tmp)
On Mon, Aug 20, 2018 at 5:00 PM MacQueen, Don via R-help
 wrote:
>
> I would like to use plotmath to annotate a plot with an expression that 
> includes a logical operator.
>
> ## works well
> tmp <- expression(x >= 3)
> plot(1)
> mtext(tmp)
>
> ## not so well
> tmp <- expression(x >= 3 &  y <= 3)
> plot(1)
> mtext(tmp)
>
> Although the text that's displayed makes sense, it won't be obvious to my 
> non-mathematical audience.
>
> I'd appreciate suggestions.
>
>
> I've found a work-around that gets the annotation to look right
>   tmpw <- expression(paste( x >= 3, " & ", y <= 3) )
>   plot(1)
>   mtext(tmpw)
>
>
> But it breaks my original purpose, illustrated by this example:
>
> df <- data.frame(x=1:5, y=1:5)
> tmp <- expression(x >= 3 & y <= 3)
> tmpw <- expression(paste( x >= 3, " & ", y <= 3) )
> with(df, eval(tmp))
> [1] FALSE FALSE  TRUE FALSE FALSE
> with(df, eval(tmpw))
> [1] "FALSE  &  TRUE" "FALSE  &  TRUE" "TRUE  &  TRUE"  "TRUE  &  FALSE" "TRUE 
>  &  FALSE"
>
> Thanks
> -Don
>
> --
> Don MacQueen
> Lawrence Livermore National Laboratory
> 7000 East Ave., L-627
> Livermore, CA 94550
> 925-423-1062
> Lab cell 925-724-7509
>
>
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Combine by columns a vector with another vector that is constant across rows

2018-07-03 Thread Gabor Grothendieck
or this variation if you don't want the first column to be named init:

 Reduce(cbind2, vec, 1:5)

On Tue, Jul 3, 2018 at 10:46 AM, Gabor Grothendieck
 wrote:
> Try Reduce:
>
>   Reduce(cbind, vec, 1:5)
>
> On Tue, Jul 3, 2018 at 9:28 AM, Viechtbauer, Wolfgang (SP)
>  wrote:
>> Hi All,
>>
>> I have one vector that I want to combine with another vector and that other 
>> vector should be the same for every row in the combined matrix. This 
>> obviously does not work:
>>
>> vec <- c(2,4,3)
>> cbind(1:5, vec)
>>
>> This does, but requires me to specify the correct value for 'n' in 
>> replicate():
>>
>> cbind(1:5, t(replicate(5, vec)))
>>
>> Other ways that do not require this are:
>>
>> t(sapply(1:5, function(x) c(x, vec)))
>> do.call(rbind, lapply(1:5, function(x) c(x, vec)))
>> t(mapply(c, 1:5, MoreArgs=list(vec)))
>>
>> I wonder if there is a simpler / more efficient way of doing this.
>>
>> Best,
>> Wolfgang
>>
>> __
>> 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.
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Combine by columns a vector with another vector that is constant across rows

2018-07-03 Thread Gabor Grothendieck
Try Reduce:

  Reduce(cbind, vec, 1:5)

On Tue, Jul 3, 2018 at 9:28 AM, Viechtbauer, Wolfgang (SP)
 wrote:
> Hi All,
>
> I have one vector that I want to combine with another vector and that other 
> vector should be the same for every row in the combined matrix. This 
> obviously does not work:
>
> vec <- c(2,4,3)
> cbind(1:5, vec)
>
> This does, but requires me to specify the correct value for 'n' in 
> replicate():
>
> cbind(1:5, t(replicate(5, vec)))
>
> Other ways that do not require this are:
>
> t(sapply(1:5, function(x) c(x, vec)))
> do.call(rbind, lapply(1:5, function(x) c(x, vec)))
> t(mapply(c, 1:5, MoreArgs=list(vec)))
>
> I wonder if there is a simpler / more efficient way of doing this.
>
> Best,
> Wolfgang
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Extract function parameters from a R expression

2018-06-20 Thread Gabor Grothendieck
If you specifically want to know which packages were loaded by the script
then using a vanilla version of R (i.e. one where only base packages are
loaded):

  vanilla_search <- search()
  source("myRprg.R")
  setdiff(search(), vanilla_search)



On Wed, Jun 20, 2018 at 4:08 AM, Sigbert Klinke
 wrote:
> Hi,
>
> I have read an R program with
>
> expr <- parse("myRprg.R")
>
> How can I extract the parameters of a specifc R command, e.g. "library"?
>
> So, if myprg.R containes the lines
>
> library("xyz")
> library("abc")
>
> then I would like to get "xyz" and "abc" back from expr.
>
> Thanks in advance
>
> Sigbert
>
> --
> https://hu.berlin/sk
>
>
>
> __
> 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.
>



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Take average of previous weeks

2018-03-25 Thread Gabor Grothendieck
There is no  `value` column in the `dput` output shown in the
question so using `tmin` instead note that the `width=` argument
of `rollapply` can be a list containing a vector of offsets (-1 is prior
value, -2 is value before that, etc.) and that we can use `rollapplyr`
with an `r` on the end to get right alignment.  See `?rollapply`

  library(dplyr)
  library(zoo)

  roll <- function(x, k) rollapplyr(x, list(-seq(1:k)), mean, fill = NA)
  df %>%
  group_by(citycode) %>%
  mutate(mean2 = roll(tmin, 2), mean3 = roll(tmin, 3), mean4 =
roll(tmin, 4)) %>%
  ungroup

(The code above has been indented 2 spaces so you can
identify inadvertent line wrapping by the email system.)


On Sun, Mar 25, 2018 at 10:48 AM, Miluji Sb  wrote:
> Dear all,
>
> I have weekly data by city (variable citycode). I would like to take the
> average of the previous two, three, four weeks (without the current week)
> of the variable called value.
>
> This is what I have tried to compute the average of the two previous weeks;
>
> df = df %>%
>   mutate(value.lag1 = lag(value, n = 1)) %>%
>   mutate(value .2.previous = rollapply(data = value.lag1,
>  width = 2,
>  FUN = mean,
>  align = "right",
>  fill = NA,
>  na.rm = T))
>
> I crated the lag of the variable first and then attempted to compute the
> average but this does not seem to to what I want. What I am doing wrong?
> Any help will be appreciated. The data is below. Thank you.
>
> Sincerely,
>
> Milu
>
> dput(droplevels(head(df, 10)))
> structure(list(year = c(1970L, 1970L, 1970L, 1970L, 1970L, 1970L,
> 1970L, 1970L, 1970L, 1970L), citycode = c(1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L), month = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
> 2L, 3L), week = c(1L, 2L, 3L, 4L, 5L, 5L, 6L, 7L, 8L, 9L), date =
> structure(c(1L,
> 2L, 3L, 4L, 5L, 5L, 6L, 7L, 8L, 9L), .Label = c("1970-01-10",
> "1970-01-17", "1970-01-24", "1970-01-31", "1970-02-07", "1970-02-14",
> "1970-02-21", "1970-02-28", "1970-03-07"), class = "factor"),
> value = c(-15.035, -20.478, -22.245, -23.576, -8.840995,
> -18.497, -13.892, -18.974, -15.919, -13.576)), .Names = c("year",
> "citycode", "month", "week", "date", "tmin"), row.names = c(NA,
> 10L), class = "data.frame")
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Equivalent of gtools::mixedsort in R base

2018-03-12 Thread Gabor Grothendieck
split any mixed columns into letter and number columns
and then order can be used on that:

  DF <- data.frame(x = c("a10", "a2", "a1"))
  o <- do.call("order", transform(DF, let = gsub("\\d", "", x),
 no =
as.numeric(gsub("\\D", "", x)),
 x = NULL))
  DF[o,, drop = FALSE ]


On Mon, Mar 12, 2018 at 12:15 AM, Sebastien Bihorel
 wrote:
> Hi,
>
> Searching for functions that would order strings that mix characters and 
> numbers in a "natural" way (ie, "a1 a2 a10" instead of "a1 a10 a2"), I found 
> the mixedsort and mixedorder from the gtools package.
>
> Problems:
> 1- mixedorder does not work in a "do.call(mixedorder, mydataframe)" call like 
> the order function does
> 2- gtools has not been updated in 2.5 years
>
> Are you aware of an equivalent of this function in base R or a another 
> contributed package (with correction of problem #1)?
>
> Thanks
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Plotting quarterly time series

2018-01-28 Thread Gabor Grothendieck
Using Achim's d this also works to generate z where FUN is a function used
to transform the index column and format is also passed to FUN.

z <- read.zoo(d, index = "time", FUN = as.yearqtr, format = "Q%q %Y")

On Sun, Jan 28, 2018 at 4:53 PM, Achim Zeileis  wrote:
> On Sun, 28 Jan 2018, p...@philipsmith.ca wrote:
>
>> I have a data set with quarterly time series for several variables. The
>> time index is recorded in column 1 of the dataframe as a character vector
>> "Q1 1961", "Q2 1961","Q3 1961", "Q4 1961", "Q1 1962", etc. I want to produce
>> line plots with ggplot2, but it seems I need to convert the time index from
>> character to date class. Is that right? If so, how do I make the conversion?
>
>
> You can use the yearqtr class in the zoo package, converting with
> as.yearqtr(..., format = "Q%q %Y"). zoo also provides an autoplot() method
> for ggplot2-based time series visualizations. See ?autoplot.zoo for various
> examples.
>
> ## example data similar to your description
> d <- data.frame(sin = sin(1:8), cos = cos(1:8))
> d$time <- c("Q1 1961", "Q2 1961", "Q3 1961", "Q4 1961", "Q1 1962",
>   "Q2 1962", "Q3 1962", "Q4 1962")
>
> ## convert to zoo series
> library("zoo")
> z <- zoo(as.matrix(d[, 1:2]), as.yearqtr(d$time, "Q%q %Y"))
>
> ## ggplot2 display
> library("ggplot2")
> autoplot(z)
>
> ## with nicer axis scaling
> autoplot(z) + scale_x_yearqtr()
>
> ## some variations
> autoplot(z, facets = Series ~ .) + scale_x_yearqtr() + geom_point()
> autoplot(z, facets = NULL) + scale_x_yearqtr() + geom_point()
>
>
>
>> __
>> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Portable R in zip file for Windows

2018-01-25 Thread Gabor Grothendieck
Can you clarify what the nature of the security restriction is?

If you can't run the R installer then how it is that you could run R?
That would still involve running an external exe even if it came
in a zip file.

Could it be that the restriction is not on running exe files but on
downloading them?

If that is it then there are obvious workarounds (rename it not
to have an exe externsion or zip it using another machine,
upload to the cloud and download onto the restricted machine)
but it might be safer to just ask the powers that be to download it
for you.  You probably don't need a new version of R more than
once a year.


On Thu, Jan 25, 2018 at 3:04 PM, Juan Manuel Truppia
<jmtrup...@gmail.com> wrote:
> What is wrong with you guys? I asked for a zip, like R Studio has for
> example. Totally clear.
> I cant execute exes. But I can unzip files.
> Thanks Gabor, I had that in mind, but can't execute the exe due to security
> restrictions.
> Geez, really, treating people who ask questions this way just makes you
> don't want to ask a single one.
>
>
> On Thu, Jan 25, 2018, 11:19 Gabor Grothendieck <ggrothendi...@gmail.com>
> wrote:
>>
>> I believe that the ordinary Windows installer for R can produce a
>> portable result by choosing the appropriate configuration options from the
>> offered screens when you run the installer  Be sure to enter the desired
>> path in the Select Destination Location screen, choose Yes on the
>> Startup options screen and ensure that all boxes are unchecked on the
>> Select additional tasks screen.
>>
>> On Wed, Jan 24, 2018 at 10:11 PM, Juan Manuel Truppia
>> <jmtrup...@gmail.com> wrote:
>> > I read a message from 2009 or 2010 where it mentioned the availability
>> > of R
>> > for Windows in a zip file, no installation required. It would be very
>> > useful for me. Is this still available somewhere?
>> >
>> > Thanks
>> >
>> > [[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.
>>
>>
>>
>> --
>> Statistics & Software Consulting
>> GKX Group, GKX Associates Inc.
>> tel: 1-877-GKX-GROUP
>> email: ggrothendieck at gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Portable R in zip file for Windows

2018-01-25 Thread Gabor Grothendieck
I believe that the ordinary Windows installer for R can produce a
portable result by choosing the appropriate configuration options from the
offered screens when you run the installer  Be sure to enter the desired
path in the Select Destination Location screen, choose Yes on the
Startup options screen and ensure that all boxes are unchecked on the
Select additional tasks screen.

On Wed, Jan 24, 2018 at 10:11 PM, Juan Manuel Truppia
 wrote:
> I read a message from 2009 or 2010 where it mentioned the availability of R
> for Windows in a zip file, no installation required. It would be very
> useful for me. Is this still available somewhere?
>
> Thanks
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] run r script in r-fiddle

2017-10-31 Thread Gabor Grothendieck
Try that source statement here -- it is running R 3.4.1:

https://www.tutorialspoint.com/execute_r_online.php

On Mon, Oct 30, 2017 at 11:14 AM, Suzen, Mehmet  wrote:
>  Note that, looks like r-fiddle runs R 3.1.2.
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Time series: xts/zoo object at annual (yearly) frequency

2017-10-06 Thread Gabor Grothendieck
Maybe one of these are close enough:

  xts(c(2, 4, 5), yearqtr(1991:1993))

  as.xts(ts(c(2, 4, 5), 1991))

of if you want only a plain year as the index then then use zoo,
zooreg or ts class:

  library(zoo)
  zoo(c(2, 4, 5), 1991:1993)

  zooreg(c(2, 4, 5), 1991)

  ts(c(2, 4, 5), 1991)

On Fri, Oct 6, 2017 at 2:56 AM, John  wrote:
> Hi,
>
>I'd like to make a time series at an annual frequency.
>
>> a<-xts(x=c(2,4,5), order.by=c("1991","1992","1993"))
> Error in xts(x = c(2, 4, 5), order.by = c("1991", "1992", "1993")) :
>   order.by requires an appropriate time-based object
>> a<-xts(x=c(2,4,5), order.by=1991:1993)
> Error in xts(x = c(2, 4, 5), order.by = 1991:1993) :
>   order.by requires an appropriate time-based object
>
>   How should I do it? I know that to do for quarterly or monthly time
> series, we use as.yearqtr or as.yearmon. What about annual?
>
>Thanks,
>
> John
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] require help

2017-09-21 Thread Gabor Grothendieck
Assuming the input data.frame, DF, is of the form  shown reproducibly
in the Note below, to convert the series to zoo or ts:

library(zoo)

# convert to zoo
z <- read.zoo(DF)

# convert to ts
as.ts(z) #


Note:

DF <- structure(list(year = c(1980, 1981, 1982, 1983, 1984), cnsm = c(174,
175, 175, 172, 173), incm = c(53.4, 53.7, 53.5, 53.2, 53.3),
with = c(60.3, 60.5, 60.2, 60.1, 60.7)), .Names = c("year",
"cnsm", "incm", "with"), row.names = c(NA, -5L), class = "data.frame")


On Sat, Sep 16, 2017 at 8:10 AM, yadav neog  wrote:
> oky.. thank you very much to all of you
>
>
> On Sat, Sep 16, 2017 at 2:06 PM, Eric Berger  wrote:
>
>> You can just use the same code that I provided before but now use your
>> dataset. Like this
>>
>> df <- read.csv(file="data2.csv",header=TRUE)
>> dates <- as.Date(paste(df$year,"-01-01",sep=""))
>> myXts <- xts(df,order.by=dates)
>> head(myXts)
>>
>> #The last command "head(myXts)" shows you the first few rows of the xts
>> object
>>year cnsmincmwlth
>> 1980-01-01 1980 173.6527 53.3635 60.3013
>> 1981-01-01 1981 175.4613 53.6929 60.4980
>> 1982-01-01 1982 174.5724 53.4890 60.2358
>> 1983-01-01 1983 171.5070 53.2223 60.1047
>> 1984-01-01 1984 173.3462 53.2851 60.6946
>> 1985-01-01 1985 171.7075 53.1596 60.7598
>>
>>
>> On Sat, Sep 16, 2017 at 9:55 AM, Berend Hasselman  wrote:
>>
>>>
>>> > On 15 Sep 2017, at 11:38, yadav neog  wrote:
>>> >
>>> > hello to all. I am working on macroeconomic data series of India, which
>>> in
>>> > a yearly basis. I am unable to convert my data frame into time series.
>>> > kindly help me.
>>> > also using zoo and xts packages. but they take only monthly
>>> observations.
>>> >
>>> > 'data.frame': 30 obs. of  4 variables:
>>> > $ year: int  1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 ...
>>> > $ cnsm: num  174 175 175 172 173 ...
>>> > $ incm: num  53.4 53.7 53.5 53.2 53.3 ...
>>> > $ wlth: num  60.3 60.5 60.2 60.1 60.7 ...
>>> > --
>>>
>>> Second try to do what you would like (I hope and think)
>>> Using Eric's sample data
>>>
>>> 
>>> zdf <- data.frame(year=2001:2010, cnsm=sample(170:180,10,replace=TRUE),
>>>  incm=rnorm(10,53,1), wlth=rnorm(10,60,1))
>>> zdf
>>>
>>> # R ts
>>> zts <- ts(zdf[,-1], start=zdf[1,"year"])
>>> zts
>>>
>>> # turn data into a zoo timeseries and an xts timeseries
>>>
>>> library(zoo)
>>> z.zoo <- as.zoo(zts)
>>> z.zoo
>>>
>>> library(xts)
>>> z.xts <- as.xts(zts)
>>> z.xts
>>> 
>>>
>>> Berend Hasselman
>>>
>>> > Yadawananda Neog
>>> > Research Scholar
>>> > Department of Economics
>>> > Banaras Hindu University
>>> > Mob. 9838545073
>>> >
>>> >   [[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/posti
>>> ng-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/posti
>>> ng-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>
>>
>
>
> --
> Yadawananda Neog
> Research Scholar
> Department of Economics
> Banaras Hindu University
> Mob. 9838545073
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] symbolic computing example with Ryacas

2017-09-19 Thread Gabor Grothendieck
Here are some more examples:

library(Ryacas)

x <- Sym("x")
yacas("x:=2")
Eval(x*x)
## [1] 4

# vignette has similar example
y <- Sym("y")
Eval(Subst(y*y, y, 3))
## [1] 9

# demo("Ryacas-Function") has similar example to this
f <- function(z) {}
body(f) <- yacas(expression(z*z))[[1]]
f(4)
## [1] 16



On Tue, Sep 19, 2017 at 2:08 PM, Vivek Sutradhara  wrote:
> Thanks for the response. Yes, I did study the vignette but did not
> understand it fully. Anyway, I have tried once again now. I am happy to say
> that I have got what I wanted.
>
> library(Ryacas)
> x <- Sym("x");U <- Sym("U");x0 <- Sym("x0");C <- Sym("C")
> my_func <- function(x,U,x0,C) {
>   return (U/(1+exp(-(x-x0)/C)))}
> FirstDeriv <- deriv(my_func(x,U,x0,C), x)
> PrettyForm(FirstDeriv)
> #slope <- yacas("Subst(x,x0),deriv(my_func(x,U,x0,C), x)")
> slope <- Subst(FirstDeriv,x,x0)
> #PrettyForm(slope) - gives errors
> PrettyForm(Simplify(slope))
>
> I was confused by the references to the yacas command.  Now, I have chosen
> to omit it. Then I get what I want.
> Thanks,
> Vivek
>
> 2017-09-19 16:04 GMT+02:00 Bert Gunter :
>
>> Have you studied the "Introduction to Ryacas" vignette that come with the
>> package?
>>
>> Cheers,
>> Bert
>>
>>
>>
>> Bert Gunter
>>
>> "The trouble with having an open mind is that people keep coming along and
>> sticking things into it."
>> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>>
>> On Tue, Sep 19, 2017 at 2:37 AM, Vivek Sutradhara 
>> wrote:
>>
>>> Hi all,
>>> I am trying to implement the following matlab code with Ryacas :
>>>
>>> syms U x x0 C
>>>
>>> d1=diff(U/(1+exp(-(x-x0)/C)),x);
>>>
>>> pretty(d1)
>>>
>>> d2=diff(U/(1+exp(-(x-x0)/C)),x,2);
>>>
>>> pretty(d2)
>>>
>>> solx2 = solve(d2 == 0, x, 'Real', true)
>>>
>>> pretty(solx2)
>>>
>>> slope2=subs(d1,solx2)
>>>
>>>
>>> I have tried the following :
>>>
>>> library(Ryacas)
>>>
>>> x <- Sym("x");U <- Sym("U");x0 <- Sym("x0");C <- Sym("C")
>>>
>>> my_func <- function(x,U,x0,C) {
>>>
>>>   return (U/(1+exp(-(x-x0)/C)))}
>>>
>>> FirstDeriv <- deriv(my_func(x,U,x0,C), x)
>>>
>>> PrettyForm(FirstDeriv)
>>>
>>> slope <- yacas("Subst(x,x0),deriv(my_func(x,U,x0,C), x)")
>>>
>>> PrettyForm(slope)
>>>
>>>
>>> I don't understand how I should use the Subst command. I want the slope of
>>> the first derivative at x=x0. How do I implement that?
>>>
>>> I would appreciate any help that I can get.
>>>
>>> Thanks,
>>>
>>> Vivek
>>>
>>> [[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/posti
>>> ng-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>
>>
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Convert data into zoo object using Performance analytics package

2017-09-18 Thread Gabor Grothendieck
Depending on how you created df maybe your code has the column names
wrong.  In any case these 4 alternatives all work.  Start a fresh R
session and then copy and paste this into it.

library(zoo)
u  <- "https://faculty.washington.edu/ezivot/econ424/sbuxPrices.csv;
fmt <- "%m/%d/%Y"

# 1
sbux1.z <- read.csv.zoo(u, FUN = as.yearmon, format = fmt)

# 2
df <- read.csv(u)
sbux2.z <- read.zoo(df, FUN = as.yearmon, format = fmt)

# 3
df <- read.csv(u)
names(head(df))
## [1] "Date"  "Adj.Close"
sbux3.z <- zoo(df$Adj.Close, as.yearmon(df$Date, fmt))

# 4
df <- read.csv(u)
sbux4.z <- zoo(df[[2]], as.yearmon(df[[1]], fmt))

On Mon, Sep 18, 2017 at 7:36 AM, Upananda Pani  wrote:
> Dear All,
>
> While i am trying convert data frame object to zoo object I am
> getting numeric(0) error in performance analytics package.
>
> The source code i am using from this website to learn r in finance:
> https://faculty.washington.edu/ezivot/econ424/returnCalculations.r
>
> # create zoo objects from data.frame objects
> dates.sbux = as.yearmon(sbux.df$Date, format="%m/%d/%Y")
> dates.msft = as.yearmon(msft.df$Date, format="%m/%d/%Y")
> sbux.z = zoo(x=sbux.df$Adj.Close, order.by=dates.sbux)
> msft.z = zoo(x=msft.df$Adj.Close, order.by=dates.msft)
> class(sbux.z)
> head(sbux.z)
>> head(sbux.z)
> Data:
> numeric(0)
>
> I will be grateful if anybody would like to guide me where i am making the
> mistake.
>
> With best regards,
> Upananda Pani
>
>
> --
>
>
> You may delay, but time will not.
>
>
> Research Scholar
> alternative mail id: up...@iitkgp.ac.in
> Department of HSS, IIT KGP
> KGP
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Case statement in sqldf

2017-09-11 Thread Gabor Grothendieck
2018-03-3 in your code should be 2018-03-31.

The line
then'201415'
needs to be fixed.

When posting please provide minimal self-contained examples. There was
no input provided and library statements not relevant to the posted
code were included.

Fixing the invalid date and bad line, getting rid of those library
statements that are unnecessary and providing some test input, it
works for me for the input shown.

(Note that it would NOT work if we omitted library(RH2) since the
default sqlite back end does not have date types and does not know
that an R date -- which is sent to sqlite as the number of days since
1970-01-01 -- corresponds to a particular character string; however,
the H2 database does have date types.  See FAQ #4 on the sqldf github
home page for more info.
https://github.com/ggrothendieck/sqldf
)

This works:

library(sqldf)
library(RH2)

cr <- data.frame(ReportDate = as.Date("2017-09-11")) # input

cr2 =  sqldf(" select ReportDate
 ,  case
   when ReportDate between  '2012-04-01'  and  '2013-03-31'
   then '2012_13'
   when  ReportDate between '2013-04-01'  and  '2014-03-31'
   then '2013_14'
   when  ReportDate between  '2014-04-01'  and  '2015-03-31'
   then '2014_15'
   when ReportDate between '2015-04-01'  and  '2016-03-31'
   then '2015_16'
   when ReportDate between '2016-04-01'  and  '2017-03-31'
   then '2016_17'
   when ReportDate between '2017-04-01'  and  '2018-03-31'
  then '2017_18' else null
 end as FY
 from cr
 where  ReportDate  >=  '2012-04-01'
 ")

giving:

  > cr2
ReportDate  FY
  1 2017-09-11 2017_18

Note that using as.yearqtr from zoo this alternative could be used:

library(zoo)
cr <- data.frame(ReportDate = as.Date("2017-09-11")) # input

fy <- as.integer(as.yearqtr(cr$ReportDate) + 3/4)
transform(cr, FY = paste0(fy-1, "_", fy %% 100))

giving:

  ReportDate  FY
1 2017-09-11 2017_18


On Mon, Sep 11, 2017 at 4:05 AM, Mangalani Peter Makananisa
 wrote:
> Hi all,
>
>
>
> I am trying to create a new  variable called Fiscal Year (FY) using case
> expression in sqldf  and I am getting a null FY , see the code below .
>
>
>> +then '2017_18' else null>> South African Revenue 
>> Service (SARS)>> Specialist: Statistical Support>> TCEI_OR (Head Office)>> 
>> Tell: +272 422 7357, Cell: +2782 456 4669>> 
>> http://www.sars.gov.za/Pages/Email-disclaimer.aspxemail: ggrothendieck at 
>> gmail.with
> Please advise me as to how I can do this mutation.
>
>
>
>   library(zoo)
>
>   library(lubridate)
>
>   library(stringr)
>
>   library(RH2)
>
>   library(sqldf)
>
>
>
> cr$ReportDate = as.Date(cr$ReportDate, format ='%Y-%m-%d')
>
>
>
>> cr2 =  sqldf(" select ReportDate
>
> +  ,  case
>
> +when ReportDate between  '2012-04-01'  and
> '2013-03-31'
>
> +then '2012_13'
>
> +when  ReportDate between '2013-04-01'  and
> '2014-03-31'
>
> +then '2013_14'
>
> +when  ReportDate between  '2014-04-01'  and
> '2015-03-31'
>
> +then'201415'
>
> +when ReportDate between '2015-04-01'  and
> '2016-03-31'
>
> +then '2015_16'
>
> +when ReportDate between '2016-04-01'  and
> '2017-03-31'
>
> +then '2016_17'
>
> +when ReportDate between '2017-04-01'  and
> '2018-03-3'
>


> +end as FY
>
> +   from cr
>
> +  where  ReportDate  >=  '2012-04-01'
>
> +  ")
>
>
>
> Thanking you in advance
>
>
>
> Kind regards,
>
>
>
> Mangalani Peter Makananisa (0005786)
>








>
>
>
>
> Disclaimer
>
> Please Note: This email and its contents are subject to our email legal
> notice which can be viewed at




-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP

__
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] Package sqldf in R and dates manipulation

2017-08-11 Thread Gabor Grothendieck
See FAQ #4 on the sqldf github home page.

On Fri, Aug 11, 2017 at 9:21 AM, Mangalani Peter Makananisa
 wrote:
> Dear all,
>
> I recently read the book " R data preperation and manipulation using sqldf 
> package"  by Djoni Darmawikarta
> However, I have a problem with manipulation of dates using this package, I do 
> not get the expected results. Do I need to install some packages  to be able 
> to subset the data by dates in sqldf?
>
> I am not getting Djoni Darmawikarta email address.
>
> Please see the practice data attached and advise,
>
> Kind regards,
>
> Mangalani Peter Makananisa (0005786)
> South African Revenue Service (SARS)
> Specialist: Statistical Support
> TCEI_OR (Head Office)
> Tell: +2712 422 7357, Cell: +2782 456 4669
>
>
> product = read.csv('D:/Users/S1033067/Desktop/sqldf prac/sqlprac.csv', 
> na.strings = '', header = T)
> head(product)
> library(sqldf)
> sqldf()
>
> # out-put
>
>> sqldf("select * from product
> + where (date(Launch_dt) >= date('01-01-2003'))
> +  ")
> [1] P_CodeP_NameLaunch_dt Price
> <0 rows> (or 0-length row.names)
>>
>
> Please Note: This email and its contents are subject to our email legal 
> notice which can be viewed at 
> http://www.sars.gov.za/Pages/Email-disclaimer.aspx
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] rollapply() produces NAs

2017-05-28 Thread Gabor Grothendieck
Maybe you want this.It computes VaRfun(r[c(i-500, i-1)] for each i for
which the argument to r makes sense.

rollapply(r, width = list(c(-500, -1)), FUN = VaRfun),

On Sat, May 27, 2017 at 5:29 PM, Sepp via R-help  wrote:
> Hello,
> I am fairly new to R and trying to calculate value at risk with exponentially 
> decreasing weights.My function works for a single vector of returns but does 
> not work with rollapply(), which is what I want to use. The function I am 
> working on should assig exponentially decreasing weights to the K most recent 
> returns and then order the returns in an ascending order. Subsequently it 
> should pick the last return for which the cumulative sum of the weights is 
> smaller or equal to a significance level. Thus, I am trying to construct a 
> cumulative distribution function and find a quantile.
> This is the function I wrote:
> VaRfun <- function(x, lambda = 0.94) {
> #create data.frame and order returns such that the lates return is the first  
> df <- data.frame(weight = c(1:length(x)), return = rev(x))  K <- nrow(df)  
> constant <- (1-lambda)/(1-lambda^(K))#assign weights to the returnsfor(i 
> in 1:nrow(df)) {df$weight[i] <- lambda^(i-1) * constant}#order 
> returns in an ascending order  df <- df[order(df$return),]
> #add the cumulative sum of the weights  df$cum.weight <- cumsum(df$weight)
> #calculate value at risk  VaR <- -tail((df$return[df$cum.weight <= .05]), 1)  
> signif(VaR, digits = 3)}
> It works for a single vector of returns but if I try to use it with 
> rollapply(), such as
> rollapply(r, width = list(-500, -1), FUN = VaRfun),
> it outputs a vector of NAs and I don't know why.
> Thank you for your help!
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

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


Re: [R] How to calculate Rolling mean for a List object

2017-05-14 Thread Gabor Grothendieck
Try this code which does not use rollapply:

w <- 3
Mean <- function(L) Reduce("+", L) / length(L)
lapply(w:length(Data), function(i) Mean(Data[seq(to = i, length = w)]))

On Sun, May 14, 2017 at 6:44 PM, Christofer Bogaso
 wrote:
> Hi again,
>
> I am looking to find a way on how to calculate Rolling average for the
> elements of a list. For example consider below object 'Data'. This is
> a list, where each elements are a Matrix. Basically, I am trying to
> get Rolling average of those Matrices with rolling window as 5.
>
> Data = structure(list(`2017-03-01` = structure(c(1.24915216491479e-06,
> -2.0209685810767e-06, -6.64165527006046e-06, -2.0209685810767e-06,
> 3.26966891657893e-06, 1.07453495291747e-05, -6.64165527006046e-06,
> 1.07453495291747e-05, 3.53132196103035e-05), .Dim = c(3L, 3L)),
> `2017-03-02` = structure(c(0.00863066441403338, -7.25585852047094e-05,
> -0.000950715788640005, -7.25585852047094e-05, 6.10004981580403e-07,
> 7.99273256915577e-06, -0.000950715788640005, 7.99273256915577e-06,
> 0.000104726642980084), .Dim = c(3L, 3L)), `2017-03-03` =
> structure(c(0.000785677680557358,
> 0.000283148300122928, 0.000170319078518317, 0.000283148300122928,
> 0.000102043066573597, 6.13808419844048e-05, 0.000170319078518317,
> 6.13808419844048e-05, 3.6921741860797e-05), .Dim = c(3L,
> 3L)), `2017-03-06` = structure(c(0.000100715163251975, 
> 1.80035062425799e-06,
> -5.05489732985851e-07, 1.80035062425799e-06, 3.21824665284709e-08,
> -9.03596565752718e-09, -5.05489732985851e-07, -9.03596565752718e-09,
> 2.53705461922188e-09), .Dim = c(3L, 3L)), `2017-03-07` =
> structure(c(0.000640065014281149,
> -0.000110994847091752, -0.000231235438845606, -0.000110994847091752,
> 1.92478198402357e-05, 4.00989612058198e-05, -0.000231235438845606,
> 4.00989612058198e-05, 8.35381203238728e-05), .Dim = c(3L,
> 3L)), `2017-03-08` = structure(c(7.72648041923266e-06,
> -2.11571338014623e-05,
> 7.82052544997182e-06, -2.11571338014623e-05, 5.79337921544145e-05,
> -2.14146538093767e-05, 7.82052544997182e-06, -2.14146538093767e-05,
> 7.91571517626794e-06), .Dim = c(3L, 3L)), `2017-03-09` =
> structure(c(4.43321118550061e-05,
> 1.90242249279913e-05, 5.68672547310199e-05, 1.90242249279913e-05,
> 8.16385953582618e-06, 2.44034267661023e-05, 5.68672547310199e-05,
> 2.44034267661023e-05, 7.29467766214148e-05), .Dim = c(3L,
> 3L)), `2017-03-10` = structure(c(0.000100081081692311, 
> 1.39245218598852e-05,
> 2.0935583168872e-05, 1.39245218598852e-05, 1.93735225227204e-06,
> 2.91281809264057e-06, 2.0935583168872e-05, 2.91281809264057e-06,
> 4.3794355057858e-06), .Dim = c(3L, 3L)), `2017-03-14` =
> structure(c(7.82185299651879e-06,
> -3.05963602958646e-05, -4.65590052688468e-05, -3.05963602958646e-05,
> 0.00011968228804236, 0.000182122586662866, -4.65590052688468e-05,
> 0.000182122586662866, 0.000277139058045361), .Dim = c(3L,
> 3L)), `2017-03-15` = structure(c(4.02156693772954e-05, 
> -2.2362610665311e-05,
> -2.08706726432905e-05, -2.2362610665311e-05, 1.24351120722764e-05,
> 1.16054944222453e-05, -2.08706726432905e-05, 1.16054944222453e-05,
> 1.08312253240602e-05), .Dim = c(3L, 3L)), `2017-03-16` =
> structure(c(2.64254966198469e-05,
> 5.78730550194069e-06, 5.0445603894268e-05, 5.78730550194069e-06,
> 1.26744656702641e-06, 1.10478196556107e-05, 5.0445603894268e-05,
> 1.10478196556107e-05, 9.62993804379875e-05), .Dim = c(3L,
> 3L)), `2017-03-17` = structure(c(0.000138433807049962, 
> 8.72005344938308e-05,
> 0.00014374477881467, 8.72005344938308e-05, 5.49282966209652e-05,
> 9.05459570205481e-05, 0.00014374477881467, 9.05459570205481e-05,
> 0.000149259504428865), .Dim = c(3L, 3L)), `2017-03-20` =
> structure(c(3.92058275846982e-05,
> 1.24332187386233e-05, -1.24235553811814e-05, 1.24332187386233e-05,
> 3.94290690251335e-06, -3.93984239286701e-06, -1.24235553811814e-05,
> -3.93984239286701e-06, 3.93678026502162e-06), .Dim = c(3L,
> 3L)), `2017-03-21` = structure(c(0.000407544227952838,
> -6.22427018306449e-05,
> 1.90596071859105e-05, -6.22427018306449e-05, 9.50609446890975e-06,
> -2.9109023406881e-06, 1.90596071859105e-05, -2.9109023406881e-06,
> 8.91360007491622e-07), .Dim = c(3L, 3L)), `2017-03-22` =
> structure(c(0.000220297355944482,
> 0.000282600064158173, 8.26030839524992e-05, 0.000282600064158173,
> 0.000362522718077154, 0.00010596421697645, 8.26030839524992e-05,
> 0.00010596421697645, 3.09729976068491e-05), .Dim = c(3L,
> 3L)), `2017-03-23` = structure(c(1.19559010537042e-05, 
> 3.56054556562106e-05,
> 5.51130473489473e-06, 3.56054556562106e-05, 0.000106035376739222,
> 1.64130261253175e-05, 5.51130473489473e-06, 1.64130261253175e-05,
> 2.54054292892148e-06), .Dim = c(3L, 3L)), `2017-03-24` =
> structure(c(0.000573948692221572,
> -7.36566239512158e-05, 5.40736580500709e-05, 

Re: [R] nlmrt problems - No confInt, NA StdErr, t-, or p-values

2017-03-21 Thread Gabor Grothendieck
You can use wrapnls from the nlmrt package to get an nls object.  Run
it instead of nlxb. It runs nlxb followed by nls so that the output is
an nls object..  Then you can use all of nls' methods.

On occiasion that fails even if nlxb succeeds since the nls
optimization can fail independently of nlxb.  Also, it does not show
the output from nlxb, only from the final nls, so you could
alternately run nlxb and then run nls2 from the nls2 package after
that.  nls2 can compute the nls object at a particular set of
coefficients so no second optimization that could fail is done. Here
is an example that uses nls2 to generate starting values for nlxb,
then runs nlxb and then uses nls2 again to get an nls object so that
it canthen  use nls methods (in this case fitted) on it.

http://stackoverflow.com/questions/42511278/nls-curve-fit-singular-matrix-error/42513058#42513058




On Tue, Mar 21, 2017 at 6:57 AM, DANIEL PRECIADO  wrote:
> Dear list,
>
> I want to use nlxb (package nlmrt) to fit different datasets to a gaussian, 
> obtain parameters (including standard error, t-and p-value) and confidence 
> intervals.
>
> nlxb generates the parameters, but very often results in NA standard 
> error,t-and p-values. Furthermore, using confint() to obtain the confidence 
> intervals generates a : Error in vcov.default(object) :   object does not 
> have variance-covariance matrix" erro.
>
> Can someone indicate why is nlxb generating NAs (when nls has no problem with 
> them) and how to obtain confidence intervals from an nlmrt object?
>
> Thanks
>
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Matrix

2017-03-07 Thread Gabor Grothendieck
Assuming that the input is x <- 1:4, try this one-liner:

> embed(c(0*x[-1], x, 0*x[-1]), 4)
 [,1] [,2] [,3] [,4]
[1,]1000
[2,]2100
[3,]3210
[4,]4321
[5,]0432
[6,]0043
[7,]0004

On Mon, Mar 6, 2017 at 11:18 AM, Peter Thuresson
 wrote:
> Hello,
>
> Is there a function in R which can transform, let say a vector:
>
> c(1:4)
>
> to a matrix where the vector is repeated but "projected" +1 one step down for 
> every (new) column.
> I want the output below from the vector above, like this:
>
> p<-c(1,2,3,4,0,0,0,0,1,2,3,4,0,0,0,0,1,2,3,4,0,0,0,0,1,2,3,4)
>
> matrix(p,7,4)
>

-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] abline with zoo series

2016-11-24 Thread Gabor Grothendieck
Recessions are typically shown by shading.  The zoo package has
xblocks for this purpose.  If app1 is your zoo object then:

plot(app1)
tt <- time(app1)
xblocks(tt, tt >= "1990-07-01" & tt <= "1991-03-31",
   col = rgb(0.7, 0.7, 0.7, 0.5)) # transparent grey

See ?xblocks for more info.



On Thu, Nov 24, 2016 at 10:03 PM, Erin Hodgess  wrote:
> Hello!  Happy Thanksgiving to those who are celebrating.
>
> I have a zoo series that I am plotting, and I would like to have some
> vertical lines at certain points, to indicate US business cycles.  Here is
> an example:
>
> app1 <- get.hist.quote(instrument="appl",
> start="1985-01-01",end="2016-08-31", quote="AdjClose", compression="m")
> #Fine
> plot(app1,main="Historical Stock Prices: Apple Corporation")
> #Still Fine
> #Now I want to use abline at July 1990 and March 1991 (as a start) for
> business cycles.  I tried v=67 and v="1990-07", no good.
>
> I have a feeling that it's really simple and I'm just not seeing it.
>
> Any help much appreciated.
>
> Thanks,
> Erin
>
>
> --
> Erin Hodgess
> Associate Professor
> Department of Mathematical and Statistics
> University of Houston - Downtown
> mailto: erinm.hodg...@gmail.com
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Split strings based on multiple patterns

2016-10-15 Thread Gabor Grothendieck
Replace newlines and colons with a space since they seem to be junk,
generate a pattern to replace the attributes with a comma and do the
replacement and finally read in what is left into a data frame using
the attributes as column names.

(I have indented each line of code below by 2 spaces so if any line
starts before that then it's been wrapped around by the email and
needs to be adjusted.)

  attributes <-
  c("Water temp", "Waterbody type", "Water pH", "Conductivity",
  "Water color", "Water turbidity", "Manmade", "Permanence", "Max water depth",
  "Primary substrate", "Evidence of cattle grazing", "Shoreline
Emergent Veg(%)",
  "Fish present", "Fish species")

  ugly2 <- gsub("[:\n]", " ", ugly)

  pat <- paste(gsub("([[:punct:]])", ".", attributes), collapse = "|")
  ugly3 <- gsub(pat, ",", ugly2)

  dd <- read.table(text = ugly3, sep = ",", strip.white = TRUE,
col.names = c("", attributes))[-1]


On Fri, Oct 14, 2016 at 7:16 PM, Joe Ceradini  wrote:
> Afternoon,
>
> I unfortunately inherited a dataframe with a column that has many fields
> smashed together. My goal is to split the strings in the column into
> separate columns based on patterns.
>
> Example of what I'm working with:
>
> ugly <- c("Water temp:14: F Waterbody type:Permanent Lake/Pond: Water
> pH:Unkwn:
> Conductivity:Unkwn: Water color: Clear: Water turbidity: clear:
> Manmade:no  Permanence:permanent:  Max water depth: <3: Primary
> substrate: Silt/Mud: Evidence of cattle grazing: none:
> Shoreline Emergent Veg(%): 1-25: Fish present: yes: Fish species: unkwn: no
> amphibians observed")
> ugly
>
> Far as I can tell, there is not a single pattern that would work for
> splitting this string. Splitting on ":" is close but not quite consistent.
> Each of these attributes should be a separate column:
>
> attributes <- c("Water temp", "Waterbody type", "Water pH", "Conductivity",
> "Water color", "Water turbidity", "Manmade", "Permanence", "Max water
> depth", "Primary substrate", "Evidence of cattle grazing", "Shoreline
> Emergent Veg(%)", "Fish present", "Fish species")
>
> So, conceptually, I want to do something like this, where the string is
> split for each of the patterns in attributes. However, strsplit only uses
> the 1st value of attributes
> strsplit(ugly, attributes)
>
> Should I loop through the values of "attributes"?
> Is there an argument in strsplit I'm missing that will do what I want?
> Different approach altogether?
>
> Thanks! Happy Friday.
> Joe
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Finding starting values for the parameters using nls() or nls2()

2016-10-09 Thread Gabor Grothendieck
If you are not tied to that model the SSasymp() model in R could be
considered and is easy to fit:

# to plot points in order
o <- order(cl$Area)
cl.o <- cl[o, ]

fm <- nls(Retention ~ SSasymp(Area, Asym, R0, lrc), cl.o)
summary(fm)

plot(Retention ~ Area, cl.o)
lines(fitted(fm) ~ Area, cl.o, col = "red")

__
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] Fwd: stfrtime function not returning proper results through sqldf package in R

2016-09-16 Thread Gabor Grothendieck
1. Convert the date from R's origin to the origin used by SQLite's
strftime function and then be sure you are using the correct SQLite
strftime syntax:

library(sqldf)
sqldf("select strftime('%m', Date + 2440588.5) month from log")

2. Alternately use the H2 backend which actually supports dates unlike
SQLite which only supports functions that interpret certain numbers
and character strings as dates.

library(RH2) # if RH2 is loaded sqldf will default to the H2 database
library(sqldf)
sqldf("select month(Date) from log")

Note that the first time you use sqldf with RH2 in a session it will
load java which is time consuming but after that it should run ok.

Note: See

1. the sqldf home page which has more info on dates and times:
https://github.com/ggrothendieck/sqldf

2. the sqldf help page:
?sqldf

3. the SQLite date and time function page which explains SQLite's
strftime function
https://www.sqlite.org/lang_datefunc.html

4. the H2 documentation:
http://www.h2database.com

On Fri, Sep 16, 2016 at 12:35 AM, Manohar Reddy  wrote:
> Hi ,
>
>
>
>   I have data something looks like below (or PFA), but when I’m extracting
> month  using *strftime*  function through *sqldf* library ,it’s returning
> below results but it’s not returning exact results ,it supposed to return
>  05,05,05,06,06,06.Can anyone please guide me how to do that with *strftime*
> function.
>
>
>
> Thanks in advance.
>
>
>
> Quiries :
>
>
> library(scales)
>
> # load data:
> log <- data.frame(Date =
> c("2013/05/25","2013/05/28","2013/05/31","2013/06/01","2013/06/02","2013/06/05","2013/06/07"),
>   Quantity = c(9,1,15,4,5,17,18))
>
>
> # convert date variable from factor to date format:
> log$Date <- as.Date(log$Date,
>   "%Y/%m/%d") # tabulate all the options here
> str(log)
>
>
>
>
>
>
>
> Manu.
>
> __
> 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.



--
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com


-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] stfrtime function not returning proper results through sqldf package in R

2016-09-16 Thread Gabor Grothendieck
To be precise it's SQLite that does not have date and time data types.
If you use an sqldf backend such as H2 that does have such types then
sqldf will pass them as such.  In the case of R's "Date" class such
objects are passed to SQLite as numbers since that is what SQLite can
understand but they are passed as dates to H2 (and other backends if
they have a date type).

Also note that after sqldf passes a date to SQLite as a number then
when SQLite passes it back to sqldf then sqldf  knows that it was
originally a date (due to the column name being the same) and coerces
it to Date class again.  In the example below sqldf passed the number
of days to 2000-01-01 since the UNIX epoch to SQLite.  SQLite then
processed it and passed it back as a number again. Then sqldf realized
that it was originally a Date because the column name is still d and
the original column d was of "Date" class and so coerces the number
from SQLite to Date.  There are a limited number of circumstances
where this heuristic works but they are sufficient that it's often
transparent even though SQLite has no date and time types.

> library(sqldf)
> DF <- data.frame(d = as.Date("2000-01-01"))
> sqldf("select d+1 as d from DF") # return next day
   d
1 2000-01-02

At the same time iif you really need to do serious date processing on
the SQL side it's much easier with an sqldf backend such as H2 that
actually supports date and time types and no heuristic is needed by
sqldf and no user workarounds on the R side are needed.


On Fri, Sep 16, 2016 at 10:15 AM, Jeff Newmiller
 wrote:
> SQLite only understands certain fundamental data types, and neither Date nor 
> POSIXct types are among them. They get stored as their internal numeric 
> representations.
>
> The internal numeric representations of Date and POSIXct are incompatible. 
> You are sending Dates to SQLite and trying to then interpret it as POSIXct by 
> handing that numeric to strftime.
>
> Note that within R the Date and POSIXct types are made sort-of compatible by 
> internal checking of class attributes that are not stored in SQLite. They are 
> still only sort-of compatible because Date has no concept of time zone and 
> always assumes GMT rather than local time when being converted.
>
> I recommend retrieving the stored Date value as a Date value into R so that 
> strftime can recognize how to interpret it. If you need to handle time as 
> well as date you may find that converting to character first before 
> converting to POSIXct with an appropriate time zone behaves with least 
> surprises.
> --
> Sent from my phone. Please excuse my brevity.
>
> On September 16, 2016 6:23:48 AM PDT, PIKAL Petr  
> wrote:
>>Hi Peter
>>
>>The devil is in detail
>>
>>Data from OP had different format and was transferred to Date object by
>>as.Date, which results in incorrect values (and NA if not transferred)
>>df <- data.frame(Date =
>>c("2013/05/25","2013/05/28","2013/05/31","2013/06/01","2013/06/02",
>>"2013/06/05","2013/06/07"), Quantity = c(9,1,15,4,5,17,18))
>>df$Date<-as.Date(df$Date)
>>cbind(df, sqldf("select strftime( '%m', Date) from df"))
>>
>>Data formatted according to your example transferred to Date object by
>>as data, again incorrect result
>>df2 <- data.frame(Date =
>>c("2013-05-25","2013-05-28","2013-05-31","2013-06-01","2013-06-02",
>>"2013-06-05","2013-06-07"), Quantity = c(9,1,15,4,5,17,18))
>>df2$Date<-as.Date(df2$Date)
>>cbind(df2, sqldf("select strftime( '%m', Date) from df2"))
>>
>>Data formatted according to your example but **not** changed to Dates,
>>correct result
>>df3 <- data.frame(Date =
>>c("2013-05-25","2013-05-28","2013-05-31","2013-06-01","2013-06-02",
>>"2013-06-05","2013-06-07"), Quantity = c(9,1,15,4,5,17,18))
>>cbind(df3, sqldf("select strftime( '%m', Date) from df3"))
>>
>>so sqldf is a bit peculiar about required input values and does not
>>know how to handle Date objects.
>>
>>Cheers
>>Petr
>>
>>
>>> -Original Message-
>>> From: peter dalgaard [mailto:pda...@gmail.com]
>>> Sent: Friday, September 16, 2016 2:45 PM
>>> To: PIKAL Petr 
>>> Cc: Manohar Reddy ; R-help >> project.org>
>>> Subject: Re: [R] stfrtime function not returning proper results
>>through sqldf
>>> package in R
>>>
>>> Presumably, sqldf does not know about Date object so passes an
>>integer that
>>> gets interpreted as who knows what...
>>>
>>> This seems to work:
>>>
>>> > df <- data.frame(date=as.character(Sys.Date()+seq(0,180,,10)))
>>> > cbind(df, sqldf("select strftime( '%m', date) from df"))
>>>  date strftime( '%m', date)
>>> 1  2016-09-1609
>>> 2  2016-10-0610
>>> 3  2016-10-2610
>>> 4  2016-11-1511
>>> 5  2016-12-0512
>>> 6  2016-12-2512
>>> 7  2017-01-1401
>>> 8  2017-02-0302
>>> 9  

Re: [R] C/C++/Fortran Rolling Window Regressions

2016-07-21 Thread Gabor Grothendieck
I would be careful about making assumptions regarding what is faster.
Performance tends to be nonintuitive.

When I ran rollapply/lm, rollapply/fastLm and roll_lm on the example
you provided rollapply/fastLm was three times faster than roll_lm.  Of
course this could change with data of different dimensions but it
would be worthwhile to do actual benchmarks before making assumptions.

I also noticed that roll_lm did not give the same coefficients as the other two.

set.seed(1)
library(zoo)
library(RcppArmadillo)
library(roll)
z <- zoo(matrix(rnorm(10), ncol = 2))
colnames(z) <- c("y", "x")

## rolling regression of width 4
library(rbenchmark)
benchmark(fastLm = rollapplyr(z, width = 4,
 function(x) coef(fastLm(cbind(1, x[, 2]), x[, 1])),
 by.column = FALSE),
   lm = rollapplyr(z, width = 4,
 function(x) coef(lm(y ~ x, data = as.data.frame(x))),
 by.column = FALSE),
   roll_lm =  roll_lm(coredata(z[, 1, drop = F]), coredata(z[, 2, drop = F]), 4,
 center = FALSE))[1:4]


 test replications elapsed relative
1  fastLm  1000.221.000
2  lm  1000.723.273
3 roll_lm  1000.642.909

On Thu, Jul 21, 2016 at 3:45 PM, jeremiah rounds
 wrote:
>  Thanks all.  roll::roll_lm was essentially what I wanted.   I think maybe
> I would prefer it to have options to return a few more things, but it is
> the coefficients, and the remaining statistics you might want can be
> calculated fast enough from there.
>
>
> On Thu, Jul 21, 2016 at 12:36 PM, Achim Zeileis 
> wrote:
>
>> Jeremiah,
>>
>> for this purpose there are the "roll" and "RcppRoll" packages. Both use
>> Rcpp and the former also provides rolling lm models. The latter has a
>> generic interface that let's you define your own function.
>>
>> One thing to pay attention to, though, is the numerical reliability.
>> Especially on large time series with relatively short windows there is a
>> good chance of encountering numerically challenging situations. The QR
>> decomposition used by lm is fairly robust while other more straightforward
>> matrix multiplications may not be. This should be kept in mind when writing
>> your own Rcpp code for plugging it into RcppRoll.
>>
>> But I haven't check what the roll package does and how reliable that is...
>>
>> hth,
>> Z
>>
>>
>> On Thu, 21 Jul 2016, jeremiah rounds wrote:
>>
>> Hi,
>>>
>>> A not unusual task is performing a multiple regression in a rolling window
>>> on a time-series.A standard piece of advice for doing in R is
>>> something
>>> like the code that follows at the end of the email.  I am currently using
>>> an "embed" variant of that code and that piece of advice is out there too.
>>>
>>> But, it occurs to me that for such an easily specified matrix operation
>>> standard R code is really slow.   rollapply constantly returns to R
>>> interpreter at each window step for a new lm.   All lm is at its heart is
>>> (X^t X)^(-1) * Xy,  and if you think about doing that with Rcpp in rolling
>>> window you are just incrementing a counter and peeling off rows (or
>>> columns
>>> of X and y) of a particular window size, and following that up with some
>>> matrix multiplication in a loop.   The psuedo-code for that Rcpp
>>> practically writes itself and you might want a wrapper of something like:
>>> rolling_lm (y=y, x=x, width=4).
>>>
>>> My question is this: has any of the thousands of R packages out there
>>> published anything like that.  Rolling window multiple regressions that
>>> stay in C/C++ until the rolling window completes?  No sense and writing it
>>> if it exist.
>>>
>>>
>>> Thanks,
>>> Jeremiah
>>>
>>> Standard (slow) advice for "rolling window regression" follows:
>>>
>>>
>>> set.seed(1)
>>> z <- zoo(matrix(rnorm(10), ncol = 2))
>>> colnames(z) <- c("y", "x")
>>>
>>> ## rolling regression of width 4
>>> rollapply(z, width = 4,
>>>   function(x) coef(lm(y ~ x, data = as.data.frame(x))),
>>>   by.column = FALSE, align = "right")
>>>
>>> ## result is identical to
>>> coef(lm(y ~ x, data = z[1:4,]))
>>> coef(lm(y ~ x, data = z[2:5,]))
>>>
>>> [[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.
>>>
>>>
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at 

Re: [R] C/C++/Fortran Rolling Window Regressions

2016-07-21 Thread Gabor Grothendieck
Just replacing lm with a faster version would speed it up.  Try lm.fit
or even faster is fastLm in the RcppArmadillo package.

On Thu, Jul 21, 2016 at 2:02 PM, jeremiah rounds
 wrote:
> Hi,
>
> A not unusual task is performing a multiple regression in a rolling window
> on a time-series.A standard piece of advice for doing in R is something
> like the code that follows at the end of the email.  I am currently using
> an "embed" variant of that code and that piece of advice is out there too.
>
> But, it occurs to me that for such an easily specified matrix operation
> standard R code is really slow.   rollapply constantly returns to R
> interpreter at each window step for a new lm.   All lm is at its heart is
> (X^t X)^(-1) * Xy,  and if you think about doing that with Rcpp in rolling
> window you are just incrementing a counter and peeling off rows (or columns
> of X and y) of a particular window size, and following that up with some
> matrix multiplication in a loop.   The psuedo-code for that Rcpp
> practically writes itself and you might want a wrapper of something like:
> rolling_lm (y=y, x=x, width=4).
>
> My question is this: has any of the thousands of R packages out there
> published anything like that.  Rolling window multiple regressions that
> stay in C/C++ until the rolling window completes?  No sense and writing it
> if it exist.
>
>
> Thanks,
> Jeremiah
>
> Standard (slow) advice for "rolling window regression" follows:
>
>
> set.seed(1)
> z <- zoo(matrix(rnorm(10), ncol = 2))
> colnames(z) <- c("y", "x")
>
> ## rolling regression of width 4
> rollapply(z, width = 4,
>function(x) coef(lm(y ~ x, data = as.data.frame(x))),
>by.column = FALSE, align = "right")
>
> ## result is identical to
> coef(lm(y ~ x, data = z[1:4,]))
> coef(lm(y ~ x, data = z[2:5,]))
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Concatenate two lists replacing elements with the same name.

2016-07-19 Thread Gabor Grothendieck
Try this:

Reduce(modifyList, list(x, y, z))

On Tue, Jul 19, 2016 at 12:34 PM, Luca Cerone  wrote:
> Dear all,
> I would like to know if there is a function to concatenate two lists
> while replacing elements with the same name.
>
> For example:
>
> x <- list(a=1,b=2,c=3)
> y <- list( b=4, d=5)
> z <- list(a = 6, b = 8, e= 7)
>
> I am looking for a function "concatfun" so that
>
> u <- concatfun(x,y,z)
>
> returns:
>
> u$a=6
> u$b=8
> u$c=3
> u$d=5
> u$e=7
>
> I.e. it combines the 3 lists, but when names have the same value it
> keeps the most recent one.
>
> Does such a function exists?
>
> Thanks for the help,
>
> Cheers,
> Luca
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] break string at specified possitions

2016-05-17 Thread Gabor Grothendieck
Here are two ways that do not use any packages:

s <- paste(letters, collapse = "") # test input

substring(s, first, last)
## [1] "abcde" "fghij" "klmnopqrs"


read.fwf(textConnection(s), last - first + 1)
## V1V2V3
## 1 abcde fghij klmnopqrs

On Wed, May 11, 2016 at 4:12 PM, Jan Kacaba  wrote:
> Dear R-help
>
> I would like to split long string at specified precomputed positions.
> 'substring' needs beginings and ends. Is there a native function which
> accepts positions so I don't have to count second argument?
>
> For example I have vector of possitions pos<-c(5,10,19). Substring
> needs input first=c(1,6,11) and last=c(5,10,19). There is no problem
> to write my own function. Just asking.
>
> Derek
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Clean method to convert date and time between time zones keeping it in POSIXct format

2016-05-09 Thread Gabor Grothendieck
This involves mucking with the internals as well but it is short:

   structure(T1, tzone = "UTC")

On Mon, May 9, 2016 at 9:24 AM, Arnaud Mosnier  wrote:
> Dear UseRs,
>
> I know two ways to convert dates and time from on time zone to another but
> I am pretty sure that there is a better (cleaner) way to do that.
>
>
> Here are the methods I know:
>
>
> ## The longest way ...
>
> T1 <- as.POSIXct("2016-05-09 10:00:00", format="%Y-%m-%d %H:%M:%S",
> tz="America/New_York")
>
> print(T1)
>
> T2 <- as.POSIXct(format(T1, tz="UTC"), tz="UTC") # format convert it to
> character, so I have to convert it back to POSIXct afterward.
>
> print(T2)
>
>
>
> ## The shortest but probably not the cleanest ...
>
> attributes(T1)$tzone <- "UTC"
>
> print(T1)
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Linear Regressions with constraint coefficients

2016-04-28 Thread Gabor Grothendieck
The nls2 package can be used to get starting values.

On Thu, Apr 28, 2016 at 8:42 AM, Aleksandrovic, Aljosa (Pfaeffikon)
<aljosa.aleksandro...@man.com> wrote:
> Hi Gabor,
>
> Thanks a lot for your help!
>
> I tried to implement your nonlinear least squares solver on my data set. I 
> was just wondering about the argument start. If I would like to force all my 
> coefficients to be inside an interval, let’s say, between 0 and 1, what kind 
> of starting values are normally recommended for the start argument (e.g. 
> Using a 4 factor model with b1, b2, b3 and b4, I tried start = list(b1 = 0.5, 
> b2 = 0.5, b3 = 0.5, b4 = 0.5))? I also tried other starting values ... Hence, 
> the outputs are very sensitive to that start argument?
>
> Thanks a lot for your answer in advance!
>
> Kind regards,
> Aljosa
>
>
>
> Aljosa Aleksandrovic, FRM, CAIA
> Quantitative Analyst - Convertibles
> aljosa.aleksandro...@man.com
> Tel +41 55 417 76 03
>
> Man Investments (CH) AG
> Huobstrasse 3 | 8808 Pfäffikon SZ | Switzerland
>
> -Original Message-
> From: Gabor Grothendieck [mailto:ggrothendi...@gmail.com]
> Sent: Dienstag, 26. April 2016 17:59
> To: Aleksandrovic, Aljosa (Pfaeffikon)
> Cc: r-help@r-project.org
> Subject: Re: [R] Linear Regressions with constraint coefficients
>
> This is a quadratic programming problem that you can solve using either a 
> quadratic programming solver with constraints or a general nonlinear solver 
> with constraints.  See https://cran.r-project.org/web/views/Optimization.html
> for more info on what is available.
>
> Here is an example using a nonlinear least squares solver and non-negative 
> bound constraints. The constraint that the coefficients sum to 1 is implied 
> by dividing them by their sum and then dividing the coefficients found by 
> their sum at the end:
>
> # test data
> set.seed(123)
> n <- 1000
> X1 <- rnorm(n)
> X2 <- rnorm(n)
> X3 <- rnorm(n)
> Y <- .2 * X1 + .3 * X2 + .5 * X3 + rnorm(n)
>
> # fit
> library(nlmrt)
> fm <- nlxb(Y ~ (b1 * X1 + b2 * X2 + b3 * X3)/(b1 + b2 + b3),
>  data = list(Y = Y, X1 = X1, X2 = X2, X3 = X3),
>  lower = numeric(3),
>  start = list(b1 = 1, b2 = 2, b3 = 3))
>
> giving the following non-negative coefficients which sum to 1 that are 
> reasonably close to the true values of 0.2, 0.3 and 0.5:
>
>> fm$coefficients / sum(fm$coefficients)
>  b1  b2  b3
> 0.18463 0.27887 0.53650
>
>
> On Tue, Apr 26, 2016 at 8:39 AM, Aleksandrovic, Aljosa (Pfaeffikon) 
> <aljosa.aleksandro...@man.com> wrote:
>> Hi all,
>>
>> I hope you are doing well?
>>
>> I’m currently using the lm() function from the package stats to fit linear 
>> multifactor regressions.
>>
>> Unfortunately, I didn’t yet find a way to fit linear multifactor regressions 
>> with constraint coefficients? I would like the slope coefficients to be all 
>> inside an interval, let’s say, between 0 and 1. Further, if possible, the 
>> slope coefficients should add up to 1.
>>
>> Is there an elegant and not too complicated way to do such a constraint 
>> regression estimation in R?
>>
>> I would very much appreciate if you could help me with my issue?
>>
>> Thanks a lot in advance and kind regards, Aljosa Aleksandrovic
>>
>>
>>
>> Aljosa Aleksandrovic, FRM, CAIA
>> Quantitative Analyst - Convertibles
>> aljosa.aleksandro...@man.com
>> Tel +41 55 417 7603
>>
>> Man Investments (CH) AG
>> Huobstrasse 3 | 8808 Pfäffikon SZ | Switzerland
>>
>>
>> -Original Message-
>> From: Kevin E. Thorpe [mailto:kevin.tho...@utoronto.ca]
>> Sent: Dienstag, 26. April 2016 14:35
>> To: Aleksandrovic, Aljosa (Pfaeffikon)
>> Subject: Re: Linear Regressions with constraint coefficients
>>
>> You need to send it to r-help@r-project.org however.
>>
>> Kevin
>>
>> On 04/26/2016 08:32 AM, Aleksandrovic, Aljosa (Pfaeffikon) wrote:
>>> Ok, will do! Thx a lot!
>>>
>>> Please find below my request:
>>>
>>> Hi all,
>>>
>>> I hope you are doing well?
>>>
>>> I’m currently using the lm() function from the package stats to fit linear 
>>> multifactor regressions.
>>>
>>> Unfortunately, I didn’t yet find a way to fit linear multifactor 
>>> regressions with constraint coefficients? I would like the slope 
>>> coefficients to be all inside an interval, let’s say, between 0 and 1. 
>>> Further, if possible, the slope coefficients should add up to 1.
>>>
>>> Is there

Re: [R] Approximate taylor series

2016-04-27 Thread Gabor Grothendieck
Regress on a multivariate polynomial:

lm(y ~ polym(x1, x2, x3, x4, degree = 3))

See ?polym

__
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] Linear Regressions with constraint coefficients

2016-04-26 Thread Gabor Grothendieck
This is a quadratic programming problem that you can solve using
either a quadratic programming solver with constraints or a general
nonlinear solver with constraints.  See
https://cran.r-project.org/web/views/Optimization.html
for more info on what is available.

Here is an example using a nonlinear least squares solver and
non-negative bound constraints. The constraint that the coefficients
sum to 1 is implied by dividing them by their sum and then dividing
the coefficients found by their sum at the end:

# test data
set.seed(123)
n <- 1000
X1 <- rnorm(n)
X2 <- rnorm(n)
X3 <- rnorm(n)
Y <- .2 * X1 + .3 * X2 + .5 * X3 + rnorm(n)

# fit
library(nlmrt)
fm <- nlxb(Y ~ (b1 * X1 + b2 * X2 + b3 * X3)/(b1 + b2 + b3),
 data = list(Y = Y, X1 = X1, X2 = X2, X3 = X3),
 lower = numeric(3),
 start = list(b1 = 1, b2 = 2, b3 = 3))

giving the following non-negative coefficients which sum to 1 that are
reasonably close to the true values of 0.2, 0.3 and 0.5:

> fm$coefficients / sum(fm$coefficients)
 b1  b2  b3
0.18463 0.27887 0.53650


On Tue, Apr 26, 2016 at 8:39 AM, Aleksandrovic, Aljosa (Pfaeffikon)
 wrote:
> Hi all,
>
> I hope you are doing well?
>
> I’m currently using the lm() function from the package stats to fit linear 
> multifactor regressions.
>
> Unfortunately, I didn’t yet find a way to fit linear multifactor regressions 
> with constraint coefficients? I would like the slope coefficients to be all 
> inside an interval, let’s say, between 0 and 1. Further, if possible, the 
> slope coefficients should add up to 1.
>
> Is there an elegant and not too complicated way to do such a constraint 
> regression estimation in R?
>
> I would very much appreciate if you could help me with my issue?
>
> Thanks a lot in advance and kind regards,
> Aljosa Aleksandrovic
>
>
>
> Aljosa Aleksandrovic, FRM, CAIA
> Quantitative Analyst - Convertibles
> aljosa.aleksandro...@man.com
> Tel +41 55 417 7603
>
> Man Investments (CH) AG
> Huobstrasse 3 | 8808 Pfäffikon SZ | Switzerland
>
>
> -Original Message-
> From: Kevin E. Thorpe [mailto:kevin.tho...@utoronto.ca]
> Sent: Dienstag, 26. April 2016 14:35
> To: Aleksandrovic, Aljosa (Pfaeffikon)
> Subject: Re: Linear Regressions with constraint coefficients
>
> You need to send it to r-help@r-project.org however.
>
> Kevin
>
> On 04/26/2016 08:32 AM, Aleksandrovic, Aljosa (Pfaeffikon) wrote:
>> Ok, will do! Thx a lot!
>>
>> Please find below my request:
>>
>> Hi all,
>>
>> I hope you are doing well?
>>
>> I’m currently using the lm() function from the package stats to fit linear 
>> multifactor regressions.
>>
>> Unfortunately, I didn’t yet find a way to fit linear multifactor regressions 
>> with constraint coefficients? I would like the slope coefficients to be all 
>> inside an interval, let’s say, between 0 and 1. Further, if possible, the 
>> slope coefficients should add up to 1.
>>
>> Is there an elegant and not too complicated way to do such a constraint 
>> regression estimation in R?
>>
>> I would very much appreciate if you could help me with my issue?
>>
>> Thanks a lot in advance and kind regards, Aljosa Aleksandrovic
>>
>>
>>
>> Aljosa Aleksandrovic, FRM, CAIA
>> Quantitative Analyst - Convertibles
>> aljosa.aleksandro...@man.com
>> Tel +41 55 417 7603
>>
>> Man Investments (CH) AG
>> Huobstrasse 3 | 8808 Pfäffikon SZ | Switzerland
>>
>>
>> -Original Message-
>> From: Kevin E. Thorpe [mailto:kevin.tho...@utoronto.ca]
>> Sent: Dienstag, 26. April 2016 14:28
>> To: Aleksandrovic, Aljosa (Pfaeffikon); r-help-ow...@r-project.org
>> Subject: Re: Linear Regressions with constraint coefficients
>>
>> I believe I approved a message with such a subject. Perhaps there was 
>> another layer that subsequently rejected it after that. I didn't notice any 
>> unusual content. Try again, making sure you send the message in plain text 
>> only.
>>
>> Kevin
>>
>> On 04/26/2016 08:16 AM, Aleksandrovic, Aljosa (Pfaeffikon) wrote:
>>> Do you know where I get help for my issue?
>>>
>>> Thanks in advance and kind regards,
>>> Aljosa
>>>
>>>
>>> Aljosa Aleksandrovic, FRM, CAIA
>>> Quantitative Analyst - Convertibles
>>> aljosa.aleksandro...@man.com
>>> Tel +41 55 417 7603
>>>
>>> Man Investments (CH) AG
>>> Huobstrasse 3 | 8808 Pfäffikon SZ | Switzerland
>>>
>>> -Original Message-
>>> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of
>>> r-help-ow...@r-project.org
>>> Sent: Dienstag, 26. April 2016 14:10
>>> To: Aleksandrovic, Aljosa (Pfaeffikon)
>>> Subject: Linear Regressions with constraint coefficients
>>>
>>> The message's content type was not explicitly allowed
>>>
>
>
> --
> Kevin E. Thorpe
> Head of Biostatistics,  Applied Health Research Centre (AHRC)
> Li Ka Shing Knowledge Institute of St. Michael's Hospital
> Assistant Professor, Dalla Lana School of Public Health
> University of Toronto
> email: kevin.tho...@utoronto.ca  Tel: 416.864.5776  Fax: 416.864.3016
>
> This email has been sent 

Re: [R] Functional programming?

2016-03-02 Thread Gabor Grothendieck
This manufactures the functions without using eval by using substitute
to substitute i-1 and a[i] into an expression for the body which is
then assigned to the body of the function:

hh <- vector("list", 5)
hh[[1]] <- f(a[1])
for(i in 2:5) {
 hh[[i]] <- hh[[1]]
 body(hh[[i]]) <- substitute(hh[[iprev]](x) * g(ai)(x), list(iprev
= i-1, ai = a[i]))
}
all.equal(h[[5]](.5), hh[[5]](.5)) # test
## [1] TRUE



This uses quote to define the body of h[[i]] as a call object and then
substitutes in the values of i-1 and a[i] assigning the result to the
body of h[[i]].

h <- vector("list", 5)
h[[1]] <- f(a[1])
for(i in 2:5) {
 h[[i]] <- h[[1]]
 body(hh[[i]]) <- do.call(substitute,
   list(quote(hh[[iprev]](x) *
g(ai)(x)),
   list(iprev = i-1, ai = a[i])))
}



On Wed, Mar 2, 2016 at 11:47 AM, Roger Koenker  wrote:
> I have a (remarkably ugly!!) code snippet  (below) that, given
> two simple functions, f and g,  generates
> a list of new functions h_{k+1} =  h_k * g, k= 1, …, K.  Surely, there are 
> vastly
> better ways to do this.  I don’t particularly care about the returned list,
> I’d be happy to have the final  h_K version of the function,
> but I keep losing my way and running into the dreaded:
>
> Error in h[[1]] : object of type 'closure' is not subsettable
> or
> Error: evaluation nested too deeply: infinite recursion / 
> options(expressions=)?
>
> Mainly I’d like to get rid of the horrible, horrible paste/parse/eval evils.  
> Admittedly
> the f,g look a bit strange, so you may have to suspend disbelief to imagine 
> that there is
> something more sensible lurking beneath this minimal (toy)  example.
>
> f <- function(u) function(x) u * x^2
> g <- function(u) function(x) u * log(x)
> set.seed(3)
> a <- runif(5)
> h <- list()
> hit <- list()
> h[[1]] <- f(a[1])
> hit[[1]] <- f(a[1])
> for(i in 2:5){
> ht <- paste("function(x) h[[", i-1, "]](x) * g(", a[i], ")(x)")
> h[[i]] <- eval(parse(text = ht))
> hit[[i]] <- function(x) {force(i); return(h[[i]] (x))}
> }
> x <- 1:99/10
> plot(x, h[[1]](x), type = "l")
> for(i in 2:5)
> lines(x, h[[i]](x), col = i)
>
> Thanks,
> Roger
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] sqldf --Warning message:

2016-02-19 Thread Gabor Grothendieck
sqldf does not use Tk so you can ignore this.

On Fri, Feb 19, 2016 at 12:32 PM, Divakar Reddy
 wrote:
> Dear R users,
>
> I'm getting Waring message while trying to load "sqldf" package in R3.2.3
> and assuming that we can ignore this as it's WARNING Message and not an
> error message.
> Can you guide me if my assumption is wrong?
>
>
>> library(sqldf);
> Loading required package: gsubfn
> Loading required package: proto
> Loading required package: RSQLite
> Loading required package: DBI
> Warning message:
> no DISPLAY variable so Tk is not available
>
>> version   _
> platform   x86_64-redhat-linux-gnu
> version.string R version 3.2.3 (2015-12-10)
>>
>
> Thanks,
> Divakar
> Phoenix,USA
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Create macro_var in R

2016-02-03 Thread Gabor Grothendieck
See

  Example 5.  Insert Variables

on the sqldf home page.

  https://github.com/ggrothendieck/sqldf


On Wed, Feb 3, 2016 at 2:16 PM, Amoy Yang via R-help
 wrote:
> First, MVAR<-c("population) should be the same as "population'". Correct?
> You use tab[[MVAR]] to refer to "population" where double [[...]] removes 
> double quotes "...", which seemingly work for r-code although it is tedious 
> in comparison direct application in SAS %let MVAR=population. But it does not 
> work for sqldef in R as shown below.
>
>> key<-"pop"
>> library(sqldf)
>> sqldf("select grade, count(*) as cnt, min(tab[[key]]) as min,
> + max(pop) as max, avg(pop) as mean, median(pop) as median,
> + stdev(pop) as stdev from tab group by grade")
> Error in sqliteSendQuery(con, statement, bind.data) :
>   error in statement: near "[[key]": syntax error
>
>
>
>
> On Wednesday, February 3, 2016 12:40 PM, "ruipbarra...@sapo.pt" 
>  wrote:
>
>
>  Hello,
>
> You can't use tab$MVAR but you can use tab[[MVAR]] if you do MVAR <- 
> "population" (no need for c()).
>
> Hope this helps,
>
> Rui Barradas
>  Citando Amoy Yang via R-help :
> population is the field-name in data-file (say, tab). MVAR<-population takes 
> data (in the column of population) rather than field-name as done in SAS:  
> %let MVAR=population;
> In the following r-program, for instance, I cannot use ... tab$MVAR...or 
> simply MVAR itself since MVAR is defined as "population" with double quotes 
> if using MVAR<-c("population")
>
>On Wednesday, February 3, 2016 11:54 AM, Duncan Murdoch 
>  wrote:
>
>
> On 03/02/2016 12:41 PM, Amoy Yang via R-help wrote:
>   There is a %LET statement in SAS: %let MVAR=population; Thus, MVAR can be 
> used through entire program.
> In R, I tried MAVR<-c("population"). The problem is that MAVR comes with 
> double quote "" that I don't need. But MVAR<-c(population) did NOT work 
> out. Any way that double quote can be removed as done in SAS when creating 
> macro_var?
> Thanks in advance for helps!
> R doesn't have a macro language, and you usually don't need one.
>
> If you are only reading the value of population, then
>
> MAVR <- population
>
> is fine.  This is sometimes the same as c(population), but in general
> it's different:  c() will remove some attributes, such as
> the dimensions on arrays.
>
> If you need to modify it in your program, it's likely more complicated.
> The normal way to go would be to put your code in a function, and have
> it return the modified version.  For example,
>
> population <- doModifications(population)
>
> where doModifications is a function with a definition like
>
> doModifications <- function(MAVR) {
> # do all your calculations on MAVR
> # then return it at the end using
> MAVR
> }
>
> Duncan Murdoch
>
>
>
> [[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.htmland provide commented, minimal, 
> self-contained, reproducible code.
>
>
>
>
> [[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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] fancy linear model and grouping

2016-02-02 Thread Gabor Grothendieck
Try the mclust package:

library(mclust)
temp.na <- na.omit(temp)
fm <- Mclust(temp.na)
g <- fm$classification
plot(temp.na, pch = g, col = g)



On Tue, Feb 2, 2016 at 6:35 AM, PIKAL Petr  wrote:
> Dear all
>
> I have data like this
>
>> dput(temp)
>
> temp <- structure(list(X1 = c(93, 82, NA, 93, 93, 79, 79, 93, 93, 85,
> 82, 93, 87, 93, 92, NA, 87, 93, 93, 93, 74, 77, 87, 93, 82, 87,
> 75, 82, 93, 92, 68, 93, 93, 73, NA, 85, 81, 79, 75, 87, 93, NA,
> 87, 87, 85, 92, 87, 92, 93, 87, 87, NA, 69, 87, 93, 87, 93, 87,
> 82, 79, 87, 93, 87, 80, 87, 87, 87, 92, 93, 69, 76, 87, 82, 93,
> 82, NA, 54, 87, 77, 73, 93, 82, 73, 93, 92, 82, 77, 93, 87, 75,
> 87, 87, 87, 60, 92, 87, 87, NA, 77, 78), X2 = c(224, 624, NA,
> 224, 224, 642, 642, 224, 224, 599, 622, 224, 239, 224, 225, NA,
> 239, 224, 224, 224, 688, 657, 239, 224, 624, 239, 672, 254, 224,
> 225, 499, 224, 224, 692, NA, 599, 627, 642, 677, 239, 224, NA,
> 239, 239, NA, 375, 239, 375, 224, 239, 239, NA, 299, 239, 224,
> 239, 224, 239, 621, 642, 239, 224, 239, 638, 239, 239, 239, 225,
> 224, 299, 672, 239, 618, 224, 620, NA, 626, 239, 657, 693, 224,
> 624, 693, 224, 225, 621, 657, 224, 239, 673, 239, 239, 239, 569,
> 224, 239, 239, NA, 657, 651)), .Names = c("X1", "X2"), row.names = c(NA,
> -100L), class = "data.frame")
>>
>
> You can see there are 3 distinct linear relationships of those 2 variables.
>
> plot(1/temp[,1], temp[,2])
>
> Is there any simple way how to evaluate such data without grouping variable? 
> I know that in case I have proper grouping variable I can evaluate it with 
> lme and get intercepts and/or slopes.
>
> My question is:
>
> Does anybody know about a way/package/function which can give me appropriate 
> grouping of such data or which can give me separate slope/intercept for each 
> set.
>
> I hope I expressed my problem clearly.
>
> Best regards
> Petr
>
>
> 
> Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou 
> určeny pouze jeho adresátům.
> Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně 
> jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze 
> svého systému.
> Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email 
> jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
> Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či 
> zpožděním přenosu e-mailu.
>
> V případě, že je tento e-mail součástí obchodního jednání:
> - vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, 
> a to z jakéhokoliv důvodu i bez uvedení důvodu.
> - a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; 
> Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany 
> příjemce s dodatkem či odchylkou.
> - trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným 
> dosažením shody na všech jejích náležitostech.
> - odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost 
> žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn nebo písemně 
> pověřen a takové pověření nebo plná moc byly adresátovi tohoto emailu 
> případně osobě, kterou adresát zastupuje, předloženy nebo jejich existence je 
> adresátovi či osobě jím zastoupené známá.
>
> This e-mail and any documents attached to it may be confidential and are 
> intended only for its intended recipients.
> If you received this e-mail by mistake, please immediately inform its sender. 
> Delete the contents of this e-mail with all attachments and its copies from 
> your system.
> If you are not the intended recipient of this e-mail, you are not authorized 
> to use, disseminate, copy or disclose this e-mail in any manner.
> The sender of this e-mail shall not be liable for any possible damage caused 
> by modifications of the e-mail or by delay with transfer of the email.
>
> In case that this e-mail forms part of business dealings:
> - the sender reserves the right to end negotiations about entering into a 
> contract in any time, for any reason, and without stating any reasoning.
> - if the e-mail contains an offer, the recipient is entitled to immediately 
> accept such offer; The sender of this e-mail (offer) excludes any acceptance 
> of the offer on the part of the recipient containing any amendment or 
> variation.
> - the sender insists on that the respective contract is concluded only upon 
> an express mutual agreement on all its aspects.
> - the sender of this e-mail informs that he/she is not authorized to enter 
> into any contracts on behalf of the company except for cases in which he/she 
> is expressly authorized to do so in writing, and such authorization or power 
> of attorney is submitted to the recipient or the person represented by the 
> recipient, or the existence of such authorization is known to the recipient 
> of the person represented by the recipient.
> __
> 

Re: [R] fancy linear model and grouping

2016-02-02 Thread Gabor Grothendieck
Try the EEV model with 3 clusters where temp is the large dataset:

   Mclust(temp, 3, modelNames = "EEV")

On Tue, Feb 2, 2016 at 8:13 AM, PIKAL Petr <petr.pi...@precheza.cz> wrote:
> Hi
>
> Thanks, it work for my example, which is actually a subset of a bigger data 
> (4000 rows) with the same characteristics. For the whole problem it does not 
> give correct clustering. I tried to set G to 3 but it did not help either.
>
> I attached the whole dataset (dput) that you can use, however after quick 
> tour through Mclust it seems to me that it is designed for slightly different 
> problem.
>
> Here is the result with whole data.
>
> fm <- Mclust(temp)
> g <- fm$classification
> plot(1/temp[,1], temp[,2], pch = g, col = g)
>
> I will go through the docs more thoroughly, to be 100% sure I did not miss 
> anything.
>
> Cheers
> Petr
>
>
>> -Original Message-
>> From: Gabor Grothendieck [mailto:ggrothendi...@gmail.com]
>> Sent: Tuesday, February 02, 2016 1:20 PM
>> To: PIKAL Petr
>> Cc: R Help R
>> Subject: Re: [R] fancy linear model and grouping
>>
>> Try the mclust package:
>>
>> library(mclust)
>> temp.na <- na.omit(temp)
>> fm <- Mclust(temp.na)
>> g <- fm$classification
>> plot(temp.na, pch = g, col = g)
>>
>>
>>
>> On Tue, Feb 2, 2016 at 6:35 AM, PIKAL Petr <petr.pi...@precheza.cz>
>> wrote:
>> > Dear all
>> >
>> > I have data like this
>> >
>> >> dput(temp)
>> >
>> > temp <- structure(list(X1 = c(93, 82, NA, 93, 93, 79, 79, 93, 93, 85,
>> > 82, 93, 87, 93, 92, NA, 87, 93, 93, 93, 74, 77, 87, 93, 82, 87, 75,
>> > 82, 93, 92, 68, 93, 93, 73, NA, 85, 81, 79, 75, 87, 93, NA, 87, 87,
>> > 85, 92, 87, 92, 93, 87, 87, NA, 69, 87, 93, 87, 93, 87, 82, 79, 87,
>> > 93, 87, 80, 87, 87, 87, 92, 93, 69, 76, 87, 82, 93, 82, NA, 54, 87,
>> > 77, 73, 93, 82, 73, 93, 92, 82, 77, 93, 87, 75, 87, 87, 87, 60, 92,
>> > 87, 87, NA, 77, 78), X2 = c(224, 624, NA, 224, 224, 642, 642, 224,
>> > 224, 599, 622, 224, 239, 224, 225, NA, 239, 224, 224, 224, 688, 657,
>> > 239, 224, 624, 239, 672, 254, 224, 225, 499, 224, 224, 692, NA, 599,
>> > 627, 642, 677, 239, 224, NA, 239, 239, NA, 375, 239, 375, 224, 239,
>> > 239, NA, 299, 239, 224, 239, 224, 239, 621, 642, 239, 224, 239, 638,
>> > 239, 239, 239, 225, 224, 299, 672, 239, 618, 224, 620, NA, 626, 239,
>> > 657, 693, 224, 624, 693, 224, 225, 621, 657, 224, 239, 673, 239, 239,
>> > 239, 569, 224, 239, 239, NA, 657, 651)), .Names = c("X1", "X2"),
>> > row.names = c(NA, -100L), class = "data.frame")
>> >>
>> >
>> > You can see there are 3 distinct linear relationships of those 2
>> variables.
>> >
>> > plot(1/temp[,1], temp[,2])
>> >
>> > Is there any simple way how to evaluate such data without grouping
>> variable? I know that in case I have proper grouping variable I can
>> evaluate it with lme and get intercepts and/or slopes.
>> >
>> > My question is:
>> >
>> > Does anybody know about a way/package/function which can give me
>> appropriate grouping of such data or which can give me separate
>> slope/intercept for each set.
>> >
>> > I hope I expressed my problem clearly.
>> >
>> > Best regards
>> > Petr
>> >
>> >
>
> 
> Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou 
> určeny pouze jeho adresátům.
> Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně 
> jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze 
> svého systému.
> Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email 
> jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
> Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či 
> zpožděním přenosu e-mailu.
>
> V případě, že je tento e-mail součástí obchodního jednání:
> - vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, 
> a to z jakéhokoliv důvodu i bez uvedení důvodu.
> - a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; 
> Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany 
> příjemce s dodatkem či odchylkou.
> - trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným 
> dosažením shody na všech jejích náležitostech.
> - odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost 
> žádné smlouvy s výjimkou pří

Re: [R] Optim() and Instability

2015-11-16 Thread Gabor Grothendieck
Since some questioned the scaling idea, here are runs first with
scaling and then without scaling.  Note how much better the solution
is in the first run (see arrows).  It is also evident from the data

> head(data, 3)
  y x1   x2   x3
1 0.660 20  7.0 1680
2 0.165  5  1.7  350
3 0.660 20  7.0 1680

> # run 1 - scaling
> str(optim(par = c(1,1, 1), min.perc_error, data = data,
+   control = list(parscale = c(1, 1, 0.0001
List of 5
 $ par: num [1:3] 0.030232 0.024411 -0.000113
 $ value  : num 0.653  <=
 $ counts : Named int [1:2] 180 NA
  ..- attr(*, "names")= chr [1:2] "function" "gradient"
 $ convergence: int 0
 $ message: NULL

> # run 2 - no scaling
> str(optim(par = c(1,1, 1), min.perc_error, data = data))
List of 5
 $ par: num [1:3] 0.6305 -0.1247 -0.0032
 $ value  : num 473  <=
 $ counts : Named int [1:2] 182 NA
  ..- attr(*, "names")= chr [1:2] "function" "gradient"
 $ convergence: int 0
 $ message: NULL

On Sat, Nov 14, 2015 at 10:32 AM, Gabor Grothendieck
<ggrothendi...@gmail.com> wrote:
> I meant the parscale parameter.
>
> On Sat, Nov 14, 2015 at 10:30 AM, Gabor Grothendieck
> <ggrothendi...@gmail.com> wrote:
>> Tyipcally the parameters being optimized should be the same order of
>> magnitude or else you can expect numerical problems.  That is what the
>> fnscale control parameter is for.
>>
>> On Sat, Nov 14, 2015 at 10:15 AM, Lorenzo Isella
>> <lorenzo.ise...@gmail.com> wrote:
>>> Dear All,
>>> I am using optim() for a relatively simple task: a linear model where
>>> instead of minimizing the sum of the squared errors, I minimize the sum
>>> of the squared relative errors.
>>> However, I notice that the default algorithm is very sensitive to the
>>> choice of the initial fit parameters, whereas I get much more stable
>>> (and therefore better?) results with the BFGS algorithm.
>>> I would like to have some feedback on this (perhaps I made a mistake
>>> somewhere).
>>> I provide a small self-contained example.
>>> You can download a tiny data set from the link
>>>
>>> https://www.dropbox.com/s/tmbj3os4ev3d4y8/data-instability.csv?dl=0
>>>
>>> whereas I paste the script I am using at the end of the email.
>>> Any feedback is really appreciated.
>>> Many thanks
>>>
>>> Lorenzo
>>>
>>> 
>>>
>>> min.perc_error <- function(data, par) {
>>>  with(data, sum(((par[1]*x1 + par[2]*x2+par[3]*x3 -
>>>  y)/y)^2))
>>>}
>>>
>>> par_ini1 <- c(.3,.1, 1e-3)
>>>
>>> par_ini2 <- c(1,1, 1)
>>>
>>>
>>> data <- read.csv("data-instability.csv")
>>>
>>> mm_def1 <-optim(par = par_ini1
>>>, min.perc_error, data = data)
>>>
>>> mm_bfgs1 <-optim(par = par_ini1
>>>, min.perc_error, data = data, method="BFGS")
>>>
>>> print("fit parameters with the default algorithms and the first seed
>>> ")
>>> print(mm_def1$par)
>>>
>>> print("fit parameters with the BFGS algorithms and the first seed  ")
>>> print(mm_bfgs1$par)
>>>
>>>
>>>
>>> mm_def2 <-optim(par = par_ini2
>>>, min.perc_error, data = data)
>>>
>>> mm_bfgs2 <-optim(par = par_ini2
>>>, min.perc_error, data = data, method="BFGS")
>>>
>>>
>>>
>>>
>>> print("fit parameters with the default algorithms and the second seed
>>> ")
>>> print(mm_def2$par)
>>>
>>> print("fit parameters with the BFGS algorithms and the second seed  ")
>>> print(mm_bfgs2$par)
>>>
>>> __
>>> 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.
>>
>>
>>
>> --
>> Statistics & Software Consulting
>> GKX Group, GKX Associates Inc.
>> tel: 1-877-GKX-GROUP
>> email: ggrothendieck at gmail.com
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Optim() and Instability

2015-11-14 Thread Gabor Grothendieck
Tyipcally the parameters being optimized should be the same order of
magnitude or else you can expect numerical problems.  That is what the
fnscale control parameter is for.

On Sat, Nov 14, 2015 at 10:15 AM, Lorenzo Isella
 wrote:
> Dear All,
> I am using optim() for a relatively simple task: a linear model where
> instead of minimizing the sum of the squared errors, I minimize the sum
> of the squared relative errors.
> However, I notice that the default algorithm is very sensitive to the
> choice of the initial fit parameters, whereas I get much more stable
> (and therefore better?) results with the BFGS algorithm.
> I would like to have some feedback on this (perhaps I made a mistake
> somewhere).
> I provide a small self-contained example.
> You can download a tiny data set from the link
>
> https://www.dropbox.com/s/tmbj3os4ev3d4y8/data-instability.csv?dl=0
>
> whereas I paste the script I am using at the end of the email.
> Any feedback is really appreciated.
> Many thanks
>
> Lorenzo
>
> 
>
> min.perc_error <- function(data, par) {
>  with(data, sum(((par[1]*x1 + par[2]*x2+par[3]*x3 -
>  y)/y)^2))
>}
>
> par_ini1 <- c(.3,.1, 1e-3)
>
> par_ini2 <- c(1,1, 1)
>
>
> data <- read.csv("data-instability.csv")
>
> mm_def1 <-optim(par = par_ini1
>, min.perc_error, data = data)
>
> mm_bfgs1 <-optim(par = par_ini1
>, min.perc_error, data = data, method="BFGS")
>
> print("fit parameters with the default algorithms and the first seed
> ")
> print(mm_def1$par)
>
> print("fit parameters with the BFGS algorithms and the first seed  ")
> print(mm_bfgs1$par)
>
>
>
> mm_def2 <-optim(par = par_ini2
>, min.perc_error, data = data)
>
> mm_bfgs2 <-optim(par = par_ini2
>, min.perc_error, data = data, method="BFGS")
>
>
>
>
> print("fit parameters with the default algorithms and the second seed
> ")
> print(mm_def2$par)
>
> print("fit parameters with the BFGS algorithms and the second seed  ")
> print(mm_bfgs2$par)
>
> __
> 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.



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Optim() and Instability

2015-11-14 Thread Gabor Grothendieck
I meant the parscale parameter.

On Sat, Nov 14, 2015 at 10:30 AM, Gabor Grothendieck
<ggrothendi...@gmail.com> wrote:
> Tyipcally the parameters being optimized should be the same order of
> magnitude or else you can expect numerical problems.  That is what the
> fnscale control parameter is for.
>
> On Sat, Nov 14, 2015 at 10:15 AM, Lorenzo Isella
> <lorenzo.ise...@gmail.com> wrote:
>> Dear All,
>> I am using optim() for a relatively simple task: a linear model where
>> instead of minimizing the sum of the squared errors, I minimize the sum
>> of the squared relative errors.
>> However, I notice that the default algorithm is very sensitive to the
>> choice of the initial fit parameters, whereas I get much more stable
>> (and therefore better?) results with the BFGS algorithm.
>> I would like to have some feedback on this (perhaps I made a mistake
>> somewhere).
>> I provide a small self-contained example.
>> You can download a tiny data set from the link
>>
>> https://www.dropbox.com/s/tmbj3os4ev3d4y8/data-instability.csv?dl=0
>>
>> whereas I paste the script I am using at the end of the email.
>> Any feedback is really appreciated.
>> Many thanks
>>
>> Lorenzo
>>
>> 
>>
>> min.perc_error <- function(data, par) {
>>  with(data, sum(((par[1]*x1 + par[2]*x2+par[3]*x3 -
>>  y)/y)^2))
>>}
>>
>> par_ini1 <- c(.3,.1, 1e-3)
>>
>> par_ini2 <- c(1,1, 1)
>>
>>
>> data <- read.csv("data-instability.csv")
>>
>> mm_def1 <-optim(par = par_ini1
>>, min.perc_error, data = data)
>>
>> mm_bfgs1 <-optim(par = par_ini1
>>, min.perc_error, data = data, method="BFGS")
>>
>> print("fit parameters with the default algorithms and the first seed
>> ")
>> print(mm_def1$par)
>>
>> print("fit parameters with the BFGS algorithms and the first seed  ")
>> print(mm_bfgs1$par)
>>
>>
>>
>> mm_def2 <-optim(par = par_ini2
>>, min.perc_error, data = data)
>>
>> mm_bfgs2 <-optim(par = par_ini2
>>, min.perc_error, data = data, method="BFGS")
>>
>>
>>
>>
>> print("fit parameters with the default algorithms and the second seed
>> ")
>> print(mm_def2$par)
>>
>> print("fit parameters with the BFGS algorithms and the second seed  ")
>> print(mm_bfgs2$par)
>>
>> __
>> 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.
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Linear regression with a rounded response variable

2015-10-21 Thread Gabor Grothendieck
This could be modeled directly using Bayesian techniques. Consider the
Bayesian version of the following model where we only observe y and X.  y0
is not observed.

   y0 <- X b + error
   y <- round(y0)

The following code is based on modifying the code in the README of the CRAN
rcppbugs R package.


library(rcppbugs)
set.seed(123)

# set up the test data - y and X are observed but not y0
NR <- 1e2L
NC <- 2L
X <- cbind(1, rnorm(10))
y0 <- X %*% 1:2
y <- round(y0)

# for comparison run a normal linear model w/ lm.fit using X and y
lm.res <- lm.fit(X,y)
print(coef(lm.res))
##x1x2
## 0.9569366 1.9170808

# RCppBugs Model
b <- mcmc.normal(rnorm(NC),mu=0,tau=0.0001)
tau.y <- mcmc.gamma(sd(as.vector(y)),alpha=0.1,beta=0.1)
y.hat <- deterministic(function(X,b) { round(X %*% b) }, X, b)
y.lik <- mcmc.normal(y,mu=y.hat,tau=tau.y,observed=TRUE)
m <- create.model(b, tau.y, y.hat, y.lik)

# run the Bayesian model based on y and X
cat("running model...\n")
runtime <- system.time(ans <- run.model(m, iterations=1e5L, burn=1e4L,
adapt=1e3L, thin=10L))
print(apply(ans[["b"]],2,mean))
## [1] 0.9882485 2.0009989


On Wed, Oct 21, 2015 at 10:53 AM, Ravi Varadhan 
wrote:

> Hi,
> I am dealing with a regression problem where the response variable, time
> (second) to walk 15 ft, is rounded to the nearest integer.  I do not care
> for the regression coefficients per se, but my main interest is in getting
> the prediction equation for walking speed, given the predictors (age,
> height, sex, etc.), where the predictions will be real numbers, and not
> integers.  The hope is that these predictions should provide unbiased
> estimates of the "unrounded" walking speed. These sounds like a measurement
> error problem, where the measurement error is due to rounding and hence
> would be uniformly distributed (-0.5, 0.5).
>
> Are there any canonical approaches for handling this type of a problem?
> What is wrong with just doing the standard linear regression?
>
> I googled and saw that this question was asked by someone else in a
> stackexchange post, but it was unanswered.  Any suggestions?
>
> Thank you,
> Ravi
>
> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
> Associate Professor,  Department of Oncology
> Division of Biostatistics & Bionformatics
> Sidney Kimmel Comprehensive Cancer Center
> Johns Hopkins University
> 550 N. Broadway, Suite -E
> Baltimore, MD 21205
> 410-502-2619
>
>
> [[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.
>



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

[[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.


Re: [R] Most appropriate function for the following optimisation issue?

2015-10-20 Thread Gabor Grothendieck
Correction.

Yes, it's the projection of S onto the subspace orthogonal to B which is:

X <- S - (B%o%B) %*% S/ sum(B*B)

and is also implied by Duncan's solution since that is what the residuals
of linear regression are.

On Tue, Oct 20, 2015 at 1:11 PM, Gabor Grothendieck <ggrothendi...@gmail.com
> wrote:

> Yes, it's the projection of S onto the subspace orthogonal to B which is:
>
> X <- S - B%*%B / sum(B*B)
>
> and is also implied by Duncan's solution since that is what the residuals
> of linear regression are.
>
> On Tue, Oct 20, 2015 at 1:00 PM, Paul Smith <phh...@gmail.com> wrote:
>
>> On Tue, Oct 20, 2015 at 11:58 AM, Andy Yuan <yuan...@gmail.com> wrote:
>> >
>> > Please could you help me to select the most appropriate/fastest
>> function to use for the following constraint optimisation issue?
>> >
>> > Objective function:
>> >
>> > Min: Sum( (X[i] - S[i] )^2)
>> >
>> > Subject to constraint :
>> >
>> > Sum (B[i] x X[i]) =0
>> >
>> > where i=1…n and S[i] and B[i] are real numbers
>> >
>> > Need to solve for X
>> >
>> > Example:
>> >
>> > Assume n=3
>> >
>> > S <- c(-0.5, 7.8, 2.3)
>> > B <- c(0.42, 1.12, 0.78)
>> >
>> > Many thanks
>>
>> I believe you can solve *analytically* your optimization problem, with
>> the Lagrange multipliers method, Andy. By doing so, you can derive
>> clean and closed-form expression for the optimal solution.
>>
>> Paul
>>
>> __
>> 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.
>>
>
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com
>



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

[[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.

Re: [R] Most appropriate function for the following optimisation issue?

2015-10-20 Thread Gabor Grothendieck
Yes, it's the projection of S onto the subspace orthogonal to B which is:

X <- S - B%*%B / sum(B*B)

and is also implied by Duncan's solution since that is what the residuals
of linear regression are.

On Tue, Oct 20, 2015 at 1:00 PM, Paul Smith  wrote:

> On Tue, Oct 20, 2015 at 11:58 AM, Andy Yuan  wrote:
> >
> > Please could you help me to select the most appropriate/fastest function
> to use for the following constraint optimisation issue?
> >
> > Objective function:
> >
> > Min: Sum( (X[i] - S[i] )^2)
> >
> > Subject to constraint :
> >
> > Sum (B[i] x X[i]) =0
> >
> > where i=1…n and S[i] and B[i] are real numbers
> >
> > Need to solve for X
> >
> > Example:
> >
> > Assume n=3
> >
> > S <- c(-0.5, 7.8, 2.3)
> > B <- c(0.42, 1.12, 0.78)
> >
> > Many thanks
>
> I believe you can solve *analytically* your optimization problem, with
> the Lagrange multipliers method, Andy. By doing so, you can derive
> clean and closed-form expression for the optimal solution.
>
> Paul
>
> __
> 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.
>



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

[[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.

Re: [R] flatten a list

2015-09-29 Thread Gabor Grothendieck
> do.call(c, lapply(temp, function(x) if (is.list(x)) x else list(x)))
[[1]]
[1] 1 2 3

[[2]]
[1] "a" "b" "c"

$duh
[1] 5 6 7 8

$zed
[1] 15 16 17


On Tue, Sep 29, 2015 at 11:00 AM, Therneau, Terry M., Ph.D. <
thern...@mayo.edu> wrote:

> I'd like to flatten a list from 2 levels to 1 level.  This has to be easy,
> but is currently opaque to me.
>
> temp <- list(1:3, list(letters[1:3], duh= 5:8),  zed=15:17)
>
> Desired result would be a 4 element list.
> [[1]] 1:3
> [[2]] "a", "b", "c"
> [[duh]] 5:8
> [[zed]] 15:17
>
> (Preservation of the names is not important)
>
> Terry T
> __
> 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.
>



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

[[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.


Re: [R] How to coerce a parameter in nls?

2015-09-22 Thread Gabor Grothendieck
You may have to do without masking and switch back to nls.  dproot2 and fo
are from prior post.

# to mask Rm6 omit it from start and set it explicitly
st <- c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65, Rm5=1.01, d50=20, c=-1)
Rm6 <- 1

fm.nls <- nls(fo, dproot2, start = st)

AIC(fm.nls)
summary(fm.nls)


On Tue, Sep 22, 2015 at 12:46 PM, Jianling Fan 
wrote:

> Hello Prof. Nash,
>
> My regression works good now. But I found another problem when I using
> nlxb. In the output, the SE, t-stat, and p-value are not available.
> Furthermore, I can't extract AIC from the output. The output looks
> like below:
>
> Do you have any suggestion for this?
>
> Thanks a lot!
>
> Regards,
>
> nlmrt class object: x
> residual sumsquares =  0.29371  on  33 observations
> after  9Jacobian and  10 function evaluations
>   namecoeff  SE   tstat  pval
> gradientJSingval
> Rm1   1.1162NA NA NA
> -3.059e-13   2.745
> Rm2  1.56072NA NA NA
> 1.417e-131.76
> Rm3  1.09775NA NA NA
> -3.179e-13   1.748
> Rm4  7.18377NA NA NA
> -2.941e-12   1.748
> Rm5  1.13562NA NA NA
> -3.305e-13   1.076
> Rm61  M NA NA NA
> 0   0.603
> d50  22.4803NA NA NA
> 4.975e-13   0.117
> c   -1.64075NA NA NA
> 4.12e-12   1.908e-17
>
>
>
> On 21 September 2015 at 13:38, ProfJCNash  wrote:
> > I've not used it for group data, and suspect that the code to generate
> > derivatives cannot cope with the bracket syntax. If you can rewrite the
> > equation without the brackets, you could get the derivatives and solve
> that
> > way. This will probably mean having a "translation" routine to glue
> things
> > together.
> >
> > JN
> >
> >
> > On 15-09-21 12:22 PM, Jianling Fan wrote:
> >>
> >> Thanks Prof. Nash,
> >>
> >> Sorry for late reply. I am learning and trying to use your nlmrt
> >> package since I got your email. It works good to mask a parameter in
> >> regression but seems does work for my equation. I think the problem is
> >> that the parameter I want to mask is a group-specific parameter and I
> >> have a "[]" syntax in my equation. However, I don't have your 2014
> >> book on hand and couldn't find it in our library. So I am wondering if
> >> nlxb works for group data?
> >> Thanks a lot!
> >>
> >> following is my code and I got a error form it.
> >>
> >>> fitdp1<-nlxb(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,
> >>
> >>  + start =c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
> >> Rm5=1.01, Rm6=1, d50=20, c=-1),
> >>  + masked=c("Rm6"))
> >>
> >> Error in deriv.default(parse(text = resexp), names(start)) :
> >>Function '`[`' is not in the derivatives table
> >>
> >>
> >> Best regards,
> >>
> >> Jianling
> >>
> >>
> >> On 20 September 2015 at 12:56, ProfJCNash  wrote:
> >>>
> >>> I posted a suggestion to use nlmrt package (function nlxb to be
> precise),
> >>> which has masked (fixed) parameters. Examples in my 2014 book on
> >>> Nonlinear
> >>> parameter optimization with R tools. However, I'm travelling just now,
> or
> >>> would consider giving this a try.
> >>>
> >>> JN
> >>>
> >>>
> >>> On 15-09-20 01:19 PM, Jianling Fan wrote:
> 
> 
>  no, I am doing a regression with 6 group data with 2 shared parameters
>  and 1 different parameter for each group data. the parameter I want to
>  coerce is for one group. I don't know how to do it. Any suggestion?
> 
>  Thanks!
> 
>  On 19 September 2015 at 13:33, Jeff Newmiller <
> jdnew...@dcn.davis.ca.us>
>  wrote:
> >
> >
> > Why not rewrite the function so that value is not a parameter?
> >
> >
> >
> ---
> > Jeff NewmillerThe .   .  Go
> > Live...
> > DCN:Basics: ##.#.   ##.#.
> 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.
> >
> > On September 18, 2015 9:54:54 PM PDT, Jianling Fan
> >  wrote:
> >>
> >>
> >> Hello, everyone,
> >>
> >> I am using a nls regression with 6 groups data. I am trying to
> coerce
> >> a parameter to 1 by using a upper and lower statement. but I 

Re: [R] How to coerce a parameter in nls?

2015-09-22 Thread Gabor Grothendieck
Just write out the 20 terms.

On Mon, Sep 21, 2015 at 10:26 PM, Jianling Fan <fanjianl...@gmail.com>
wrote:

> Hello, Gabor,
>
> Thanks again for your suggestion. And now I am trying to improve the
> code by adding a function to replace the express "Rm1 * ref.1 + Rm2 *
> ref.2 + Rm3 * ref.3 + Rm4 * ref.4 + Rm5 * ref.5 + Rm6 * ref.6" because
> I have some other dataset need to fitted to the same model but with
> more groups (>20).
>
> I tried to add the function as:
>
> denfun<-function(i){
>for(i in 1:6){
>  Rm<-sum(Rm[i]*ref.i)
>  return(Rm)}
> }
>
> but I got another error when I incorporate this function into my
> regression:
>
> >fitdp1<-nlxb(den ~ denfun(6)/(1+(depth/d50)^c),
>data = dproot2,
>  start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
> Rm5=1.01, Rm6=1, d50=20, c=-1),
> masked = "Rm6")
>
> Error in deriv.default(parse(text = resexp), names(start)) :
>   Function 'denfun' is not in the derivatives table
>
> I think there must be something wrong with my function. I tried some
> times but am not sure how to improve it because I am quite new to R.
>
> Could anyone please give me some suggestion.
>
> Thanks a lot!
>
>
> Jianling
>
>
> On 22 September 2015 at 00:43, Gabor Grothendieck
> <ggrothendi...@gmail.com> wrote:
> > Express the formula in terms of simple operations like this:
> >
> > # add 0/1 columns ref.1, ref.2, ..., ref.6
> > dproot2 <- do.call(data.frame, transform(dproot, ref = outer(dproot$ref,
> > seq(6), "==") + 0))
> >
> > # now express the formula in terms of the new columns
> > library(nlmrt)
> > fitdp1<-nlxb(den ~ (Rm1 * ref.1 + Rm2 * ref.2 + Rm3 * ref.3 + Rm4 *
> ref.4 +
> > Rm5 * ref.5 + Rm6 * ref.6)/(1+(depth/d50)^c),
> >  data = dproot2,
> >  start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65, Rm5=1.01,
> Rm6=1,
> > d50=20, c=-1),
> >  masked = "Rm6")
> >
> > where we used this input:
> >
> > Lines <- "   depth   den ref
> > 1 20 0.573   1
> > 2 40 0.780   1
> > 3 60 0.947   1
> > 4 80 0.990   1
> > 5100 1.000   1
> > 6 10 0.600   2
> > 7 20 0.820   2
> > 8 30 0.930   2
> > 9 40 1.000   2
> > 1020 0.480   3
> > 1140 0.734   3
> > 1260 0.961   3
> > 1380 0.998   3
> > 14   100 1.000   3
> > 1520 3.2083491   4
> > 1640 4.9683383   4
> > 1760 6.2381133   4
> > 1880 6.5322348   4
> > 19   100 6.5780660   4
> > 20   120 6.6032064   4
> > 2120 0.614   5
> > 2240 0.827   5
> > 2360 0.950   5
> > 2480 0.995   5
> > 25   100 1.000   5
> > 2620 0.4345774   6
> > 2740 0.6654726   6
> > 2860 0.8480684   6
> > 2980 0.9268951   6
> > 30   100 0.9723207   6
> > 31   120 0.9939966   6
> > 32   140 0.9992400   6"
> >
> > dproot <- read.table(text = Lines, header = TRUE)
> >
> >
> >
> > On Mon, Sep 21, 2015 at 12:22 PM, Jianling Fan <fanjianl...@gmail.com>
> > wrote:
> >>
> >> Thanks Prof. Nash,
> >>
> >> Sorry for late reply. I am learning and trying to use your nlmrt
> >> package since I got your email. It works good to mask a parameter in
> >> regression but seems does work for my equation. I think the problem is
> >> that the parameter I want to mask is a group-specific parameter and I
> >> have a "[]" syntax in my equation. However, I don't have your 2014
> >> book on hand and couldn't find it in our library. So I am wondering if
> >> nlxb works for group data?
> >> Thanks a lot!
> >>
> >> following is my code and I got a error form it.
> >>
> >> > fitdp1<-nlxb(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,
> >> + start =c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
> >> Rm5=1.01, Rm6=1, d50=20, c=-1),
> >> + masked=c("Rm6"))
> >>
> >> Error in deriv.default(parse(text = resexp), names(start)) :
> >>   Function '`[`' is not in the derivatives table
> >>
> >>
> >> Best regards,
> >>
> >> Jianling
> >>
> >>
> >> On 20 September 2015 at 12:56, ProfJCNash <profjcn...@gmail.com

Re: [R] How to coerce a parameter in nls?

2015-09-22 Thread Gabor Grothendieck
Or if you really can't bear to write out 20 terms have R do it for you:

# number of terms is the number of unique values in ref column
nterms <- length(unique(dproot$ref))

dproot2 <- do.call(data.frame, transform(dproot, ref =
outer(dproot$ref, seq(nterms),
"==") + 0))

# construct the formula as a string
terms <- paste( sprintf("Rm%d*ref.%d", 1:nterms, 1:nterms), collapse = "+")
fo <- sprintf("den ~ (%s)/(1+(depth/d50)^c)", terms)

library(nlmrt)
fm <- nlxb(fo, data = dproot2, masked = "Rm6",
 start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65, Rm5=1.01, Rm6=1,
d50=20, c=-1))


On Tue, Sep 22, 2015 at 7:04 AM, Gabor Grothendieck <ggrothendi...@gmail.com
> wrote:

> Just write out the 20 terms.
>
> On Mon, Sep 21, 2015 at 10:26 PM, Jianling Fan <fanjianl...@gmail.com>
> wrote:
>
>> Hello, Gabor,
>>
>> Thanks again for your suggestion. And now I am trying to improve the
>> code by adding a function to replace the express "Rm1 * ref.1 + Rm2 *
>> ref.2 + Rm3 * ref.3 + Rm4 * ref.4 + Rm5 * ref.5 + Rm6 * ref.6" because
>> I have some other dataset need to fitted to the same model but with
>> more groups (>20).
>>
>> I tried to add the function as:
>>
>> denfun<-function(i){
>>for(i in 1:6){
>>  Rm<-sum(Rm[i]*ref.i)
>>  return(Rm)}
>> }
>>
>> but I got another error when I incorporate this function into my
>> regression:
>>
>> >fitdp1<-nlxb(den ~ denfun(6)/(1+(depth/d50)^c),
>>data = dproot2,
>>  start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
>> Rm5=1.01, Rm6=1, d50=20, c=-1),
>> masked = "Rm6")
>>
>> Error in deriv.default(parse(text = resexp), names(start)) :
>>   Function 'denfun' is not in the derivatives table
>>
>> I think there must be something wrong with my function. I tried some
>> times but am not sure how to improve it because I am quite new to R.
>>
>> Could anyone please give me some suggestion.
>>
>> Thanks a lot!
>>
>>
>> Jianling
>>
>>
>> On 22 September 2015 at 00:43, Gabor Grothendieck
>> <ggrothendi...@gmail.com> wrote:
>> > Express the formula in terms of simple operations like this:
>> >
>> > # add 0/1 columns ref.1, ref.2, ..., ref.6
>> > dproot2 <- do.call(data.frame, transform(dproot, ref = outer(dproot$ref,
>> > seq(6), "==") + 0))
>> >
>> > # now express the formula in terms of the new columns
>> > library(nlmrt)
>> > fitdp1<-nlxb(den ~ (Rm1 * ref.1 + Rm2 * ref.2 + Rm3 * ref.3 + Rm4 *
>> ref.4 +
>> > Rm5 * ref.5 + Rm6 * ref.6)/(1+(depth/d50)^c),
>> >  data = dproot2,
>> >  start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65, Rm5=1.01,
>> Rm6=1,
>> > d50=20, c=-1),
>> >  masked = "Rm6")
>> >
>> > where we used this input:
>> >
>> > Lines <- "   depth   den ref
>> > 1 20 0.573   1
>> > 2 40 0.780   1
>> > 3 60 0.947   1
>> > 4 80 0.990   1
>> > 5100 1.000   1
>> > 6 10 0.600   2
>> > 7 20 0.820   2
>> > 8 30 0.930   2
>> > 9 40 1.000   2
>> > 1020 0.480   3
>> > 1140 0.734   3
>> > 1260 0.961   3
>> > 1380 0.998   3
>> > 14   100 1.000   3
>> > 1520 3.2083491   4
>> > 1640 4.9683383   4
>> > 1760 6.2381133   4
>> > 1880 6.5322348   4
>> > 19   100 6.5780660   4
>> > 20   120 6.6032064   4
>> > 2120 0.614   5
>> > 2240 0.827   5
>> > 2360 0.950   5
>> > 2480 0.995   5
>> > 25   100 1.000   5
>> > 2620 0.4345774   6
>> > 2740 0.6654726   6
>> > 2860 0.8480684   6
>> > 2980 0.9268951   6
>> > 30   100 0.9723207   6
>> > 31   120 0.9939966   6
>> > 32   140 0.9992400   6"
>> >
>> > dproot <- read.table(text = Lines, header = TRUE)
>> >
>> >
>> >
>> > On Mon, Sep 21, 2015 at 12:22 PM, Jianling Fan <fanjianl...@gmail.com>
>> > wrote:
>> >>
>> >> Thanks Prof. Nash,
>> >>
>> >> Sorry for late reply. I am learning and trying to use your nlmrt
>> >> package since I got your email. It

Re: [R] How to coerce a parameter in nls?

2015-09-21 Thread Gabor Grothendieck
Express the formula in terms of simple operations like this:

# add 0/1 columns ref.1, ref.2, ..., ref.6
dproot2 <- do.call(data.frame, transform(dproot, ref = outer(dproot$ref,
seq(6), "==") + 0))

# now express the formula in terms of the new columns
library(nlmrt)
fitdp1<-nlxb(den ~ (Rm1 * ref.1 + Rm2 * ref.2 + Rm3 * ref.3 + Rm4 * ref.4 +
Rm5 * ref.5 + Rm6 * ref.6)/(1+(depth/d50)^c),
 data = dproot2,
 start = c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65, Rm5=1.01, Rm6=1,
d50=20, c=-1),
 masked = "Rm6")

where we used this input:

Lines <- "   depth   den ref
1 20 0.573   1
2 40 0.780   1
3 60 0.947   1
4 80 0.990   1
5100 1.000   1
6 10 0.600   2
7 20 0.820   2
8 30 0.930   2
9 40 1.000   2
1020 0.480   3
1140 0.734   3
1260 0.961   3
1380 0.998   3
14   100 1.000   3
1520 3.2083491   4
1640 4.9683383   4
1760 6.2381133   4
1880 6.5322348   4
19   100 6.5780660   4
20   120 6.6032064   4
2120 0.614   5
2240 0.827   5
2360 0.950   5
2480 0.995   5
25   100 1.000   5
2620 0.4345774   6
2740 0.6654726   6
2860 0.8480684   6
2980 0.9268951   6
30   100 0.9723207   6
31   120 0.9939966   6
32   140 0.9992400   6"

dproot <- read.table(text = Lines, header = TRUE)



On Mon, Sep 21, 2015 at 12:22 PM, Jianling Fan 
wrote:

> Thanks Prof. Nash,
>
> Sorry for late reply. I am learning and trying to use your nlmrt
> package since I got your email. It works good to mask a parameter in
> regression but seems does work for my equation. I think the problem is
> that the parameter I want to mask is a group-specific parameter and I
> have a "[]" syntax in my equation. However, I don't have your 2014
> book on hand and couldn't find it in our library. So I am wondering if
> nlxb works for group data?
> Thanks a lot!
>
> following is my code and I got a error form it.
>
> > fitdp1<-nlxb(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,
> + start =c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
> Rm5=1.01, Rm6=1, d50=20, c=-1),
> + masked=c("Rm6"))
>
> Error in deriv.default(parse(text = resexp), names(start)) :
>   Function '`[`' is not in the derivatives table
>
>
> Best regards,
>
> Jianling
>
>
> On 20 September 2015 at 12:56, ProfJCNash  wrote:
> > I posted a suggestion to use nlmrt package (function nlxb to be precise),
> > which has masked (fixed) parameters. Examples in my 2014 book on
> Nonlinear
> > parameter optimization with R tools. However, I'm travelling just now, or
> > would consider giving this a try.
> >
> > JN
> >
> >
> > On 15-09-20 01:19 PM, Jianling Fan wrote:
> >>
> >> no, I am doing a regression with 6 group data with 2 shared parameters
> >> and 1 different parameter for each group data. the parameter I want to
> >> coerce is for one group. I don't know how to do it. Any suggestion?
> >>
> >> Thanks!
> >>
> >> On 19 September 2015 at 13:33, Jeff Newmiller  >
> >> wrote:
> >>>
> >>> Why not rewrite the function so that value is not a parameter?
> >>>
> >>>
> ---
> >>> Jeff NewmillerThe .   .  Go
> >>> Live...
> >>> DCN:Basics: ##.#.   ##.#.  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.
> >>>
> >>> On September 18, 2015 9:54:54 PM PDT, Jianling Fan
> >>>  wrote:
> 
>  Hello, everyone,
> 
>  I am using a nls regression with 6 groups data. I am trying to coerce
>  a parameter to 1 by using a upper and lower statement. but I always
>  get an error like below:
> 
>  Error in ifelse(internalPars < upper, 1, -1) :
>    (list) object cannot be coerced to type 'double'
> 
>  does anyone know how to fix it?
> 
>  thanks in advance!
> 
>  My code is below:
> 
> 
> 
> > dproot
> 
> depth   den ref
>  1 20 0.573   1
>  2 40 0.780   1
>  3 60 0.947   1
>  4 80 0.990   1
>  5100 1.000   1
>  6 10 0.600   2
>  7 20 0.820   2
>  8 30 0.930   2
>  9 40 1.000   2
>  1020 0.480   3
>  1140 0.734   3
>  1260 0.961   3
>  1380 0.998   3
>  14   100 1.000   3
>  1520 3.2083491   4
>  1640 4.9683383   4
>  1760 6.2381133   4
>  18

Re: [R] How to read CSV from web?

2015-09-14 Thread Gabor Grothendieck
Your read.csv call works for me under Windows on "R version 3.2.2 Patched
(2015-08-25 r69180)"  but not on "R version 3.1.3 Patched (2015-03-16
r68169)".

Suggest you upgrade your R installation and try again.

If you are on Windowsw and don't want to upgrade right now an alternative
is to issue this command first: setInternet2()

Also note that header = TRUE and sep = "," are the defaults for read.csv so
those arguments can be optionally omitted.



On Wed, Jul 29, 2015 at 6:37 AM, jpara3 
wrote:

> data<-read.csv("
> https://raw.githubusercontent.com/sjkiss/Survey/master/mlogit.out.csv
> ",header=T,sep=",")
>
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/How-to-read-CSV-from-web-tp4710502p4710513.html
> Sent from the R help mailing list archive at Nabble.com.
>
> __
> 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.
>



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

[[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.


Re: [R] extracting every nth character from a string...

2015-09-06 Thread Gabor Grothendieck
This uses a regular expression but is shorter:

> gsub("(.).", "\\1", "ABCDEFG")
[1] "ACEG"

It replaces each successive pair of characters with the first of that
pair.  If there is an odd number of characters then the last character is
not matched and therefore kept -- thus it works properly for both even and
odd.


On Sat, Sep 5, 2015 at 4:59 PM, Evan Cooch  wrote:

> Suppose I had the following string, which has length of integer multiple
> of some value n. So, say n=2, and the example string has a length of  (2x4)
> = 8 characters.
>
> str <- "ABCDEFGH"
>
> What I'm trying to figure out is a simple, base-R coded way (which I
> heuristically call StrSubset in the following) to extract every nth
> character from the string, to generate a new string.
>
> So
>
> str <- "ABCDEFGH"
>
> new_str <- StrSubset(str);
>
> print(new_str)
>
> which would yield
>
> "ACEG"
>
>
> Best I could come up with is something like the following, where I extract
> every odd character from the string:
>
> StrSubset <- function(string)
>   {
> paste(unlist(strsplit(string,""))[seq(1,nchar(string),2)],collapse="") }
>
>
> Anything more elegant come to mind? Trying to avoid regex if possible
> (harder to explain to end-users), but if that meets the 'more elegant'
> sniff test, happy to consider...
>
> Thanks in advance...
>
> __
> 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.
>



-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

[[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.


Re: [R] : Ramanujan and the accuracy of floating point computations - using Rmpfr in R

2015-07-04 Thread Gabor Grothendieck
You can do that with bc if you pass the entire expression to bc(...) in
quotes but in that case you will have to use bc notation, not R notation
so, for example, exp is e and atan is a.

 library(bc)
 bc(e(sqrt(163)*4*a(1)))
[1]
262537412640768743.2500725971981856888793538563373369908627075374103782106479101186073116295306145602054347


On Fri, Jul 3, 2015 at 3:01 PM, Ravi Varadhan ravi.varad...@jhu.edu wrote:

 Thank you all.  I did think about declaring `pi' as a special constant,
 but for some reason didn't actually try it.
 Would it be easy to have the mpfr() written such that its argument is
 automatically of extended precision? In other words, if I just called:
 mpfr(exp(sqrt(163)*pi, 120), then all the constants, 163, pi, are
 automatically of 120 bits precision.

 Is this easy to do?

 Best,
 Ravi
  
 From: David Winsemius dwinsem...@comcast.net
 Sent: Friday, July 3, 2015 2:06 PM
 To: John Nash
 Cc: r-help; Ravi Varadhan
 Subject: Re: [R] : Ramanujan and the accuracy of floating point
 computations - using Rmpfr in R

  On Jul 3, 2015, at 8:08 AM, John Nash john.n...@uottawa.ca wrote:
 
 
 
 
  Third try -- I unsubscribed and re-subscribed. Sorry to Ravi for extra
  traffic.
 
  In case anyone gets a duplicate, R-help bounced my msg from non-member
 (I changed server, but not email yesterday,
  so possibly something changed enough). Trying again.
 
  JN
 
  I got the Wolfram answer as follows:
 
  library(Rmpfr)
  n163 - mpfr(163, 500)
  n163
  pi500 - mpfr(pi, 500)
  pi500
  pitan - mpfr(4, 500)*atan(mpfr(1,500))
  pitan
  pitan-pi500
  r500 - exp(sqrt(n163)*pitan)
  r500
  check - 262537412640768743.25007259719818568887935385...
  savehistory(jnramanujan.R)
 
  Note that I used 4*atan(1) to get pi.

 RK got it right by following the example in the help page for mpfr:

 Const(pi, 120)

 The R `pi` constant is not recognized by mpfr as being anything other than
 another double .


 There are four special values that mpfr recognizes.

 —
 Best;
 David


  It seems that may be important,
  rather than converting.
 
  JN
 
  On 15-07-03 06:00 AM, r-help-requ...@r-project.org wrote:
 
  Message: 40
  Date: Thu, 2 Jul 2015 22:38:45 +
  From: Ravi Varadhan ravi.varad...@jhu.edu
  To: 'Richard M. Heiberger' r...@temple.edu, Aditya Singh
   aps...@yahoo.com
  Cc: r-help r-help@r-project.org
  Subject: Re: [R] : Ramanujan and the accuracy of floating point
   computations - using Rmpfr in R
  Message-ID:
   14ad39aaf6a542849bbf3f62a0c2f...@dom-eb1-2013.win.ad.jhu.edu
  Content-Type: text/plain; charset=utf-8
 
  Hi Rich,
 
  The Wolfram answer is correct.
  http://mathworld.wolfram.com/RamanujanConstant.html
 
  There is no code for Wolfram alpha.  You just go to their web engine
 and plug in the expression and it will give you the answer.
  http://www.wolframalpha.com/
 
  I am not sure that the precedence matters in Rmpfr.  Even if it does,
 the answer you get is still wrong as you showed.
 
  Thanks,
  Ravi
 
  -Original Message-
  From: Richard M. Heiberger [mailto:r...@temple.edu]
  Sent: Thursday, July 02, 2015 6:30 PM
  To: Aditya Singh
  Cc: Ravi Varadhan; r-help
  Subject: Re: [R] : Ramanujan and the accuracy of floating point
 computations - using Rmpfr in R
 
  There is a precedence error in your R attempt.  You need to convert
  163 to 120 bits first, before taking
  its square root.
 
  exp(sqrt(mpfr(163, 120)) * mpfr(pi, 120))
  1 'mpfr' number of precision  120   bits
  [1] 262537412640768333.51635812597335712954
 
  ## just the last four characters to the left of the decimal point.
  tmp - c(baseR=8256, Wolfram=8744, Rmpfr=8333, wrongRmpfr=7837)
  tmp-tmp[2]
  baseRWolfram  Rmpfr wrongRmpfr
   -488  0   -411   -907
  You didn't give the Wolfram alpha code you used.  There is no way of
 verifying the correct value from your email.
  Please check that you didn't have a similar precedence error in that
 code.
 
  Rich
 
 
  On Thu, Jul 2, 2015 at 2:02 PM, Aditya Singh via R-help 
 r-help@r-project.org wrote:
  Ravi
 
  I am a chemical engineer by training. Is there not something like law
 of corresponding states in numerical analysis?
 
  Aditya
 
 
 
  --
  On Thu 2 Jul, 2015 7:28 AM PDT Ravi Varadhan wrote:
 
  Hi,
 
  Ramanujan supposedly discovered that the number, 163, has this
 interesting property that exp(sqrt(163)*pi), which is obviously a
 transcendental number, is real close to an integer (close to 10^(-12)).
 
  If I compute this using the Wolfram alpha engine, I get:
  262537412640768743.25007259719818568887935385...
 
  When I do this in R 3.1.1 (64-bit windows), I get:
  262537412640768256.
 
  The absolute error between the exact and R's value is 488, with a
 relative error of about 1.9x10^(-15).
 
  In order to replicate Wolfram Alpha, I tried doing this in Rmfpr
 but I am unable to get accurate results:
 
  

Re: [R] Spline Graphs

2015-06-29 Thread Gabor Grothendieck
Sorry if this appears twice but I am not sure the first attempt got through.

This is not much to go on but here is a short self contained example which
creates a longitudinal data frame L in long form from the built in data
frame BOD and then plots it as points and splines using lattice:

L - rbind(cbind(Id = 1, BOD), cbind(Id = 2, BOD))
library(lattice)
xyplot(demand ~ Time | Id, L, type = c(p, spline))

On Mon, Jun 29, 2015 at 2:28 AM, deva d devazresea...@gmail.com wrote:

 I wish to analyse longitudinal data and fit spline graphs to it looking to
 the data pattern.

 can someone suggest some starting point, and package in R to be used for
 it.

 what would be the requirement for structuring the raw data.

 **

 *Deva*

 [[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.




-- 
Statistics  Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

[[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.


Re: [R] help for lay person assisting R user with disability

2015-06-18 Thread Gabor Grothendieck
On Thu, Jun 18, 2015 at 10:32 AM, Courtney Bryant cbry...@andrew.cmu.edu
wrote:

 Good Morning,
 I am currently working with a disabled R user who is a student here at
 CMU.  The student has both sight and mobility issues.  The student has
 asked for an assistant who is well versed in R to enter data for her, which
 we are having a hard time finding.  I would like information from R
 developers/users about how/how well R interfaces with Excel (an easier
 skill set to find!)   In your opinion, could it be as easy as uploading
 data from excel into R?

 Also, do you know of a way to enlarge the R interface or otherwise assist
 in making the program accessible to a low vision person?  My  limited
 understanding leads me to believe that screen magnifiers like zoom text
 don't work particularly well.  If you have information on that, I would
 very much appreciate it.

 Thanks for your help and for bearing with me!
 Courtney


1. If the data file is in the form of a rectangular table with rows and
columns and the first row is a header row then if, in Excel, it is saved as
a .csv file it can be read into R like this:

   DF - read.csv(/Users/JoeDoe/myspreadsheet.csv)

2. The openxlsx, readxl (and a number of other packages) can alternetely be
used to directly read in an xls or xlsx file, e.g.

  install.packages(readxl)
  library(readxl)
  DF - read_excel(/Users/JoeDoe/myspreadsheet.xlsx)

3. The Windows magnifier that comes with Windows does work with R.

-- 
Statistics  Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

[[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.


Re: [R] Toronto CRAN mirror 403 error?

2015-05-29 Thread Gabor Grothendieck
On Fri, May 29, 2015 at 10:12 PM, Mark Drummond m...@markdrummond.ca wrote:
 I've been getting a 403 when I try pulling from the Toronto CRAN mirror
 today.

 http://cran.utstat.utoronto.ca/

 Is there a contact list for mirror managers?


See the cran_mirrors.csv file in
  R.home(doc)
of your R distribution.

__
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] Extract month and year as one column

2015-03-13 Thread Gabor Grothendieck
On Fri, Mar 13, 2015 at 11:36 AM, Kumsaa waddee...@gmail.com wrote:
 How could I extract both month and year from a date? I know how to
 separately extract both using lubridate package:

 df$month - month(df$date)
 df$year- year(df$date)

 I wish to extract year and month as one column

 dput(mydf)
 structure(list(date = structure(c(14975, 14976, 14977, 14978,
 14979, 14980, 14981, 14982, 14983, 14984, 14985, 14986, 14987,
 14988, 14989, 15340, 15341, 15342, 15343, 15344, 15345, 15346,
 15347, 15348, 15349, 15350, 15351, 15352, 15353, 15354), class = Date),
 temp = c(6.5140544091004, 3.69073712745884, 3.04839429519466,
 9.16988228171461, -1.17176248610603, 2.88216040747883,
 4.98853844809017,
 4.07520306701834, 9.82902813943658, 2.79305715971987, 8.04721677924611,
 7.50667729759095, 2.91055000121842, 1.65559895014064, 4.8019596483372,
 16.2567986804179, 13.3352908067145, 16.6955807821108, 6.28373374879922,
 6.97181051627531, 5.74282686202818, 4.37018386569785, 12.5725962512824,
 4.6583055309578, 8.76457542037641, 10.7070862034423, 12.84023567151,
 5.78620621848167, 5.98643374478599, 13.0993210289842)), .Names =
 c(date,
 temp), row.names = c(NA, 30L), class = data.frame)


The zoo package has a yearmon class that represents dates as year
and month with no day:

 library(zoo)
 transform(mydf, yearmon = as.yearmon(date))
date  temp  yearmon
1 2011-01-01  6.514054 Jan 2011
2 2011-01-02  3.690737 Jan 2011
3 2011-01-03  3.048394 Jan 2011
4 2011-01-04  9.169882 Jan 2011
5 2011-01-05 -1.171762 Jan 2011
6 2011-01-06  2.882160 Jan 2011
etc.

__
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] Substring replacement in string

2015-02-28 Thread Gabor Grothendieck
On Fri, Feb 27, 2015 at 5:19 PM, Alrik Thiem alrik.th...@gmail.com wrote:
 I would like to replace all lower-case letters in a string that are not part
 of certain fixed expressions. For example, I have the string:

 pmin(pmax(pmin(x1, X2), pmin(X3, X4)) == Y, pmax(Z1, z1))

 Where I would like to replace all lower-case letters that do not belong to
 the functions pmin and pmax by 1 - toupper(...) to get

 pmin(pmax(pmin(1 - X1, X2), pmin(X3, X4)) == Y, pmax(Z1, 1 - Z1))


Assuming x is the input string:

gsub((\\b[a-oq-z][a-z0-9]+), 1-\\U\\1, x, perl = TRUE)
## [1] pmin(pmax(pmin(1-X1, X2), pmin(X3, X4)) == Y, pmax(Z1, 1-Z1))



-- 
Statistics  Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] Substring replacement in string

2015-02-28 Thread Gabor Grothendieck
Replace the + (i.e. 1 or more) in the pattern with a * (i.e. 0 or more):

   x - pmin(pmax(pmin(a,B),pmin(a,C,d))==Y,pmax(E,e))

   gsub((\\b[a-oq-z][a-z0-9]*), 1-\\U\\1, x, perl = TRUE)

giving:

   [1] pmin(pmax(pmin(1-A,B),pmin(1-A,C,1-D))==Y,pmax(E,1-E))

Here is a visualization of the regular expression:

   https://www.debuggex.com/i/5ByOCQS2zIdPEf-f.png


On Sat, Feb 28, 2015 at 8:16 AM, Alrik Thiem alrik.th...@gmail.com wrote:
 Dear Gabor,

 Many thanks. Works like a charm, but I can't get it to work with

 pmin(pmax(pmin(a,B),pmin(a,C,d))==Y,pmax(E,e))

 i.e., with strings where there're no integers following the components in the 
 pmin/pmax functions. Could this be generalized to handle both cases?

 Best wishes,
 Alrik

 -Ursprüngliche Nachricht-
 Von: Gabor Grothendieck [mailto:ggrothendi...@gmail.com]
 Gesendet: Samstag, 28. Februar 2015 13:35
 An: Alrik Thiem
 Cc: r-help@r-project.org
 Betreff: Re: [R] Substring replacement in string

 On Fri, Feb 27, 2015 at 5:19 PM, Alrik Thiem alrik.th...@gmail.com wrote:
 I would like to replace all lower-case letters in a string that are not part
 of certain fixed expressions. For example, I have the string:

 pmin(pmax(pmin(x1, X2), pmin(X3, X4)) == Y, pmax(Z1, z1))

 Where I would like to replace all lower-case letters that do not belong to
 the functions pmin and pmax by 1 - toupper(...) to get

 pmin(pmax(pmin(1 - X1, X2), pmin(X3, X4)) == Y, pmax(Z1, 1 - Z1))


 Assuming x is the input string:

 gsub((\\b[a-oq-z][a-z0-9]+), 1-\\U\\1, x, perl = TRUE)
 ## [1] pmin(pmax(pmin(1-X1, X2), pmin(X3, X4)) == Y, pmax(Z1, 1-Z1))



 --
 Statistics  Software Consulting
 GKX Group, GKX Associates Inc.
 tel: 1-877-GKX-GROUP
 email: ggrothendieck at gmail.com




-- 
Statistics  Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com

__
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] sqldf() difference between R 3.1.2 and 3.0.1

2015-02-11 Thread Gabor Grothendieck
On Wed, Feb 11, 2015 at 9:45 AM, Doran, Harold hdo...@air.org wrote:
 I have a function written and tested using R 3.0.1 and sqldf_0.4-7.1 that 
 works perfectly. However, using this same code with R 3.1.2 and sqldf_0.4-10 
 yields the error below that I am having a difficult time deciphering. Hence, 
 same code behaves differently on different versions of R and sqldf().

 Error in sqliteSendQuery(con, statement, bind.data) :
   error in statement: no such column: V1


 Reproducible example below as well as complete sessionInfo all provided below.


 My function and code using the function are below.

 dorReader - function(dorFile, layout, sepChar = '\n'){
 sepChar - as.character(sepChar)
 dorFile - as.character(dorFile)
 layout$type2 - ifelse(layout$type == 'C', 'character',
 
 ifelse(layout$type == 'N', 'numeric', 'Date'))
 dor - file(dorFile)
 attr(dor, file.format) - list(sep = sepChar)
 getVars - paste(select,
paste(substr(V1, , layout$Start, , ,
  layout$Length, ) ', layout$Variable.Name, ', 
 collapse = , ), from dor)
 dat - sqldf(getVars)

 classConverter - function(obj, types){
 out - lapply(1:length(obj),FUN = 
 function(i){FUN1 - switch(types[i],character = as.character,numeric = 
 as.numeric,factor = as.factor, Date = as.character); FUN1(obj[,i])})
 names(out) - colnames(obj)
 as.data.frame(out)
 }
 dat - classConverter(dat, layout$type2)
 names(dat) - layout$Variable.Name
 dat
 }

 ### contents of fwf file 'sample.txt'
 1234567
 1234567
 1234567
 1234567
 1234567
 1234567
 1234567
 1234567

 layout - data.frame(Variable.Name =c('test1', 'test2'), Length = c(3,4), 
 Start =c(1,4), End = c(3,7), type = c('N', 'N'))

 tmp - dorReader('sample.txt', layout)

sqldf is documented to use the sqliteImportFile defaults for
file.format components.  It may be that RSQLite 1.0 has changed the
default for header in sqliteImportFile.

Try replacing your statement that sets file.format with this:

attr(dor, file.format) - list(sep = sepChar, header = FALSE)

__
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.


  1   2   3   4   5   6   7   8   9   10   >