Re: [R] Optimization failed in fitdistr (Weibull distribution)

[peter dalgaard]

On 28 Oct 2013, at 13:07 , kmmoon100  wrote:

> Hello everyone,
> 
> This is Kangmin.
> 
> I am trying to produce shape and scale of my wind data. My data is based on
> wind speed frequency with 1km/hr increment. data is described below.
> 
> Windspeed (km/h)Frequency
> 1 351
> 2 147
> 3 317
> 4 378
> 5 527
> 6 667
> 7 865
> 8 970
> 9 987
> 10907
> 11905
> 12642
> 131000
> 14983
> 15847
> 16842
> 17757
> 18698
> 19632
> 20626
> 21599
> 22529
> 23325
> 24391
> 25356
> 26267
> 27230
> 28223
> 29190
> 30142
> 31124
> 32104
> 3397
> 3437
> 3562
> 3646
> 3742
> 3824
> 399
> 4013
> 419
> 425
> 432
> 
> R codes to calculate shape and scale are described below:
> 
> Pine.windfrequency.4weeks<-read.table("C:/Users/kmoon/Documents/Pine_frequency_4weeks.csv",header=TRUE,sep=",")
> fitdistr(Pine.windfrequency.4weeks$Frequency, densfun="weibull")
> 
> I have got an error message when I was using 'fitdistr' function
> 
> "Error in fitdistr(Pine.windfrequency.4weeks$Frequency, densfun = "weibull")
> : 
>  optimization failed"
> 
> Please help me calculating shape and scale of weibull distribution.
> 
> And please understand that I am not an user familiar with R program but I am
> really trying to make my analysis work on R!

There really is no substitute for knowledge and understanding! Did it not occur 
to you that the Windspeed column needs to enter into your analysis? 

I suppose you wanted the following:

> tt<-read.delim("/tmp/foo")
> summary(tt)
 Windspeed..km.h.   Frequency 
 Min.   : 1.0 Min.   :   2.0  
 1st Qu.:11.5 1st Qu.: 100.5  
 Median :22.0 Median : 351.0  
 Mean   :22.0 Mean   : 415.7  
 3rd Qu.:32.5 3rd Qu.: 682.5  
 Max.   :43.0 Max.   :1000.0  
> x <- rep(tt$Windspeed..km.h., tt$Frequency)
> library(MASS) 
> fitdistr(x, densfun="weibull")
  shape scale   
   1.99900495   16.43640142 
 ( 0.01174133) ( 0.06468371)
Warning messages:
1: In densfun(x, parm[1], parm[2], ...) : NaNs produced
2: In densfun(x, parm[1], parm[2], ...) : NaNs produced
3: In densfun(x, parm[1], parm[2], ...) : NaNs produced
4: In densfun(x, parm[1], parm[2], ...)


> 
> Thank you!!!
> 
> Kangmin.
> 
> 
> 
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Optimization-failed-in-fitdistr-Weibull-distribution-tp4679167.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

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

Re: [R] Optimization failed in fitdistr (Weibull distribution)

[Berend Hasselman]

On 29-10-2013, at 00:35, kmmoon100  wrote:

> Hi Berend,
> 
> Thank you for your reply.
> How can I use dput function for this type of data?
> I looked up the description of the function but I still can't understand how
> to use it for solving my error.
> 

You don't use dput() to solve your error.
You use dput to get a plain text representation of your data.

Select the output in any editor, copy and paste it in a mail to the list.
If the output is huge then of course you can do dput(head(….,…)) to create a 
reproducible example.
And then maybe people on the list can suggest solutions.

Berend

> 
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Optimization-failed-in-fitdistr-Weibull-distribution-tp4679178p4679232.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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

Re: [R] Optimization failed in fitdistr (Weibull distribution)

[kmmoon100]

Hi Berend,

Thank you for your reply.
How can I use dput function for this type of data?
I looked up the description of the function but I still can't understand how
to use it for solving my error.

Regards,

Kangmin.



--
View this message in context: 
http://r.789695.n4.nabble.com/Optimization-failed-in-fitdistr-Weibull-distribution-tp4679178p4679232.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] need some help with this example

[Glenn Schultz]

I have a class Details that contains information needed by FirstSet to do some 
calculations then a super class that returns Details and FirstSet.  The problem 
seems to be in FirstSet where I use the function getAnumber(id).  

setClass("Details",
         representation(
           ID = "character",
           Anumber = "numeric"))

setGeneric("Details",
           def = function(object){standardGeneric("Details")})

setMethod("initialize",
          signature(.Object = "Details"),
          function(.Object, ID = "character", Anumber = numeric()){
            .Object@ID = ID
            .Object@Anumber = 2
            return(.Object)
          })

# getter for A number
setGeneric("getAnumber", function(object){standardGeneric("getAnumber")})

setMethod("getAnumber","Details",
          function(object){return(Object@Anumber)})

setClass("FirstSet",
         representation(
           Anothernumber = "numeric"))

setGeneric(
  name = "FirstSet",
  def = function(object){standardGeneric("FirstSet")}
)
setMethod("initialize",
          signature(.Object = "FirstSet"),
          function (.Object, id = "character", multiplier = numeric())
          { x = getAnumber(id)
            y = x * multiplier
            .Object@Anothernumber = y
            return(.Object)
          }
)

setClass("Super", contains = c("Details", "FirstSet"))

setGeneric("Super",
           def = function(object){standardGeneric("Super")})

setMethod("initialize",
          signature(.Object = "Super"),
          function(.Object, id = "character", Anumber = Anumber()){
            Details<- new("Details", ID = id, Anumber = number)
            FirstSet <- new("FirstSet", Anothernumber = Anothernumber)
            Super <- new("Super", Details, FirstSet)
            return(.Object)
          })

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

Re: [R] Optimize function in R: unable to find maximum of logistic function

[Shantanu MULLICK]

Hello,

Thanks a lot for the replies.

I tried to use the "optim" function in R, and chose the "L-BFGS-B" method
of optimization. Now the result is very sensitive to the starting values. I
use the same function, but I give the line of  code I use for "optim".

optim(0.25, f, method= "L-BFGS-B", lower = 0, upper=5). Here 0.25 is the
starting value.

Result: $ par = 0, $Value = 1 (this is correct). $par indicates that the
maximum value of the function is at 0, and $ Value says that the maximum
value is 1.

However, I produced a quick table for the results with different starting
values:
Starting Value $par $Value
 Result
0.25 0  1
  Correct
0.5  0.4999  4.1223e-09
Wrong
11   0
   Wrong
3 3   0
  Wrong

I wanted to optimise the logistic function, as the logistic function
approximates an indicator function. The indicator function is defined as
follows:  I(u), where I(u) = 1, where u = 0, and I(u) = 0, for all other
 u. u is strictly positive.

My objective function contains a indicator function.As the indicator
function is non-continuous, the optimization would become non-convex. To
bypass this limitation, I used the logistic function to approximate the
indicator function.

A simple example: f(u) = I(u) + u/6;Maximise f(u) where u is between (0
to 5).
Correct Answer: f(u) =1, where u=0.
Code:
f= function (k) {
T_s = 20
result = (2- 2/(1+ exp(-2*T_s*k))) + (k/6)
return(result)
}
optim(0.15, f, method= "L-BFGS-B", lower = 0, upper=5,control =
list(fnscale=-1) )
## Here 0.15 is the starting value. I produce a similar table as above:

Starting Value $par $Value
 Result
0.15 0  1
  Correct
0.16 5   0.833
 Wrong
0.17 5   0.833
 Wrong
2   5   0.833
   Wrong

Questions:
1. If my objective function contains an Indicator function, is there an
efficient way to implement it?
2. Does R have any quick non-convex optimiser functions?

Best Regards,
Shantanu Mullick
PhD Candidate
ESSEC Business School

On 28 October 2013 22:01, Rolf Turner  wrote:

>
> This could be described as a bug, perhaps.  Or it could be described as
> an indication that numerical optimization is inevitably tricky. Notice that
> if you narrow down your search interval from [0,5] to [0,0.5] you get the
> right answer:
>
> > optimize(f, c(0, 0.5), maximum= TRUE,tol=1e-10)
> $maximum
> [1] 4.192436e-11
>
> $objective
> [1] 1
>
> I guess there's a problem with finding a gradient that is effectively
> (numerically) zero when "k" is equal to 5.
>
>
> cheers,
>
> Rolf Turner
>
>
>
> On 10/29/13 06:00, Shantanu MULLICK wrote:
>
>> Hello Everyone,
>>
>> I want to perform a 1-D optimization by using the optimize() function. I
>> want to find the maximum value of a "logistic" function. The optimize()
>> function gives the wrong result.
>>
>> My code:
>> f= function (k) {
>> T_s = 20
>> result = (2- 2/(1+ exp(-2*T_s*k)))
>> return(result)
>> }
>> optimize(f, c(0, 5), tol = 0.1, maximum= TRUE)
>>
>> The maximum value for the function happens at k=0, and the maximum value
>> is
>> 1. Yet  the optimise function, says that the maximum value happens at k=
>> 4.9995, and the maximum value is 0.
>>
>>

[[alternative HTML version deleted]]

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

Re: [R] Optimize function in R: unable to find maximum of logistic function

[William Dunlap]

> One can do better yet by using the logspace_add(logx, logy) function
> available in the C API to R, but there is currently not an R-language 
> interface
> to that (or logspace_subtract(logx,logy)).  They calculate log(exp(logx) +- 
> exp(logy))
> with minimal roundoff error.
> 
> You could probably also use the built-in plogis() function, with its log=TRUE 
> and
> perhaps lower.tail=FALSE arguments.

Here are functions, "f2" and "f3", which compute log(f()), via the above 
algorithms,
along with your original "f" and my first suggestion "f1".  f2 and f3 give 
almost
identical results over a wide range of inputs (pos. and neg.).  f1 agrees with 
them
for positive inputs but falls apart for larger negative inputs.  Any of f1, f2, 
or f3 will
do your optimization, although I don't know why you want to optimize a monotonic
function.

f <-
function (k, T_s = 20) {
2 - 2/(1 + exp(-2*T_s*k))
}
f1 <-
function(k, T_s = 20) {
# log(f(k))
log(2) - 2*T_s*k - log1p(exp(-2*T_s*k))
}
f2 <-
function(k, T_s = 20) {
# log(f(k))
log(2) - logspace_add(2*T_s*k, 0)
}
logspace_add <-
function(logx, logy)
{
# log(exp(logx) + exp(logy))
mn <- pmin(logx, logy)
mx <- pmax(logx, logy)
mx + log1p(exp(mn - mx))
}
f3 <-
function(k, T_s = 20) {
# log(f(k))
log(2) + plogis(k, lower.tail=FALSE, log=TRUE, scale=1/(T_s*2))
}

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
> Behalf
> Of William Dunlap
> Sent: Monday, October 28, 2013 3:03 PM
> To: Rolf Turner; Shantanu MULLICK
> Cc: r-help@r-project.org
> Subject: Re: [R] Optimize function in R: unable to find maximum of logistic 
> function
> 
> Note that f(x) returns exactly zero for all x above 1, hence its estimated
> derivative is 0 everywhere in that region.   You are asking optimize to find
> the non-zero part of a function which is 0 in more than 80% of its domain.
> 
> You can put a trace on f() to see where optimize() looks:
>> trace(f, quote(cat("k=", deparse(k), "\n")), print=FALSE)
>   [1] "f"
>   > optimize(f, c(0, 5), tol = 1e-14, maximum= TRUE)
>   k= 1.90983005625053
>   k= 3.09016994374947
>   k= 3.81966011250105
>   k= 4.27050983124842
>   k= 4.54915028125263
>   k= 4.72135954999579
>   k= 4.82779073125683
>   k= 4.89356881873896
>   ... eventually homing in on 5 ...
> It looks in a few places, but f() is 0 in all of them so it figures that
> is the mininum and the maximum value it ever takes.  I suppose
> you could ask that optimize look in more places for a place with
> a non-zero derivative but it would be very expensive to look in all
> (c. 2^50) places.
> 
> You probably should have it maximize the log of your objective
> function, which gives you more range to play with, and you will
> have to do some numerical analysis to minimize underflow and
> roundoff error.  E.g., one could use
> f1 <- function (k)  {
> # log(f(k))
> T_s <- 20
> log(2) - 2 * T_s * k - log1p(exp(-2 * T_s * k))
>}
> where log1p(x) calculates log(1+x) more accurately than a computer's
> 'log' and '+' can.  (Remember that 1+10^-17 exactly equals 1 here.)
> 
> One can do better yet by using the logspace_add(logx, logy) function
> available in the C API to R, but there is currently not an R-language 
> interface
> to that (or logspace_subtract(logx,logy)).  They calculate log(exp(logx) +- 
> exp(logy))
> with minimal roundoff error.
> 
> You could probably also use the built-in plogis() function, with its log=TRUE 
> and
> perhaps lower.tail=FALSE arguments.
> 
> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com
> 
> 
> > -Original Message-
> > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
> > Behalf
> > Of Rolf Turner
> > Sent: Monday, October 28, 2013 2:01 PM
> > To: Shantanu MULLICK
> > Cc: r-help@r-project.org
> > Subject: Re: [R] Optimize function in R: unable to find maximum of logistic 
> > function
> >
> >
> > This could be described as a bug, perhaps.  Or it could be described as
> > an indication that numerical optimization is inevitably tricky. Notice that
> > if you narrow down your search interval from [0,5] to [0,0.5] you get the
> > right answer:
> >
> >  > optimize(f, c(0, 0.5), maximum= TRUE,tol=1e-10)
> > $maximum
> > [1] 4.192436e-11
> >
> > $objective
> > [1] 1
> >
> > I guess there's a problem with finding a gradient that is effectively
> > (numerically) zero when "k" is equal to 5.
> >
> >
> >  cheers,
> >
> >  Rolf Turner
> >
> >
> > On 10/29/13 06:00, Shantanu MULLICK wrote:
> > > Hello Everyone,
> > >
> > > I want to perform a 1-D optimization by using the optimize() function. I
> > > want to find the maximum value of a "logistic" function. The optimize()
> > > function gives the wrong result.
> > >
> > > My code:
> > > f= function (k) {
> > > T_s = 20
> > > result = (2- 2/(1+ exp(-2*T_s*k)))
> > > return(result)
> > > }
> 

Re: [R] pie graphs in log scale axis

[jwd]

On Mon, 28 Oct 2013 15:03:36 -0700
jwd  wrote:

> On Mon, 28 Oct 2013 13:40:51 +
> "Ortiz, John"  wrote:
...
> > 
> It is not clear what you want.  Why is there an 'x = "log"' term in
> plot()?  If you want to plot log(x) try using:

Got that backward typing, should be log = "x"

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

Re: [R] pie graphs in log scale axis

[jwd]

On Mon, 28 Oct 2013 13:40:51 +
"Ortiz, John"  wrote:

> 
> Dear list users,
> 
> 
> I'm doing a plot integrating Grid output with Base Graphics output
> (gridBase, Murrell 2012).
> 
> My goal is to produce a xy plot where each point is represented by a
> pie.  I could get it using the attach code, but now I want to change
> the X axis to log scale.
> 
> when I introduce the log="x" parameter in the line:
> 
> plot(x, y, xlim=c(0.1, 1.2), log="x", ylim=c(0.1, 1.2), type="n")
> 
> I get this warning!
> 
> vps <- baseViewports()
> Warning message:
> In baseViewports() : viewport scales NOT set to user coordinates
> 
It is not clear what you want.  Why is there an 'x = "log"' term in
plot()?  If you want to plot log(x) try using:

plot(log(x), y, xlim=c(0.1, 1.2), ylim=c(0.1, 1.2), type="n")

That generates an x/y plot with pies at various points.

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

Re: [R] Optimize function in R: unable to find maximum of logistic function

[William Dunlap]

Note that f(x) returns exactly zero for all x above 1, hence its estimated
derivative is 0 everywhere in that region.   You are asking optimize to find
the non-zero part of a function which is 0 in more than 80% of its domain.

You can put a trace on f() to see where optimize() looks:
   > trace(f, quote(cat("k=", deparse(k), "\n")), print=FALSE)
  [1] "f"
  > optimize(f, c(0, 5), tol = 1e-14, maximum= TRUE)
  k= 1.90983005625053 
  k= 3.09016994374947 
  k= 3.81966011250105 
  k= 4.27050983124842 
  k= 4.54915028125263 
  k= 4.72135954999579 
  k= 4.82779073125683 
  k= 4.89356881873896
  ... eventually homing in on 5 ...
It looks in a few places, but f() is 0 in all of them so it figures that
is the mininum and the maximum value it ever takes.  I suppose
you could ask that optimize look in more places for a place with
a non-zero derivative but it would be very expensive to look in all
(c. 2^50) places.

You probably should have it maximize the log of your objective
function, which gives you more range to play with, and you will
have to do some numerical analysis to minimize underflow and
roundoff error.  E.g., one could use
f1 <- function (k)  {
# log(f(k))
T_s <- 20
log(2) - 2 * T_s * k - log1p(exp(-2 * T_s * k))
   }
where log1p(x) calculates log(1+x) more accurately than a computer's
'log' and '+' can.  (Remember that 1+10^-17 exactly equals 1 here.)

One can do better yet by using the logspace_add(logx, logy) function
available in the C API to R, but there is currently not an R-language interface
to that (or logspace_subtract(logx,logy)).  They calculate log(exp(logx) +- 
exp(logy))
with minimal roundoff error.

You could probably also use the built-in plogis() function, with its log=TRUE 
and
perhaps lower.tail=FALSE arguments.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
> Behalf
> Of Rolf Turner
> Sent: Monday, October 28, 2013 2:01 PM
> To: Shantanu MULLICK
> Cc: r-help@r-project.org
> Subject: Re: [R] Optimize function in R: unable to find maximum of logistic 
> function
> 
> 
> This could be described as a bug, perhaps.  Or it could be described as
> an indication that numerical optimization is inevitably tricky. Notice that
> if you narrow down your search interval from [0,5] to [0,0.5] you get the
> right answer:
> 
>  > optimize(f, c(0, 0.5), maximum= TRUE,tol=1e-10)
> $maximum
> [1] 4.192436e-11
> 
> $objective
> [1] 1
> 
> I guess there's a problem with finding a gradient that is effectively
> (numerically) zero when "k" is equal to 5.
> 
> 
>  cheers,
> 
>  Rolf Turner
> 
> 
> On 10/29/13 06:00, Shantanu MULLICK wrote:
> > Hello Everyone,
> >
> > I want to perform a 1-D optimization by using the optimize() function. I
> > want to find the maximum value of a "logistic" function. The optimize()
> > function gives the wrong result.
> >
> > My code:
> > f= function (k) {
> > T_s = 20
> > result = (2- 2/(1+ exp(-2*T_s*k)))
> > return(result)
> > }
> > optimize(f, c(0, 5), tol = 0.1, maximum= TRUE)
> >
> > The maximum value for the function happens at k=0, and the maximum value is
> > 1. Yet  the optimise function, says that the maximum value happens at k=
> > 4.9995, and the maximum value is 0.
> >
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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

Re: [R] Revo R for Arima Implementation

[Rolf Turner]

Your question still makes no sense at all to me.  You provide no example 
code

and no data.  You do not specify, in any comprehensible way, what model you
are trying to fit.

And as Jeff Newmiller said, if you are concerned with "Revo R" why are 
you posting

to "r-help"?

Please get your act together or stop cluttering up the r-help list.

cheers,

Rolf Turner

On 10/29/13 02:26, Anindita Chattopadhyay wrote:

Hi Rolf,

Thanks for the response. I have re-phrased the problem. Hope this will help.

We're working on a project where our model tries to predict the value of 
bookings as a time series of past bookings as well as some external flags.
The data was found to be stationary using the Dicky-Fuller test.
We have only a non-seasonal component for the AR and MA terms (small p & q 
only).
No differencing was done.

Then, we used values of p = 1, 7, 8 and q = 1, 7.
The main issue is that the p values are additive and q values are 
multiplicative.
We need to understand how we can implement this in Revo R.

Just to add, since we have external flags , we cannot make use of SARIMA 
function.
Please let me know if any other questions.

Thank you in advance!

Regards,
Anindita Chattopadhyay | +919886800606 | www.mu-sigma.com |

-Original Message-
From: Rolf Turner [mailto:r.tur...@auckland.ac.nz]
Sent: Saturday, October 26, 2013 1:56 AM
To: Anindita Chattopadhyay
Cc: r-help@r-project.org; Harish K
Subject: Re: [R] Revo R for Arima Implementation


Your question is pretty well totally opaque to me.  Describe the model you want 
to fit in mathematical terms, rather than referring to The Package That Must 
Not Be Named.  This is the ***R***  list.

It is possible that you might want to make use of the "seasonal"
argument to the
arima() function.

  cheers,

  Rolf Turner

On 10/26/13 01:19, Anindita Chattopadhyay wrote:

Hello There,

We have used ARIMA(multiplicative MA and additive AR) model in SAS to come to 
our results. Please let me know how could that be implemented in R.

Like,  In SAS we can pass the variable  as AR= (1,7) and MA=(1)(7) which is 
additive and multiplicative respectively.

In R we have the option as : order = c(0, 0, 0) which is (p,d,q).Is there any 
specific way/code such that we specify the additive  AR and multiplicative MA 
simultaneously?


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

Re: [R] Optimize function in R: unable to find maximum of logistic function

[Rolf Turner]

This could be described as a bug, perhaps.  Or it could be described as
an indication that numerical optimization is inevitably tricky. Notice that
if you narrow down your search interval from [0,5] to [0,0.5] you get the
right answer:

> optimize(f, c(0, 0.5), maximum= TRUE,tol=1e-10)
$maximum
[1] 4.192436e-11

$objective
[1] 1

I guess there's a problem with finding a gradient that is effectively
(numerically) zero when "k" is equal to 5.


cheers,

Rolf Turner


On 10/29/13 06:00, Shantanu MULLICK wrote:

Hello Everyone,

I want to perform a 1-D optimization by using the optimize() function. I
want to find the maximum value of a "logistic" function. The optimize()
function gives the wrong result.

My code:
f= function (k) {
T_s = 20
result = (2- 2/(1+ exp(-2*T_s*k)))
return(result)
}
optimize(f, c(0, 5), tol = 0.1, maximum= TRUE)

The maximum value for the function happens at k=0, and the maximum value is
1. Yet  the optimise function, says that the maximum value happens at k=
4.9995, and the maximum value is 0.



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

Re: [R] Optimize function in R: unable to find maximum of logistic function

[Rui Barradas]

Hello,

This seems like a bug, removing 'maximum = TRUE' has no effect except 
that the returned value says it's a minimum, not a maximum.


Rui Barradas

Em 28-10-2013 17:00, Shantanu MULLICK escreveu:

Hello Everyone,

I want to perform a 1-D optimization by using the optimize() function. I
want to find the maximum value of a "logistic" function. The optimize()
function gives the wrong result.

My code:
f= function (k) {
T_s = 20
result = (2- 2/(1+ exp(-2*T_s*k)))
return(result)
}
optimize(f, c(0, 5), tol = 0.1, maximum= TRUE)

The maximum value for the function happens at k=0, and the maximum value is
1. Yet  the optimise function, says that the maximum value happens at k=
4.9995, and the maximum value is 0.

Thanks in advance!

Warm Regards,
Shantanu

[[alternative HTML version deleted]]

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



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

Re: [R] Create Time Lists with a for loop

[Rui Barradas]

Hello,

Instead of string manipulation you should think of POSIXt objects. And 
maybe ?seq.POSIXt will help.



first <- strptime("2011-10-12 10:59:00","%Y-%m-%d %H:%M:%S")
last <- strptime("2011-10-12 11:13:00","%Y-%m-%d %H:%M:%S")

seq(first, last, by = "5 min")


Hope this helps,

Rui Barradas

Em 28-10-2013 18:42, Alaios escreveu:

Dear all,
in my code I have written the following list

TimeFramesShort <-list(c(strptime("2011-10-12 10:59:00","%Y-%m-%d 
%H:%M:%S"),strptime("2011-10-13 11:02:00","%Y-%m-%d %H:%M:%S"))  # rest items were 
truncated


I was wondering if could somehow  take the first element of the list and create 
element that are only 5 minutes apart something like
list(c(strptime("2011-10-12 10:59:00","%Y-%m-%d %H:%M:%S"),strptime("2011-10-12 
11:04:00","%Y-%m-%d %H:%M:%S"),
c(strptime("2011-10-12 11:04:00","%Y-%m-%d %H:%M:%S"),strptime("2011-10-12 
11:09:00","%Y-%m-%d %H:%M:%S"),
c(strptime("2011-10-12 11:09:00","%Y-%m-%d %H:%M:%S"),strptime("2011-10-12 
11:13:00","%Y-%m-%d %H:%M:%S")
...
and so on.

The problem I see is that pure string manipulation will not work as after one 
hour that should also make changes in the hours.

Regards
Alex
[[alternative HTML version deleted]]



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



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

Re: [R] Chess Playing Software Written in R

[Jeff Newmiller]

Please avail yourself of search tools before posting here, and reference why 
they did not answer your question.

RSiteSearch( "chess" )

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

Rob Bernard  wrote:
>Have there been attempts to build a chess engine in R?  If so, is it
>available anywhere?
>
>I recognize that R really isn�t the right language for a chess engine,
>but
>I was more curious if it had been attempted.
>
>Thank you.
>
>Rob Bernard

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] fitdistr: was Heteroscedasity...

[Collin Lynch]

Hello again, first off thank you for your suggestion Mr. Rigby, I'll take
a look at the GAMLSS package.

I have a (slightly) related followup question regarding the 'fitdistr'
function.  I was examining my data, a sample of which is attached, using
this function and I am confused about the interpretation of the loglik
result.

Based upon my experience I understand log-likelyhood value returned from
this function should be in the range -Inf - 0 with values closer to 0
being the best choice.  However when I test this data I get the following
results:

fitdistr(E, "normal")$loglik
[1] 11.15125

fitdistr(E + 1, "lognormal")$loglik
[1] -0.8575117

fitdistr(E, "exponential")$loglik
[1] -73.18107

Is this a sign of error, in my system or on my part?  And if not am I
correct in interpreting lognormal as the best choice?

Thank you in advance,
Collin Lynch.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 line in x, y plot

[Jim Lemon]

On 10/29/2013 02:30 AM, Ahmed Attia wrote:

Hi

I have a question about drawing a linear line in x, y plot. I usually
use the following code, but for this time the x values are to small
(-0.08 to -0.02)

I wrote the following code, but r does not draw the line. However, it
does not give an error when it takes the code.


reg1<- lm(CWSI~NWI, data=Ahmed)

NWI<- seq(-0.08, -0.02, len = -0.02)
lines(NWI, predict(reg1, list(NWI =
NWI)),xlim=c(-0.08,-0.02),ylim=c(0,1),pch=1,col=1,lwd=2, lty=1)

When I wrote the following code

abline(reg1,pch=2,col=2,lwd=2, lty=2,xlim=c(-0.06,-0.02),ylim=c(0.3,0.8))

the line shows up, but I can not control the xlim or ylim. It bascally
goes across the figure


Hi Ahmed,
Have a look at ablineclip in the plotrix package.

Jim

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

Re: [R] indicating significant differences in boxplots

[Jim Lemon]

On 10/29/2013 01:51 AM, Johannes Radinger wrote:

Hi,

I'd like to follow up an older discussion on adding significance stars to
boxplots.
(
http://r.789695.n4.nabble.com/indicating-significant-differences-in-boxplots-td862335.html
)
As suggested by Jim Lemon, there is following solution for a single boxplot:

boxplot(count ~ spray, data = InsectSprays, col = "lightgray")

#and we'll assume that the second and third boxplots require a star:
par(xpd=TRUE)
yrange<-par("usr")[3:4]
ypos<-yrange[2]+diff(yrange)/40
segments(2,ypos,3,ypos)
text(2.5,ypos+diff(yrange)/40,"*",cex=2)
par(xpd=FALSE)

But what if one wants to plot two boxplots (e.g. with each only two groups
which are significantly different) and
add such stars to each of the boxplots. I usually use par(mfrow(1,2)) to
plot all boxplots in a row.
How can I also add segments and stars on top of each plot?


Hi johannes,
The general strategy above can be applied. Fortunately, boxplots are 
centered on successive integers (i.e. 1,2,...) by default. The 
calculations for ypos yield 1/40th of the vertical span of the plot 
above the top of the plot. Just add the segments and stars after each 
plot, as each boxplot advances to the next frame in your plot matrix.


Jim

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

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

[Arnaud Mosnier]

Thanks a lot Simon, that's useful.
I will take a look at this process-pinning things.

Arnaud


2013/10/28 Simon Zehnder 

> First,
>
> use only the number of cores as a number of thread - i.e. I would not use
> hyper threading, etc.. Each core has its own caches and it is always
> fortunate if a process has enough memory; hyper threads use all the same
> cache on the core they are running on. detectCores() gives me for example 4
> - I know I have 2. I would therefore call makeCluster() with nnode = 2.
> mcaffinity() lets you perform a technique called process-pinning (see
> process affinity) and is only possible if the OS supports it. It makes
> sometimes sense to assign certain processes to certain CPUs such that each
> process has enough memory in caches (e.g. for a 16 Core machine using 8
> processes on CPUs 1, 3, 5, 7, 9, 11, 13 and 15; so each process has the
> cache of two CPUs).
> A lot of functions though are not available for Windows.
>
> At first it comes always the problem you want to solve and then you look
> how much memory will be used in a process and how much you have (more often
> the memory bandwidth is the bottleneck and not the computing power). Look
> at the architecture of your chips (how much L1 Cache, how much L2 cache).
> Then you decide how many cores to use and if it makes sense to pin
> processes to certain cores.
>
> There are no general recipes for parallel computing - each problem is
> different. Some problems are even not scalable.
>
> Simon
>
>
> On 28 Oct 2013, at 17:51, Arnaud Mosnier  wrote:
>
> > Thanks Simon,
> >
> > I already read the parallel vignette but I did not found what I wanted.
> > May be you can be more specific on a part of the document that can
> provide me hints !
> >
> > Arnaud
> >
> >
> > 2013/10/28 Simon Zehnder 
> > See library(help = "parallel”)
> >
> >
> > On 28 Oct 2013, at 17:19, Arnaud Mosnier  wrote:
> >
> > > Hi all,
> > >
> > > I am quite new in the world of parallelization and I wonder if there
> is a
> > > way to increase the speed of creation of a parallel socket cluster. The
> > > time spend to include threads increase exponentially with the number of
> > > thread considered and I use of computer with two 8 cores CPU and thus
> > > showing a total of 32 threads in windows 7.
> > > Currently, I use the default parameters (type = "PSOCK"), but is there
> any
> > > fine tuning parameters that I can use to take advantage of this system
> ?
> > >
> > > Thanks in advance for your help !
> > >
> > > Arnaud
> > >
> > > R version 3.0.1 (2013-05-16)
> > > Platform: x86_64-w64-mingw32/x64 (64-bit)
> > >
> > >   [[alternative HTML version deleted]]
> > >
> > > __
> > > R-help@r-project.org mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> >
> >
>
>

[[alternative HTML version deleted]]

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

[R] gmmBoost creating huge 500+ gb vectors

[Luke Miner]

I'm trying to run bGLMM on a large dataset and it is failing after trying to 
create huge vectors with the error:
Error: cannot allocate vector of size 553.7 Gb
The dataset that I'm using is fairly large, around 1.5 million observations. 
The command is as follows:
boost1 <- bGLMM(A~ B + C, rnd = list(D=~1), data = train, family =binomial)
A is binary, B and C are continuous and D is a factor variable with 250,000 
levels that I'm trying to model as a random effect.
A similar model runs fine on lme4. Is my dataset simply much too large for this 
package or am I doing something obviously wrong?
  
[[alternative HTML version deleted]]

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

[R] make system() output invisible on mac

[John Denton]

Hi,

I'm trying to figure out a way to make the output of a call to system() 
invisible on the R mac GUI. Is there a mac equivalent for the calls 
invisible=TRUE, show.output.on.console=FALSE, intern=TRUE for the mac 
environment? Thanks.

~John
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Optimize function in R: unable to find maximum of logistic function

[Shantanu MULLICK]

Hello Everyone,

I want to perform a 1-D optimization by using the optimize() function. I
want to find the maximum value of a "logistic" function. The optimize()
function gives the wrong result.

My code:
f= function (k) {
T_s = 20
result = (2- 2/(1+ exp(-2*T_s*k)))
return(result)
}
optimize(f, c(0, 5), tol = 0.1, maximum= TRUE)

The maximum value for the function happens at k=0, and the maximum value is
1. Yet  the optimise function, says that the maximum value happens at k=
4.9995, and the maximum value is 0.

Thanks in advance!

Warm Regards,
Shantanu

[[alternative HTML version deleted]]

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

[R] Chess Playing Software Written in R

[Rob Bernard]

Have there been attempts to build a chess engine in R?  If so, is it
available anywhere?

I recognize that R really isn’t the right language for a chess engine, but
I was more curious if it had been attempted.

Thank you.

Rob Bernard

-- 
--
Robert N. Bernard
minus...@alumni.princeton.edu

[[alternative HTML version deleted]]

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

[R] high cross-validation errors ( xerror ) in regression tree?

[Lin, Tin-Chi]

Hi All,

This is my first time to seek help from the R-forum (and also conduct a formal 
data-mining analysis). I searched the archive a bit but didn't get responses 
that could totally address my question. Any comments would be highly 
appreciated.

I am using the rpart function to analyze factors that might contribute to 
heightened injury rate; the outcome is a continuous variable. After fitting the 
initial tree and pruning it, the final tree has five terminal nodes with 
cross-validation errors shown below,

CP nsplit rel error xerrorxstd
1 0.139141  0   1.0 1.0033 0.26163
2 0.128314  1   0.86086 1.2752 0.28481
3 0.036021  3   0.60423 1.4315 0.29652
4 0.022675  4   0.56821 1.5142 0.29749
5 0.02  5   0.54554 1.4615 0.28818

My questions are:

 (1) Is this pruned tree even valid? The cross-validation error is exceedingly 
high, well above 1.00

(2) What contributes to the high cross-validation errors (xerror), and why did 
it go up but then down a little bit?

My guess was that the data is quite noisy; therefore, the splitting was pretty 
much based on random noise, resulting in poor prediction. I've found that it 
helps a bit to increase the minimum data points required at splitting and 
terminal nodes (at the expense of the overall R-square of the tree, 
inevitably), but still the problem lingers.

Any thoughts?


The initial tree command and the last five lines of the output is
tree1 <- rpart(overexertion ~, method = "anova", data = data, xval = 10,  
minbucket=4, minsplit=10, cp=0)
Root node error: 502364/347 = 1447.7
n=347 (179 observations deleted due to missingness)

> tree1$cptable[dim(tree1$cptable)[1] - 5:0, ]
 CP nsplit rel error   xerror  xstd
43 9.769565e-05 54 0.2926771 1.641099 0.2626767
44 5.053530e-05 55 0.2925794 1.640780 0.2626771
45 4.314452e-05 56 0.2925288 1.640926 0.2626727
46 2.960797e-05 57 0.2924857 1.640925 0.2626727
47 1.570814e-05 58 0.2924561 1.640941 0.2626724
48 0.00e+00 59 0.2924404 1.640906 0.2626728

The pruning command is
fit9 <- prune(tree1, cp=.02)  #set the cost-complexity parameter at 0.2




Tin-chi Lin
Liberty Mutual Research Institute for Safety
71 Frankland Rd, Hopkinton, MA 01748

Ext: 732-7466
Phone: (508)4970266
Email: tin-chi@libertymutual.com






[[alternative HTML version deleted]]

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

[R] Create Time Lists with a for loop

[Alaios]

Dear all, 
in my code I have written the following list

TimeFramesShort <-list(c(strptime("2011-10-12 10:59:00","%Y-%m-%d 
%H:%M:%S"),strptime("2011-10-13 11:02:00","%Y-%m-%d %H:%M:%S"))  # 
rest items were truncated


I was wondering if could somehow  take the first element of the list and create 
element that are only 5 minutes apart something like
list(c(strptime("2011-10-12 10:59:00","%Y-%m-%d %H:%M:%S"),strptime("2011-10-12 
11:04:00","%Y-%m-%d %H:%M:%S"),
c(strptime("2011-10-12 11:04:00","%Y-%m-%d %H:%M:%S"),strptime("2011-10-12 
11:09:00","%Y-%m-%d %H:%M:%S"),
c(strptime("2011-10-12 11:09:00","%Y-%m-%d %H:%M:%S"),strptime("2011-10-12 
11:13:00","%Y-%m-%d %H:%M:%S")
...
and so on.

The problem I see is that pure string manipulation will not work as after one 
hour that should also make changes in the hours.

Regards
Alex 
[[alternative HTML version deleted]]

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

Re: [R] Split a data.frame

[Christofer Bogaso]

Hi Arun,

Sorry for late reply.

I would like to thank you for your pointer. This is what I wanted. Thanks
for your help.

Thanks and regards,


On Sun, Oct 27, 2013 at 10:33 PM, arun  wrote:

> Hi,
> DF$Col2
> # [1] a e b b a b c e b a c c e e a c d c c e
> #Levels: a b c d e
>
>
> "b" is not found in any of the "Grs".  Also Gr3 is not presnt in DF$Col2
>
> So, I am not sure whether this works for you.
>
> indx <- 1+ 2*DF$Col2 %in% Gr1 + 4*DF$Col2 %in% Gr2 + 8*DF$Col2 %in% Gr3
>  indx <- indx[indx>1]
>  split(DF[DF$Col2%in% c(Gr1,Gr2,Gr3),],indx)
>
> A.K.
>
>
>
>
>
> On Sunday, October 27, 2013 11:19 AM, Christofer Bogaso <
> bogaso.christo...@gmail.com> wrote:
> Hi again,
>
> Let say I have following DF:
>
> DF <- structure(list(Col1 = 1:20, Col2 = structure(c(1L, 5L, 2L, 2L,
> 1L, 2L, 3L, 5L, 2L, 1L, 3L, 3L, 5L, 5L, 1L, 3L, 4L, 3L, 3L, 5L
> ), .Label = c("a", "b", "c", "d", "e"), class = "factor")), .Names =
> c("Col1",
> "Col2"), row.names = c(NA, -20L), class = "data.frame")
>
> DF
>
>
> Now I create 3 groups like:
>
> Gr1 <- c('a', 'c', 'd')
> Gr2 <- c('e')
> Gr3 <- c('f', 'x')
>
>
> My goal is to split DF according to these groups. And to generate NULL (or
> something like that) if a particular group contains no value.
>
> So far I tried to split DF according to split() function. However it looks
> to me like, this split() does not offer this kind of customization.
>
>
> Can someone here help me how to split my data.frame according to this
> criteria?
>
> Thanks and regards,
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>

[[alternative HTML version deleted]]

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

Re: [R] for loop help, repeat a function multiple times

[arun]

Hi,
You may try:
x <- 1:5 
 set.seed(49)
 mat1 <- do.call(rbind,lapply(1:1000,function(y) sample(x,3)))

#or
mat2 <- matrix(0,ncol=3,nrow=1000)
set.seed(49)
 for(i in seq_len(nrow(mat2))) mat2[i,] <- sample(x,3)
all.equal(mat1,mat2)
#[1] TRUE

#or
set.seed(49)

mat3 <- t(replicate(1000,sample(x,3)))
 all.equal(mat2,mat3)
#[1] TRUE



A.K.





Hi, I have written a program that runs a series of calculations using 
random numbers, so that each time that I run the program I get a vector 
with different results. I need to repeat this program multiple times and 
combine all of the vectors into a single data frame or matrix. I am 
fairly certain that a "for loop" is what I need, but after digging 
around online and in multiple books, I am very confused about how to set up the 
loop. Just as an example, suppose that I want to randomly select 3 numbers from 
1 to 5. I could do it using the following two lines: 
x <- 1:5 
y <- sample(x,3) 

now suppose that I want to do this 1000 times, how do I set up a
 loop that will give me 1000 iterations of y (i.e. it will select 1000 
sets of 3 values ranging from 1:5)? In other words, how do I loop y a 
specified number of times? 
Thanks for the help

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

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

[Prof Brian Ripley]

On 28/10/2013 16:19, Arnaud Mosnier wrote:

Hi all,

I am quite new in the world of parallelization and I wonder if there is a
way to increase the speed of creation of a parallel socket cluster. The
time spend to include threads increase exponentially with the number of


It increases linearly in my tests (or a decent OS).  But really if 
parallel computing is worthwhile you will be doing minutes of work on 
each worker process and the startup time will not be signifcant.



thread considered and I use of computer with two 8 cores CPU and thus
showing a total of 32 threads in windows 7.


The first way to speed things up: use a decent OS:  forking clusters is 
much faster.



Currently, I use the default parameters (type = "PSOCK"), but is there any
fine tuning parameters that I can use to take advantage of this system ?

Thanks in advance for your help !

Arnaud

R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)

[[alternative HTML version deleted]]



--
Brian D. Ripley,  rip...@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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

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

[Simon Zehnder]

First,

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

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

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

Simon


On 28 Oct 2013, at 17:51, Arnaud Mosnier  wrote:

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

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

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

[Arnaud Mosnier]

Thanks Simon,

I already read the parallel vignette but I did not found what I wanted.
May be you can be more specific on a part of the document that can provide
me hints !

Arnaud


2013/10/28 Simon Zehnder 

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

[[alternative HTML version deleted]]

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

Re: [R] pie graphs in log scale axis

[David Winsemius]

On Oct 28, 2013, at 6:40 AM, Ortiz, John wrote:



Dear list users,


I'm doing a plot integrating Grid output with Base Graphics output   
(gridBase, Murrell 2012).


My goal is to produce a xy plot where each point is represented by a  
pie.  I could get it using the attach code,

but now I want to change the X axis to log scale.

when I introduce the log="x" parameter in the line:

plot(x, y, xlim=c(0.1, 1.2), log="x", ylim=c(0.1, 1.2), type="n")

I get this warning!

vps <- baseViewports()
Warning message:
In baseViewports() : viewport scales NOT set to user coordinates


This is the code to produce it without introduce log scale yet:

library(grid)
library(gridBase)
x <- c(0.88, 1.00, 0.67, 0.34)
y <- c(0.87, 0.43, 0.24, 0.94)
z <- matrix(runif(4*2), ncol=2)

oldpar <- par(no.readonly=TRUE)

plot(x, y, xlim=c(0.1, 1.2), ylim=c(0.1, 1.2), type="n")

vps <- baseViewports()
pushViewport(vps$inner, vps$figure, vps$plot)

for (i in 1:4) {
pushViewport(viewport(x=unit(x[i], "native"),
y=unit(y[i], "native"),
width=0.1,
height=0.1))
par(plt=gridPLT(), new=TRUE)
pie(z[i,], radius=1, labels=rep("", 2))
popViewport()
}

popViewport(2)
par(oldpar)

Thanks by some advice!!


I haven't figured out the ratinale for this effort but this does  
result in a base log="xy" plot with pies


 library(grid)
library(gridBase)
x <- c(0.88, 1.00, 0.67, 0.34)
y <- c(0.87, 0.43, 0.24, 0.94)
z <- matrix(runif(4*2), ncol=2)

oldpar <- par(no.readonly=TRUE)

plot(x, y, xlim=c(0.01, 1.2), ylim=c(0.01, 1.2), type="n", log="xy")

# Note need to modify limits

vps <- baseViewports()
pushViewport(vps$inner, vps$figure, vps$plot)

for (i in 1:4) {
pushViewport(viewport(x=unit(log(x[i]), "native"),
# this seemed the sensible strategy
y=unit(log(y[i]), "native"),
# Didn't think log(units(...)) would be as successful
width=0.1,
height=0.1))
par(plt=gridPLT(), new=TRUE)
pie(z[i,], radius=1, labels=rep("", 2))
popViewport()
}

popViewport(2)
par(oldpar)



Regards,


John Ortiz

Geologist
Smithsonian Tropical Research Institute Panama

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


David Winsemius, MD
Alameda, CA, USA

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

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

[Simon Zehnder]

See library(help = "parallel”)


On 28 Oct 2013, at 17:19, Arnaud Mosnier  wrote:

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

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

Re: [R] linear line in x, y plot

[David Winsemius]

On Oct 28, 2013, at 8:30 AM, Ahmed Attia wrote:


Hi

I have a question about drawing a linear line in x, y plot. I usually
use the following code, but for this time the x values are to small
(-0.08 to -0.02)


That is not the problem.



I wrote the following code, but r does not draw the line. However, it
does not give an error when it takes the code.


reg1<- lm(CWSI~NWI, data=Ahmed)

NWI <- seq(-0.08, -0.02, len = -0.02)


>  seq(-0.08, -0.02, len = -0.02)
integer(0)

The length argument to seq doesn't make any sense. Please review:

?seq

--  
David

lines(NWI, predict(reg1, list(NWI =
NWI)),xlim=c(-0.08,-0.02),ylim=c(0,1),pch=1,col=1,lwd=2, lty=1)

When I wrote the following code

abline(reg1,pch=2,col=2,lwd=2,  
lty=2,xlim=c(-0.06,-0.02),ylim=c(0.3,0.8))


the line shows up, but I can not control the xlim or ylim. It bascally
goes across the figure

I appreciate your help
--
Ahmed M. Attia


Research Assistant
Dept. Of Soil&Crop Sciences
Texas A&M University
ahmed.at...@ag.tamu.edu
Cell phone: 001-979-248-5215

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


David Winsemius, MD
Alameda, CA, USA

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

[R] Automatically Remove Aliased Terms from a Model

[Thaler,Thorn,LAUSANNE,Applied Mathematics]

Dear all,

I am trying to implement a function which removes aliased terms from a model. 
The challenge I am facing is that with "alias" I get the aliased coefficients 
of the model, which I have to translate into the terms from the model formula. 
What I have tried so far:

--8<--
d <- expand.grid(a = 0:1, b=0:1)
d$c <- (d$a + d$b)  %% 2
d$y <- rnorm(4)
d <- within(d, {a <- factor(a); b <- factor(b); c <- factor(c)})
l <- lm(y ~ a * b + c, d)

removeAliased <- function(mod) {
  ## Retrieve all terms in the model
  X <- attr(mod$terms, "term.label")
  ## Get the aliased coefficients  
  rn <- rownames(alias(mod)$Complete)
  ## remove factor levels from coefficient names to retrieve the terms
  regex.base <- unique(unlist(lapply(mod$model[, sapply(mod$model, is.factor)], 
levels)))
  aliased <- gsub(paste(regex.base, "$", sep = "", collapse = "|"),  "", 
gsub(paste(regex.base, ":", sep = "", collapse = "|"), ":", rn))
  uF <- formula(paste(". ~ .", paste(aliased, collapse = "-"), sep = "-"))
  update(mod, uF)
}

removeAliased(l)
-->8--

This function works in principle, but this workaround with removing the factor 
levels is just, well, a workaround which could cause problems in some 
circumstances (when the name of a level matches the end of another variable, 
when I use a different contrast and R names the coefficients differently etc. - 
and I am not sure which other cases I am overlooking).

So my question is whether there are some more intelligent ways of doing what I 
want to achieve? Is there a function to translate a coefficient of a LM back to 
the term, something like:

termFromCoef("a1") ## a1
termFromCoef("a1:b1") ## a:b

With this I could simply translate the rownames from alias into the terms 
needed for the model update.

Thanks for your help.

Kind Regards,

Thorn Thaler 
NRC Lausanne
Applied Mathematics

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

[R] speed of makeCluster (package parallel)

[Arnaud Mosnier]

Hi all,

I am quite new in the world of parallelization and I wonder if there is a
way to increase the speed of creation of a parallel socket cluster. The
time spend to include threads increase exponentially with the number of
thread considered and I use of computer with two 8 cores CPU and thus
showing a total of 32 threads in windows 7.
Currently, I use the default parameters (type = "PSOCK"), but is there any
fine tuning parameters that I can use to take advantage of this system ?

Thanks in advance for your help !

Arnaud

R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)

[[alternative HTML version deleted]]

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

Re: [R] Error: Data is not atomic

[Berend Hasselman]

On 28-10-2013, at 17:16, Baro  wrote:

> thank you so much. I will write my questions as you mentioned
> 
> 

Final advice. 
Please reply to the list only. Everybody can then see your reply and respond  
if necessary.

Forwarded to the list.

Berend

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

Re: [R] Error: Data is not atomic

[Berend Hasselman]

On 28-10-2013, at 17:01, Baro  wrote:

> Hi experts,
> 
> I want to user haar wavelet transform. If I am using this simple command,
> it works nice:
> 
> k<-c(1,2,3,4,5,6,7,8)
> ywd<-wd(k,filter.number=1,family="DaubExPhase")
> 
> but if the K is a list like this:
> 
> *[[1]]*
> *[1] 401*
> *
> *
> *[[2]]*
> *[1] 481*
> *
> *
> *[[3]]*
> *[1] 480*
> *
> *
> *[[4]]*
> *[1] 482*
> *
> *
> *[[5]]*
> *[1] 395*
> *...*
> *
> *
> it doesnt work and I get this error:
> 
> Error in wd(k, filter.number = 1, family = "DaubExPhase") :
> * Data is not atomic*
> *
> *
> how can I solve this problem?how can I convert this list to a normal list?
> 

Please do not post in html as requested by the posting guide.
And don't post the same question twice within a short period of time.
Please do NOT put code or Routput in bold or whatever. Plain text only.
In plain text that is displayed as:

*[1] 395*

which is obviously nonsense.

Have a look at unlist().

Berend


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

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

Re: [R] linear line in x, y plot

[Adams, Jean]

A few problems ...

This statement doesn't make sense.
seq(-0.08, -0.02, len = -0.02)
Perhaps you meant
seq(-0.08, -0.02, by = 0.02)

The xlim= and ylim= are arguments to higher level plot functions, like
plot(), and won't work for functions lines() or abline().

Are you trying to limit the range of the predicted line on the graph?  If
so, perhaps the following simple example will help.

# randomly generate some fake data
CWSI <- rnorm(20, sd=0.1)
NWI <- rnorm(20, sd=0.1)
# linear regression
reg1 <- lm(CWSI~NWI)
# show predicted line for small segment
newNWI <- seq(-0.08, -0.02, by = 0.02)
plot(NWI, CWSI)
lines(newNWI, predict(reg1, list(NWI = newNWI)))

Jean



On Mon, Oct 28, 2013 at 10:30 AM, Ahmed Attia  wrote:

> Hi
>
> I have a question about drawing a linear line in x, y plot. I usually
> use the following code, but for this time the x values are to small
> (-0.08 to -0.02)
>
> I wrote the following code, but r does not draw the line. However, it
> does not give an error when it takes the code.
>
>
> reg1<- lm(CWSI~NWI, data=Ahmed)
>
> NWI <- seq(-0.08, -0.02, len = -0.02)
> lines(NWI, predict(reg1, list(NWI =
> NWI)),xlim=c(-0.08,-0.02),ylim=c(0,1),pch=1,col=1,lwd=2, lty=1)
>
> When I wrote the following code
>
> abline(reg1,pch=2,col=2,lwd=2, lty=2,xlim=c(-0.06,-0.02),ylim=c(0.3,0.8))
>
> the line shows up, but I can not control the xlim or ylim. It bascally
> goes across the figure
>
> I appreciate your help
> --
> Ahmed M. Attia
>
>
> Research Assistant
> Dept. Of Soil&Crop Sciences
> Texas A&M University
> ahmed.at...@ag.tamu.edu
> Cell phone: 001-979-248-5215
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

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

[R] Error: Data is not atomic

[Baro]

Hi experts,

I want to user haar wavelet transform. If I am using this simple command,
it works nice:

k<-c(1,2,3,4,5,6,7,8)
ywd<-wd(k,filter.number=1,family="DaubExPhase")

but if the K is a list like this:

*[[1]]*
*[1] 401*
*
*
*[[2]]*
*[1] 481*
*
*
*[[3]]*
*[1] 480*
*
*
*[[4]]*
*[1] 482*
*
*
*[[5]]*
*[1] 395*
*...*
*
*
it doesnt work and I get this error:

Error in wd(k, filter.number = 1, family = "DaubExPhase") :
 * Data is not atomic*
*
*
how can I solve this problem?how can I convert this list to a normal list?

[[alternative HTML version deleted]]

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

Re: [R] Code Book from SPSS Data

[Frank Harrell]

Not certain about .por but this works with ordinary SPSS files:

require(Hmisc)
dat <- spss.get(...)  # gets variable labels, etc.
contents(dat)
html(contents(dat), ...)

The last command produces a hyperlinked data dictionary, e.g., for each 
variable the number of levels is given and you click on that number to 
see the levels.  Variables having the same levels are combined in the 
latter part.


Frank

--
Frank E Harrell Jr Professor and Chairman  School of Medicine
   Department of Biostatistics Vanderbilt University

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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: Undeliverable mail

[hill0093 University of Minnesota]

-- Forwarded message --
From: 
Date: Mon, Oct 28, 2013 at 10:53 AM
Subject: Undeliverable mail
To: hill0...@umn.edu
Cc: postmas...@vs-m.tc.umn.edu


Your message was not delivered to the following recipients:

   r-h...@r-proj.org: No such host

Original-Recipient: rfc822;r-h...@r-proj.org
Final-Recipient: rfc822;r-h...@r-proj.org
Action: failed
Status: 5.1.2


-- Forwarded message --
From: hill0093 University of Minnesota 
To: r-h...@r-proj.org
Cc:
Date: Mon, 28 Oct 2013 10:53:16 -0500
Subject: Fwd: Labeling with a function of the scale value of an axis


-- Forwarded message --
From: hill0093 University of Minnesota 
Date: Fri, Oct 25, 2013 at 12:05 PM
Subject: Labeling with a function of the scale value of an axis
To: r-c...@r-project.org


I am sending this to your email address because
it is listed as the maintainer of graphics containing axis.
Please deliver it to an appropriate person.

This is an inappropriate R run:
> xvalue<-c(5.2,1.3,9.7,2.8,8.1,4.7,6.6,7.4)
> yvalue<-c(9,3,4,7,2,5,3,6)
> plot(xvalue,yvalue)
> axis(1,at=NULL,labels=1/xvalue,digits=5)
Error in axis(1, at = NULL, labels = 1/xvalue, digits = 5) :
  'labels' is supplied and not 'at'
In addition: Warning message:
In axis(1, at = NULL, labels = 1/xvalue, digits = 5) :
  "digits" is not a graphical parameter
>
I want a way (not changing the axis scale) of
labeling (labels = the numbers at the tick marks)
the axis with a mathematical function of what
would be the scaled value at the tick mark.
A commonly useful function is label=1/scaledValue(not 1/xvalue),
but there are are other functions (user supplied) like
changing units or format of the labels.
For quick-plots, it would just use positions
automatically computed for the scaled value,
and it could (not necessarily) automatically
adjust depending on the new sizes of the labels.
There might also be parameters specifying
the label's total number of digits and/or
total width and/or number of digits after the point.
It should be fairly easy for the writer of the
R axis code to supply the extremely useful labels of
at least the 1/value at the tick mark originally for
value as a useful starter. Where the user wants the
labels would be separate as it is now.
I cannot modify axis or any other R code because
of the complexity of the R (foreign) language which
I cannor understand. Do I instead have to write
subroutines that do this job using the current
axis routine?
Dewayne Hillman
hill0...@umn.edu

[[alternative HTML version deleted]]

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

Re: [R] Revo R for Arima Implementation

[Jeff Newmiller]

You still don't seem to get the hint that support that is SPECIFIC to 
Revolution R is off- topic on this mailing list. Since you have identified your 
subject so clearly, we can only safely respond that we cannot answer you.

If you change your mind and decide to ask what can be done in R (which should 
also work in Revolution R but may not be optimal), you should also read the 
Posting Guide and post a reproducible example that demonstrates where you are 
in your analysis. You should also use the RSiteSearch function to look up 
packages that might address your needs, since this is not a statistics theory 
discussion forum either.
---
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.

Anindita Chattopadhyay  wrote:
>Hi Rolf,
>
>Thanks for the response. I have re-phrased the problem. Hope this will
>help.
>
>We're working on a project where our model tries to predict the value
>of bookings as a time series of past bookings as well as some external
>flags.
>The data was found to be stationary using the Dicky-Fuller test.
>We have only a non-seasonal component for the AR and MA terms (small p
>& q only).
>No differencing was done.
>
>Then, we used values of p = 1, 7, 8 and q = 1, 7.
>The main issue is that the p values are additive and q values are
>multiplicative.
>We need to understand how we can implement this in Revo R.
>
>Just to add, since we have external flags , we cannot make use of
>SARIMA function.
>Please let me know if any other questions.
>
>Thank you in advance!
>
>Regards,
>Anindita Chattopadhyay | +919886800606 | www.mu-sigma.com |
>
>-Original Message-
>From: Rolf Turner [mailto:r.tur...@auckland.ac.nz]
>Sent: Saturday, October 26, 2013 1:56 AM
>To: Anindita Chattopadhyay
>Cc: r-help@r-project.org; Harish K
>Subject: Re: [R] Revo R for Arima Implementation
>
>
>Your question is pretty well totally opaque to me.  Describe the model
>you want to fit in mathematical terms, rather than referring to The
>Package That Must Not Be Named.  This is the ***R***  list.
>
>It is possible that you might want to make use of the "seasonal"
>argument to the
>arima() function.
>
> cheers,
>
> Rolf Turner
>
>On 10/26/13 01:19, Anindita Chattopadhyay wrote:
>> Hello There,
>>
>> We have used ARIMA(multiplicative MA and additive AR) model in SAS to
>come to our results. Please let me know how could that be implemented
>in R.
>>
>> Like,  In SAS we can pass the variable  as AR= (1,7) and MA=(1)(7)
>which is additive and multiplicative respectively.
>>
>> In R we have the option as : order = c(0, 0, 0) which is (p,d,q).Is
>there any specific way/code such that we specify the additive  AR and
>multiplicative MA simultaneously?
>>
>This email message may contain proprietary, private and confidential
>information. The information transmitted is intended only for the
>person(s) or entities to which it is addressed. Any review,
>retransmission, dissemination or other use of, or taking of any action
>in reliance upon, this information by persons or entities other than
>the intended recipient is prohibited and may be illegal. If you
>received this in error, please contact the sender and delete the
>message from your system. Mu Sigma takes all reasonable steps to ensure
>that its electronic communications are free from viruses. However,
>given Internet accessibility, the Company cannot accept liability for
>any virus introduced by this e-mail or any attachment and you are
>advised to use up-to-date virus checking software.
>
>__
>R-help@r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide
>http://www.R-project.org/posting-guide.html
>and provide commented, minimal, self-contained, reproducible code.

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

Re: [R] How to extract residuals from multiple regressions from a loop

[Adams, Jean]

Ben,

It helps us help you if you provide a simple example with code that anyone
can run.  When I created a simplified version of you situation, the NAs in
the response variable yielded NAs in the residuals just as the help for
na.exclude indicates.

?na.exclude

Can you replicate the problem (and the error message) you are having using
reproducible code?

Jean


# number of observations
n1 <- 10
# number of response variables
n2 <- 5
# randomly generate the independent variable
x <- rnorm(n1)
# create a matrix to store the residuals in
keep <- matrix(NA, nrow=n1, ncol=n2)

# fit a linear model for each response variable
for(i in 1:n2) {
 # randomly generate the response
y <- rnorm(n1)
# randomly put in some missing values
 y[sample(n1, 2)] <- NA
# fit the data
fit <- lm(y ~ x, na.action=na.exclude)
 # keep the residuals
keep[, i] <- resid(fit)
}
# look at the resulting matrix of residuals
keep

[,1][,2]   [,3][,4]   [,5]
 [1,] -1.1435507  0.29915778 -0.4593465 -0.61984029  0.8691960
 [2,]  1.3410409  0.51701634 NA  0.65691397 NA
 [3,]  1.0896517  NA  0.2239847 -0.63233644 -1.0831747
 [4,] NA -1.80344171 -1.2984848 -0.26543679  0.4486482
 [5,] -2.0836253  1.05313477 NA -0.02031142 -0.1059559
 [6,]  1.4498942  0.61520388 -0.2172015 -0.90647457  0.4935462
 [7,] -0.6265764 -0.08396366  0.5153020  NA -0.2501000
 [8,] NA  NA  0.9337658 -0.10289794 NA
 [9,] -0.4217575 -0.42633169  0.1070141  1.89038347 -0.6588342
[10,]  0.3949231 -0.17077571  0.1949662  NA  0.2866744



On Mon, Oct 28, 2013 at 8:19 AM, Ben Ammar  wrote:

>
>Dear all
>
>I've  got  the following problem, I want to extract the residuals from
>regression loops. The problem here is that some columns include NA's at
> the
>beginning and end (i.e. each time series of stocks starts at different
>points  in  time and ends at different points in time). When I want to
>transfer these residuals into a matrix to determine the residual
> matrix, I
>get the error message ("number of items to replace is not a multiple of
>replacement length"). I tried it with na.action=na.exclude but that
> doesn't
>work because that command doesn't actually change the vector length.
> With a
>loop I came this far:
>Number of stocks is 50 and maximum time period is 258 months:
>
>for (i in 1:50) {CAPM.res[,i] <- residuals(lm(timeseries[,i]~exc.mkt),
>na.action=na.exclude)}
>
>as I said it doesn't work because of the different column length in the
>matrix "timeseries". So right now I'm doing kind of manually which works
>perfectly but is quite intensive and looks like that:
>test.1 <- lm(timeseries[,1]~exc.mkt, na.action=na.exclude)
>residual.test.1 <- residuals(test.1)
>CAPM.res[,1] <- residual.Life.1
>
>test.2 <- lm(timeseries[,2]~exc.mkt, na.action=na.exclude)
>residual.test.2 <- residuals(test.2)
>CAPM.res[,2] <- residual.test.2
>
>and so on for the remaining 49 stocks. When I look at that I
> obviously
>see that this must be done with a loop but in the end I can't put in the
>matrix because of the different lengths. So far I got this:
>test<-matrix(0,50,258)
>residual.test<-matrix(0,50,258)
>for (i in 1:50) {lm(timeseries[,i]~exc.mkt, na.action=na.exclude)
>{residual.test[i] <- residuals(test[i])
>{CAPM.res[,i] <- residual.test[i]
>}}}
>
>but here I get the error message: "Error: $ operator is invalid for
> atomic
>vectors"
>and I don't think "test"  and "residual.test" is defined correctly
> because I
>don't know where to look for the residuals.
>
>Does anyone have an idea how to extract the residuals and put them in a
>258x50 matrix?
>Any help would be very much appreciated!
>
>Cheers,
>Ben
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

[[alternative HTML version deleted]]

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

[R] Error:Data is not atomic

[Babak Bastan]

Hi experts,

I want to user haar wavelet transform. If I am using this simple command,
it works nice:

k<-c(1,2,3,4,5,6,7,8)
ywd<-wd(k,filter.number=1,family="DaubExPhase")

but if the K is a list like this:

*[[1]]*
*[1] 401*
*
*
*[[2]]*
*[1] 481*
*
*
*[[3]]*
*[1] 480*
*
*
*[[4]]*
*[1] 482*
*
*
*[[5]]*
*[1] 395*
*...*
*
*
it doesnt work and I get this error:

Error in wd(k, filter.number = 1, family = "DaubExPhase") :
 * Data is not atomic*
*
*
how can I solve this problem?how can I convert this list to a normal list?

[[alternative HTML version deleted]]

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

[R] linear line in x, y plot

[Ahmed Attia]

Hi

I have a question about drawing a linear line in x, y plot. I usually
use the following code, but for this time the x values are to small
(-0.08 to -0.02)

I wrote the following code, but r does not draw the line. However, it
does not give an error when it takes the code.


reg1<- lm(CWSI~NWI, data=Ahmed)

NWI <- seq(-0.08, -0.02, len = -0.02)
lines(NWI, predict(reg1, list(NWI =
NWI)),xlim=c(-0.08,-0.02),ylim=c(0,1),pch=1,col=1,lwd=2, lty=1)

When I wrote the following code

abline(reg1,pch=2,col=2,lwd=2, lty=2,xlim=c(-0.06,-0.02),ylim=c(0.3,0.8))

the line shows up, but I can not control the xlim or ylim. It bascally
goes across the figure

I appreciate your help
-- 
Ahmed M. Attia


Research Assistant
Dept. Of Soil&Crop Sciences
Texas A&M University
ahmed.at...@ag.tamu.edu
Cell phone: 001-979-248-5215

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 extract residuals from multiple regressions from a loop

[William Dunlap]

Can you trim down your example to a size where you can show us
the data (using dump() or dput()) and the commands you used so
one can just copy it from your mail and paste it into R to reproduce
the problem?

I don't see your problem when I made up data similar to what you described:

> timeseries <- ts(cbind(c(NA, 2:9, NA), c(1,NA,3,NA,5:10), c(rep(NA,9), 10)))
> exc.mkt <- log(1:10)
> for(i in 1:ncol(timeseries)) 
> print(length(residuals(lm(timeseries[,i]~exc.mkt, na.action=na.exclude
[1] 10
[1] 10
[1] 10
> # without na.action=na.exclude residuals do vary in length
> for(i in 1:ncol(timeseries)) 
> print(length(residuals(lm(timeseries[,i]~exc.mkt
[1] 8
[1] 8
[1] 1

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
> Behalf
> Of Ben Ammar
> Sent: Monday, October 28, 2013 6:20 AM
> To: r-help@r-project.org
> Subject: [R] How to extract residuals from multiple regressions from a loop
> 
> 
>Dear all
> 
>I've  got  the following problem, I want to extract the residuals from
>regression loops. The problem here is that some columns include NA's at the
>beginning and end (i.e. each time series of stocks starts at different
>points  in  time and ends at different points in time). When I want to
>transfer these residuals into a matrix to determine the residual matrix, I
>get the error message ("number of items to replace is not a multiple of
>replacement length"). I tried it with na.action=na.exclude but that doesn't
>work because that command doesn't actually change the vector length. With a
>loop I came this far:
>Number of stocks is 50 and maximum time period is 258 months:
> 
>for (i in 1:50) {CAPM.res[,i] <- residuals(lm(timeseries[,i]~exc.mkt),
>na.action=na.exclude)}
> 
>as I said it doesn't work because of the different column length in the
>matrix "timeseries". So right now I'm doing kind of manually which works
>perfectly but is quite intensive and looks like that:
>test.1 <- lm(timeseries[,1]~exc.mkt, na.action=na.exclude)
>residual.test.1 <- residuals(test.1)
>CAPM.res[,1] <- residual.Life.1
> 
>test.2 <- lm(timeseries[,2]~exc.mkt, na.action=na.exclude)
>residual.test.2 <- residuals(test.2)
>CAPM.res[,2] <- residual.test.2
> 
>and so on for the remaining 49 stocks. When I look at that I obviously
>see that this must be done with a loop but in the end I can't put in the
>matrix because of the different lengths. So far I got this:
>test<-matrix(0,50,258)
>residual.test<-matrix(0,50,258)
>for (i in 1:50) {lm(timeseries[,i]~exc.mkt, na.action=na.exclude)
>{residual.test[i] <- residuals(test[i])
>{CAPM.res[,i] <- residual.test[i]
>}}}
> 
>but here I get the error message: "Error: $ operator is invalid for atomic
>vectors"
>and I don't think "test"  and "residual.test" is defined correctly because 
> I
>don't know where to look for the residuals.
> 
>Does anyone have an idea how to extract the residuals and put them in a
>258x50 matrix?
>Any help would be very much appreciated!
> 
>Cheers,
>Ben
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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

Re: [R] Optimization failed in fitdistr (Weibull distribution)

[Berend Hasselman]

On 28-10-2013, at 16:07, Rui Barradas  wrote:

> Hello,
> 
> I can't reproduce your error:
> 
> windfreq <-
> c(1351L, 2147L, 3317L, 4378L, 5527L, 6667L, 7865L, 8970L, 9987L,
> 10907L, 11905L, 12642L, 131000L, 14983L, 15847L, 16842L, 17757L,
> 18698L, 19632L, 20626L, 21599L, 22529L, 23325L, 24391L, 25356L,
> 26267L, 27230L, 28223L, 29190L, 30142L, 31124L, 32104L, 3397L,
> 3437L, 3562L, 3646L, 3742L, 3824L, 399L, 4013L, 419L, 425L, 432L
> 
> library(MASS)
> 
> fitdistr(windfreq, "weibull")

You seem to have collapsed the two columns of data into a single column.
That's why it works for you.

The OP should use the dput function to show the data.
If I take the data as presented I get the same error (but I don't use sep="," 
since the data displayed don't contain a ,).
And show the actual code used such as library(MASS).

Berend

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

Re: [R] Optimization failed in fitdistr (Weibull distribution)

[Rui Barradas]

Hello,

I can't reproduce your error:

windfreq <-
c(1351L, 2147L, 3317L, 4378L, 5527L, 6667L, 7865L, 8970L, 9987L,
10907L, 11905L, 12642L, 131000L, 14983L, 15847L, 16842L, 17757L,
18698L, 19632L, 20626L, 21599L, 22529L, 23325L, 24391L, 25356L,
26267L, 27230L, 28223L, 29190L, 30142L, 31124L, 32104L, 3397L,
3437L, 3562L, 3646L, 3742L, 3824L, 399L, 4013L, 419L, 425L, 432L

library(MASS)

fitdistr(windfreq, "weibull")


Hope this helps,

Rui Barradas

Em 28-10-2013 12:07, kmmoon100 escreveu:

Hello everyone,

This is Kangmin.

I am trying to produce shape and scale of my wind data. My data is based on
wind speed frequency with 1km/hr increment. data is described below.

Windspeed (km/h)Frequency
1   351
2   147
3   317
4   378
5   527
6   667
7   865
8   970
9   987
10  907
11  905
12  642
13  1000
14  983
15  847
16  842
17  757
18  698
19  632
20  626
21  599
22  529
23  325
24  391
25  356
26  267
27  230
28  223
29  190
30  142
31  124
32  104
33  97
34  37
35  62
36  46
37  42
38  24
39  9
40  13
41  9
42  5
43  2

R codes to calculate shape and scale are described below:

Pine.windfrequency.4weeks<-read.table("C:/Users/kmoon/Documents/Pine_frequency_4weeks.csv",header=TRUE,sep=",")
fitdistr(Pine.windfrequency.4weeks$Frequency, densfun="weibull")

I have got an error message when I was using 'fitdistr' function

"Error in fitdistr(Pine.windfrequency.4weeks$Frequency, densfun = "weibull")
:
   optimization failed"

Please help me calculating shape and scale of weibull distribution.

And please understand that I am not an user familiar with R program but I am
really trying to make my analysis work on R!

Thank you!!!

Kangmin.



--
View this message in context: 
http://r.789695.n4.nabble.com/Optimization-failed-in-fitdistr-Weibull-distribution-tp4679167.html
Sent from the R help mailing list archive at Nabble.com.

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



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

[R] Countour for netCDF-like meshes

[Julio Sergio Santana]

Hi,
A model gives me as output a netCDF file which I read with ncdf package. The 
output consists of three similar matrices: one contains the variable value, 
and the other two, the longitude and latitude coordinates. I need to do a 
contour graph with such information, however, the R contour function 
expects, 
two vectors X, Y, and a matrix Z, as its input. X and Y, being the same 
length as the number of rows and columns of Z matrix, respectively. This 
implies that the Z values in the same row, all have the same abscissa, and 
the Z values in the same colum, all have the same ordinate.

I know, in systems  similar to matlab, ferret, ncl, and the like, you just 
give the three matrices, and the system performs all the necessary 
interpolations to display the contour. Do you know if there is a way to do 
the same in R?

Though not exactly an output in a netCDF, I'm presenting here a very 
simplified example of the kind of output I need to plot.

  val <- matrix(c(rep(15,5),rep(c(15:17,16,15),2),rep(15,5)),nrow=4,byrow=T)
  v1 <- 1:5
  v2 <- 10:13
  xx <- cbind(v2, v2+.1, v2+.2, v2+.1, v2)
  yy <- rbind(v1, v1+0.1, v1+0.2, v1+0.1)

val: is the matrix with the values I want to plot.
xx: is the corresponding abscissa values.
yy: is the corresponging ordinate values.


Thanks,

  -Sergio.

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

Re: [R] indicating significant differences in boxplots

[Johannes Radinger]

Hi,

I'd like to follow up an older discussion on adding significance stars to
boxplots.
(
http://r.789695.n4.nabble.com/indicating-significant-differences-in-boxplots-td862335.html
)
As suggested by Jim Lemon, there is following solution for a single boxplot:

boxplot(count ~ spray, data = InsectSprays, col = "lightgray")

#and we'll assume that the second and third boxplots require a star:
par(xpd=TRUE)
yrange<-par("usr")[3:4]
ypos<-yrange[2]+diff(yrange)/40
segments(2,ypos,3,ypos)
text(2.5,ypos+diff(yrange)/40,"*",cex=2)
par(xpd=FALSE)

But what if one wants to plot two boxplots (e.g. with each only two groups
which are significantly different) and
add such stars to each of the boxplots. I usually use par(mfrow(1,2)) to
plot all boxplots in a row.
How can I also add segments and stars on top of each plot?

/johannes

[[alternative HTML version deleted]]

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

[R] "Optimization fail" error from fitdistr (Weibull distribution)

[kmmoon100]

Hello everyone, 

This is Kangmin. 

I am trying to produce shape and scale of my wind data. My data is based on
wind speed frequency with 1km/hr increment. data is described below.
 
Windspeed (km/h)Frequency 
1 351 
2 147 
3 317 
4 378 
5 527 
6 667 
7 865 
8 970 
9 987 
10 907 
11 905 
12 642 
13 1000 
14 983 
15 847 
16 842 
17 757 
18 698 
19 632 
20 626 
21 599 
22 529 
23 325 
24 391 
25 356 
26 267 
27 230 
28 223 
29 190 
30 142 
31 124 
32 104 
33 97 
34 37 
35 62 
36 46 
37 42 
38 24 
39 9 
40 13 
41 9 
42 5 
43 2 

R codes to calculate shape and scale are described below: 

Pine.windfrequency.4weeks<-read.table("C:/Users/kmoon/Documents/Pine_frequency_4weeks.csv",header=TRUE,sep=",")
 fitdistr(Pine.windfrequency.4weeks$Frequency, densfun="weibull") 

I have got an error message when I was using 'fitdistr' function 

"Error in fitdistr(Pine.windfrequency.4weeks$Frequency, densfun = "weibull")
: 
  optimization failed" 

Please help me calculating shape and scale of weibull distribution. 

And please understand that I am not an user familiar with R program but I am
really trying to make my analysis work on R!
 
Thank you!!! 

Kangmin. 



--
View this message in context: 
http://r.789695.n4.nabble.com/Optimization-fail-error-from-fitdistr-Weibull-distribution-tp4679169.html
Sent from the R help mailing list archive at Nabble.com.

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

Re: [R] Problems with lme random slope+intercept model

[anord]

Dear Ben, 
Thank you, I agree. I have forwarded this to the r-sig-mixed-models list.

Best,
Andreas



--
View this message in context: 
http://r.789695.n4.nabble.com/Problems-with-lme-random-slope-intercept-model-tp4679104p4679168.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] Optimization failed in fitdistr (Weibull distribution)

[kmmoon100]

Hello everyone,

This is Kangmin.

I am trying to produce shape and scale of my wind data. My data is based on
wind speed frequency with 1km/hr increment. data is described below.

Windspeed (km/h)Frequency
1   351
2   147
3   317
4   378
5   527
6   667
7   865
8   970
9   987
10  907
11  905
12  642
13  1000
14  983
15  847
16  842
17  757
18  698
19  632
20  626
21  599
22  529
23  325
24  391
25  356
26  267
27  230
28  223
29  190
30  142
31  124
32  104
33  97
34  37
35  62
36  46
37  42
38  24
39  9
40  13
41  9
42  5
43  2

R codes to calculate shape and scale are described below:

Pine.windfrequency.4weeks<-read.table("C:/Users/kmoon/Documents/Pine_frequency_4weeks.csv",header=TRUE,sep=",")
fitdistr(Pine.windfrequency.4weeks$Frequency, densfun="weibull")

I have got an error message when I was using 'fitdistr' function

"Error in fitdistr(Pine.windfrequency.4weeks$Frequency, densfun = "weibull")
: 
  optimization failed"

Please help me calculating shape and scale of weibull distribution.

And please understand that I am not an user familiar with R program but I am
really trying to make my analysis work on R!

Thank you!!!

Kangmin.



--
View this message in context: 
http://r.789695.n4.nabble.com/Optimization-failed-in-fitdistr-Weibull-distribution-tp4679167.html
Sent from the R help mailing list archive at Nabble.com.

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

Re: [R] coxph: how to define interaction terms?

[rm]

The output is as follows. My question is how to include only the interaction
terms where x is equal to "Maintained".

(x has two possible values "Maintained" and "Nonmaintained".)





--
View this message in context: 
http://r.789695.n4.nabble.com/coxph-how-to-define-interaction-terms-tp4679162p4679176.html
Sent from the R help mailing list archive at Nabble.com.

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

Re: [R] R-help Digest, Vol 128, Issue 30

[Lorenz, David]

Pavlos,
  There are several ways to evaluate how well new data fit an old
regression.Part of the answer depends on what you are concerned about. For
example, if you are concerned about bias, you can test whether the mean of
the new data is within the expected range of the mean of that many new
values. The equations for these prediction intervals should be in good
texts on linear regression.
Dave

Date: Sun, 27 Oct 2013 13:36:12 +0200
From: Pavlos Pavlidis 
To: r-help 
Subject: Re: [R] how well *new data* fit a pre-computed model
Message-ID:

Content-Type: text/plain

Here is a link to a plot that illustrates the question:
bio.lmu.de/~pavlidis/pg1.pdf

the question is how to evaluate whether the blue points fit the curve well
enough. The curve has been produced from the black points

best
pavlos


On Sun, Oct 27, 2013 at 1:30 PM, Pavlos Pavlidis wrote:

> Hi all,
> I have fitted polynomial models to a bunch of data (time-course analysis).
> The experiment is "the expression value of gene A under condition K over
> time". The data points that have been used to fit the model are about 200
> (dataset A). Furthermore I have a few data (dataset B; about 10 points)
for
> "the expression values of gene A under condition G over time". The
question
> is:
>
> how can I evaluate how well the dataset B fits the model generated by
> dataset A?
>
> kind regards,
> pavlos
>
> --
>
> Pavlos Pavlidis, PhD
>
> Foundation for Research and Technology - Hellas
> Institute of Molecular Biology and Biotechnology
> Íikolaou Plastira 100, Vassilika Vouton
> GR - 711 10, Heraklion, Crete, Greece
>



--

Pavlos Pavlidis, PhD

Foundation for Research and Technology - Hellas
Institute of Molecular Biology and Biotechnology
Íikolaou Plastira 100, Vassilika Vouton
GR - 711 10, Heraklion, Crete, Greece

[[alternative HTML version deleted]]



--

[[alternative HTML version deleted]]

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

[R] How to extract residuals from multiple regressions from a loop

[Ben Ammar]

   Dear all

   I've  got  the following problem, I want to extract the residuals from
   regression loops. The problem here is that some columns include NA's at the
   beginning and end (i.e. each time series of stocks starts at different
   points  in  time and ends at different points in time). When I want to
   transfer these residuals into a matrix to determine the residual matrix, I
   get the error message ("number of items to replace is not a multiple of
   replacement length"). I tried it with na.action=na.exclude but that doesn't
   work because that command doesn't actually change the vector length. With a
   loop I came this far:
   Number of stocks is 50 and maximum time period is 258 months:

   for (i in 1:50) {CAPM.res[,i] <- residuals(lm(timeseries[,i]~exc.mkt),
   na.action=na.exclude)}

   as I said it doesn't work because of the different column length in the
   matrix "timeseries". So right now I'm doing kind of manually which works
   perfectly but is quite intensive and looks like that:
   test.1 <- lm(timeseries[,1]~exc.mkt, na.action=na.exclude)
   residual.test.1 <- residuals(test.1)
   CAPM.res[,1] <- residual.Life.1

   test.2 <- lm(timeseries[,2]~exc.mkt, na.action=na.exclude)
   residual.test.2 <- residuals(test.2)
   CAPM.res[,2] <- residual.test.2

   and so on for the remaining 49 stocks. When I look at that I obviously
   see that this must be done with a loop but in the end I can't put in the
   matrix because of the different lengths. So far I got this:
   test<-matrix(0,50,258)
   residual.test<-matrix(0,50,258)
   for (i in 1:50) {lm(timeseries[,i]~exc.mkt, na.action=na.exclude)
   {residual.test[i] <- residuals(test[i])
   {CAPM.res[,i] <- residual.test[i]
   }}}

   but here I get the error message: "Error: $ operator is invalid for atomic
   vectors"
   and I don't think "test"  and "residual.test" is defined correctly because I
   don't know where to look for the residuals.

   Does anyone have an idea how to extract the residuals and put them in a
   258x50 matrix?
   Any help would be very much appreciated!

   Cheers,
   Ben
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] pie graphs in log scale axis

[Ortiz, John]

Dear list users,


I'm doing a plot integrating Grid output with Base Graphics output  (gridBase, 
Murrell 2012).

My goal is to produce a xy plot where each point is represented by a pie.  I 
could get it using the attach code,
but now I want to change the X axis to log scale.

when I introduce the log="x" parameter in the line:

plot(x, y, xlim=c(0.1, 1.2), log="x", ylim=c(0.1, 1.2), type="n")

I get this warning!

vps <- baseViewports()
Warning message:
In baseViewports() : viewport scales NOT set to user coordinates


This is the code to produce it without introduce log scale yet:

library(grid)
library(gridBase)
x <- c(0.88, 1.00, 0.67, 0.34)
y <- c(0.87, 0.43, 0.24, 0.94)
z <- matrix(runif(4*2), ncol=2)

oldpar <- par(no.readonly=TRUE)

plot(x, y, xlim=c(0.1, 1.2), ylim=c(0.1, 1.2), type="n")

vps <- baseViewports()
pushViewport(vps$inner, vps$figure, vps$plot)

for (i in 1:4) {
 pushViewport(viewport(x=unit(x[i], "native"),
 y=unit(y[i], "native"),
 width=0.1,
 height=0.1))
 par(plt=gridPLT(), new=TRUE)
 pie(z[i,], radius=1, labels=rep("", 2))
 popViewport()
 }

popViewport(2)
par(oldpar)

Thanks by some advice!!

Regards,


John Ortiz

Geologist
Smithsonian Tropical Research Institute Panama

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

Re: [R] Revo R for Arima Implementation

[Anindita Chattopadhyay]

Hi Rolf,

Thanks for the response. I have re-phrased the problem. Hope this will help.

We're working on a project where our model tries to predict the value of 
bookings as a time series of past bookings as well as some external flags.
The data was found to be stationary using the Dicky-Fuller test.
We have only a non-seasonal component for the AR and MA terms (small p & q 
only).
No differencing was done.

Then, we used values of p = 1, 7, 8 and q = 1, 7.
The main issue is that the p values are additive and q values are 
multiplicative.
We need to understand how we can implement this in Revo R.

Just to add, since we have external flags , we cannot make use of SARIMA 
function.
Please let me know if any other questions.

Thank you in advance!

Regards,
Anindita Chattopadhyay | +919886800606 | www.mu-sigma.com |

-Original Message-
From: Rolf Turner [mailto:r.tur...@auckland.ac.nz]
Sent: Saturday, October 26, 2013 1:56 AM
To: Anindita Chattopadhyay
Cc: r-help@r-project.org; Harish K
Subject: Re: [R] Revo R for Arima Implementation


Your question is pretty well totally opaque to me.  Describe the model you want 
to fit in mathematical terms, rather than referring to The Package That Must 
Not Be Named.  This is the ***R***  list.

It is possible that you might want to make use of the "seasonal"
argument to the
arima() function.

 cheers,

 Rolf Turner

On 10/26/13 01:19, Anindita Chattopadhyay wrote:
> Hello There,
>
> We have used ARIMA(multiplicative MA and additive AR) model in SAS to come to 
> our results. Please let me know how could that be implemented in R.
>
> Like,  In SAS we can pass the variable  as AR= (1,7) and MA=(1)(7) which is 
> additive and multiplicative respectively.
>
> In R we have the option as : order = c(0, 0, 0) which is (p,d,q).Is there any 
> specific way/code such that we specify the additive  AR and multiplicative MA 
> simultaneously?
>
This email message may contain proprietary, private and confidential 
information. The information transmitted is intended only for the person(s) or 
entities to which it is addressed. Any review, retransmission, dissemination or 
other use of, or taking of any action in reliance upon, this information by 
persons or entities other than the intended recipient is prohibited and may be 
illegal. If you received this in error, please contact the sender and delete 
the message from your system. Mu Sigma takes all reasonable steps to ensure 
that its electronic communications are free from viruses. However, given 
Internet accessibility, the Company cannot accept liability for any virus 
introduced by this e-mail or any attachment and you are advised to use 
up-to-date virus checking software.

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

Re: [R] bug in dummy.coef?

[Martin Maechler]

> Rolf Turner 
> on Sat, 26 Oct 2013 09:34:14 +1300 writes:

> On 10/25/13 19:24, Alexandra Kuznetsova wrote:
>> Thank you for the reply!  I need to use it in the package
>> for constructing the contrast matrices.
>  Really?
>> Renaming a variable like y2 <- y^2 will not be a solution
>> for me...

>  With sufficient effort and diligence (and programming
> skill) I am sure that this technique could be made to work
> in your context.
>> Will try to look whether the function dummy.coef can be
>> modified to work this issue.

>  As I said before, the problem is really buried inside
> the all.vars() function, whose workings are buried in the
> R internals.  This would not be at all easy (and would
> probably be dangerous) to modify.

Well, I've now spent much too much time on this now, trying to
explore possible changes.
It's definitely an order more complicated than you think, Rolf.

A change to all.vars() would break other cases and other code.
But I now strongly believe it's not even that part of
dummy.coef.lm() that needs to be amended.

Consider  example(dummy.coef) 's examples,
and also the test case from R/tests/reg-tests-1a.R, i.e., the
following block

  ## PR 1048 bug in dummy.coef.lm, Adrian Baddeley, 2001-08-10
  ## modified to give a sensible test
  old <- getOption("contrasts")
  options(contrasts=c("contr.helmert", "contr.poly"))
  DF <- data.frame(x=1:20,y=rnorm(20),z=factor(1:20 <= 10))
  dummy.coef(lm(y ~ z * I(x), data=DF))
  dummy.coef(lm(y ~ z * poly(x,1), data=DF))
  ## failed in 1.3.0.  Second one warns: deficiency of the method.
  options(contrasts=old)

where interestingly, there is no warning anymore, even though
the comment says there would be one, and I think I see the warning
in dummy.coef.lm's code ... which would been triggered somewhere in the
past ...

Many years ago, partly in S / S+ times, we used to advocate the
use of model.tables() over dummy.coef(), but as you can see on
its help page, it is neither implemented for all cases. 

Martin Maechler, ETH Zurich

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

Re: [R] Problem with coordinates when trying to draw lines into a raster (image) file

[Adams, Jean]

Rick,

I see a couple errors in my code ... you will need to attach the MASS
package for the eqscplot() function and that function's first arguments
should be any coordinates, I use (1, 1)

library(jpeg)
library(MASS)
# download image from internet for use in example
# http://www.tootsie.com/wallpaper/wp_concord2_1344x1008.jpg
img <- readJPEG("c:/temp/wp_concord2_1344x1008.jpg")
# define location of two yellow gumballs, could use locator() to get these
coordinates
px.x <- c(172, 576)
px.y <- c(433, 657)
# set up plot
par(xaxs="i", yaxs="i")
eqscplot(1, 1, type="n", xlim=c(0, 1334), ylim=c(0, 1008))
# add image
rasterImage(img, 0, 0, 1334, 1008)
# add line connecting two yellow gumballs
lines(px.x, px.y, lwd=5, col="yellow")

Jean



On Mon, Oct 28, 2013 at 4:12 AM, Rick Turner wrote:

> I forgot to mention that I get an error from your code - the 'lines()'
> command fails with the error 'object px.x not found'.  I am not sure why
> (yet) but need to dig into it... sigh.
>
> Rick
>
>

[[alternative HTML version deleted]]

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

Re: [R] plotting multiple variables

[Duncan Mackay]

Hi Hansol,

If you know the lattice system 

library(lattice)
library(latticeExtra)

useOuterStrips(
xyplot(y ~ x|class, data = dat, ...)
)

See also ? combineLimits

Otherwise

mfrow(10,10)

and loop your data  to plot

This is about as far as I can go - I have the same comment as Jim

Regards

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mac...@northnet.com.au

   

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Jim Lemon
Sent: Monday, 28 October 2013 16:25
To: Hansol Yu
Cc: r-help@r-project.org
Subject: Re: [R] plotting multiple variables

On 10/28/2013 04:53 PM, Hansol Yu wrote:
> When my data has 50 rows, 100 columns and class column. How can I plot 
> this data and show classes? Do I have to draw 100 lines?
>
Hi Hansol,
I doubt that anyone on the R help list really knows the answer to your
question unless your lecturer is one of us.

Jim

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

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

[R] Package for PAV and MPAV algorithms

[Babak Bastan]

Hi experts,

Is there any packag for PAV and MPAV algorithms (Pattern Anomaly Value)

[[alternative HTML version deleted]]

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

[R] coxph: how to define interaction terms?

[rm]

I’m trying to set up Cox Proptional Hazard model with interactions between
time and the covariates (which are categorical). The problem that I face is
that how to define the interactions, i.e. “x+cutStart:x”, properly.

The code below illustrates the problem. R gives the error message ” X matrix
deemed to be singular”, because “x+cutStart:x” includes too many
combinations of the dummies in the model. 

Any help is much appreciated!

aml2=survSplit(aml,cut=c(10,20,30),end='time',event='status',start='start')
aml2$cutStart=as.factor(aml2$start)
coxph(Surv(start,time,status)~x+cutStart:x,data=aml2)





--
View this message in context: 
http://r.789695.n4.nabble.com/coxph-how-to-define-interaction-terms-tp4679162.html
Sent from the R help mailing list archive at Nabble.com.

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

Re: [R] Code Book from SPSS Data

[Gesmann, Markus]

Peter,

There is a function called codebook as part of the memisc package that probably 
does what you want.
See also the package vignette: 
http://cran.r-project.org/web/packages/memisc/vignettes/anes48.pdf

library(memisc)
fn <- "NILT12w2.por"
dat <- spss.portable.file(fn)
codebook(dat)
names(dat)
description(dat)

I hope this helps.

Markus



-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of Peter Maclean
Sent: 26 October 2013 21:37
To: r-help@r-project.org
Subject: Re: [R] Code Book from SPSS Data

I do not have SPSS and I would like to create a code book in a data frame 
format using R. I am reading the SPSS file using "memisc package". The script 
is:
#Data for 2012 available at 
https://urldefense.proofpoint.com/v1/url?u=http://www.ark.ac.uk/nilt/datasets/&k=VTIXiGvdT7U4yPSpeHcrHQ%3D%3D%0A&r=dUkLGPeM%2BYkyyiRRq50yGs%2BmEf8kG%2FyCNQPwZn%2FaQD0%3D%0A&m=GVKjmjUiiJIiVcY%2FZ%2B35UGhRSRtdr9OJQhLQUx7Ax7Q%3D%0A&s=4dfa4233895c7e2f635380b0f8392a989a28848049001908965f058667a15143
#Also attached
ibrary(memisc)
## change the working directory
getwd()
setwd('')
data <- spss.portable.file("NILT12w2.por")
Get names
names(data) 
#Get Variable Lebels
des <- as.data.frame(description(data))
#Descriptive Statistics & Code Book
#Results are very long for printing
codebook(data)
#How could I extract a codebook (without Summary statistics for printing)?



Peter Maclean
Department of Economics
UDSM

--
The information in this E-Mail and in any attachments is...{{dropped:19}}

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

Re: [R] add a color band

[Alaios]

Hi Jim and thanks for your answer... I might be too tired with my new born or 
just exhausted.

I am sharing for everyone a small data snipset that you can load
https://www.dropbox.com/s/fh8jhwujgunmtrb/DataToPlotAsImage.Rdata



load("DataToPlotAsImage.Rdata")
require(plotrix)
browser()
test<-data
# this transforms the values of "test" into red->yellow
color2D.matplot(test,axes="F",xlab="",ylab="",main="color.scale",
  extremes=c("#FF","#00"),show.legend=FALSE)
  
axis(1,at=seq(1,ncol(test),length.out=10),labels=seq(201,300,length.out=10))
color.legend(104,30,112,70,seq(-110,-30,length=11),
  align="rb",rect.col=color.scale(1:30,1,c(0,1),0),gradient="y")

as you can see I have problems where the legend appears. My par("usr"  returned 
me
par("usr")
# [1]   0 351   0 200

but I am not sure how to read that to place the legend at a useful place. 
second I am not sure why the image is so full with black rows..

What I want is to have the legend visible
and later on customize the x axis to write custom string of different 
size... First I need though to fix the more severe problems as I have 
described

RegardsAlex



On Sunday, October 27, 2013 12:25 PM, Jim Lemon  wrote:
 
On 10/27/2013 08:39 AM, Alaios wrote:
> Hi Jim and thanks for your answer... I might be too tired with my new
> born or just exhausted.
>
> I am attaching for everyone a small data snipset that you can load
>
>
> load("DataToPlotAsImage.Rdata")
> require(plotrix)
> browser()
> test<-data
> # this transforms the values of "test" into red->yellow
> color2D.matplot(test,axes="F",xlab="",ylab="",main="color.scale",
> extremes=c("#FF","#00"),show.legend=FALSE)
>
> axis(1,at=seq(1,ncol(test),length.out=10),labels=seq(201,300,length.out=10))
> color.legend(104,30,112,70,seq(-110,-30,length=11),
> align="rb",rect.col=color.scale(1:30,1,c(0,1),0),gradient="y")
>
> as you can see I have problems where the legend appears. My par("usr"
> returned me
> par("usr")
> # [1] 0 351 0 200
>
> but I am not sure how to read that to place the legend at a useful place.
> second I am not sure why the image is so full with black rows..
>
> What I want is to have the legend visible
> and later on customize the x axis to write custom string of different
> size... First I need though to fix the more severe problems as I have
> described
>
Hi Alex,
I'm not sure why you have created a copy of "data" to plot it. I can get 
quite a sensible plot using this:

par(mar=c(5,4,4,5))
color2D.matplot(data,1,c(0,1),0,xlab="",ylab="",
  main="color.scale",xrange=c(-110,-50),border=NA)
color.legend(357,30,370,100,seq(-110,-50,length.out=6),
  align="rb",rect.col=color.scale(1:6,1,c(0,1),0),
  gradient="y")

Notice several things. First, when you have a large number of cells in a 
plot like this, setting the border to NA means that you don't get mostly 
borders (default = black) in the plot. The second thing is that your 
data range is -107.18150 to -54.07662. In order to get rounded numbers 
in your legend, I have set the xrange argument to -110 to -50. This 
gives a neat looking legend that spans your data, a bit like the
"pretty" function would do. It also means that the color mapping is to 
that range and is the same in the legend as in the plot. I have left 
enough space on the right of the plot to fit in the legend, as that was 
where you said you wanted it in your last email. What par("usr") tells 
you is the dimensions of the plot in user units. Here it is x=0 at the 
left, x=351 at the right, y=0 at the bottom and y=200 at the top.


Jim
[[alternative HTML version deleted]]

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