[R] linear model solving

[Ragia Ibrahim]

Dear group
IF I had an objective function and some constrains formed in linear model form. 
is there a way,..library in R that helps me to solve such amodel and find the 
unknown variable in it?

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

Re: [R] Ranking

[Bert Gunter]

It is perhaps worth mentioning that the OP's desire to do the
conversion to numeric to calculate won-lost percentages is completely
unnecessary and indicates that he/she could benefit by spending some
additional time learning R. See, e.g. ?tapply, ?table, ?prop.table,
and friends.

Cheers,
Bert


Bert Gunter

"Data is not information. Information is not knowledge. And knowledge
is certainly not wisdom."
   -- Clifford Stoll


On Sun, Nov 15, 2015 at 8:28 PM, David L Carlson  wrote:
> I used your code but deleted sep="\t" since there were no tabs in your email 
> and added the fill= argument I mentioned before.
>
>
>
> David
>
>
>
>  Original message 
> From: Ashta 
> Date: 11/14/2015 6:40 PM (GMT-06:00)
> To: David L Carlson 
> Cc: R help 
> Subject: Re: [R] Ranking
>
> Thank you David,
>
> My intention was if I change the status column  to numeric
> 0= Lost and 1 Won, then I can use this numeric variables  to calculate
> the  Percent game Won by each country.
> how did you read the data first?
> That was my problem.   The actual data is in a file have to be read or laded.
>
> Thank you !
>
>
>
>
>
>
> On Sat, Nov 14, 2015 at 6:10 PM, David L Carlson  wrote:
>> It is always good to read the manual page for a function, but especially 
>> when it is not working as you expected. In this case if you look at the 
>> arguments for read.table(), you will find one called fill=TRUE that is 
>> useful in this case.
>>
>> Based on your ifelse(), you seem to be assuming that a blank is not missing 
>> data but a lost game. You may also discover that in your example wins are 
>> coded as w and W.  Since character variables get converted to factors by 
>> default, you could use something like:
>>
>>> levels(test$STATUS) <- c("L", "W", "W")
>>> addmargins(xtabs(~Country+STATUS, test), 2)
>>STATUS
>> Country L W Sum
>> FRA 2 3   5
>> GER 1 3   4
>> SPA 2 1   3
>> UNK 1 2   3
>> USA 1 2   3
>>
>> I'll let you figure out how to get the last column.
>>
>> David L. Carlson
>> Department of Anthropology
>> Texas A&M University
>>
>> -Original Message-
>> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Ashta
>> Sent: Saturday, November 14, 2015 4:28 PM
>> To: R help 
>> Subject: [R] Ranking
>>
>> Hi all,
>>
>> I have the following raw data some records  don't have the second variable.
>>
>> test <- read.table(textConnection(" Country  STATUS
>> USA
>> USAW
>> USAW
>> GER
>> GERW
>> GERw
>> GERW
>> UNKW
>> UNK
>> UNKW
>> FRA
>> FRA
>> FRAW
>> FRAW
>> FRAW
>> SPA
>> SPAW
>> SPA  "),header = TRUE,  sep= "\t")
>> test
>>
>> It is not reading it correctly.
>>
>> Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings,  :
>>   line 17 did not have 2 elements
>>
>>
>>
>> After reading   I want change the status column  to numeric so that I
>> can use the table function
>>
>> test$STATUS <- ifelse(is.na(test$STATUS), 0,  1)
>>
>> at the end I want the following table (Country, Won, Lost , Number of
>> games played and % of score ) and pick the top 3 countries.
>>
>> COUNTRY   Won   Lost   NG%W
>>  USA 21 3  (2/3)*100
>>  GER 31 4  (3/4)*100
>>  UNK 21 3  (2/3)*100
>>  FRA 3 25  (3/5)*100
>>  SPA 1 2 3  (1/3)*100
>>
>> Thank you in  advance
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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

Re: [R] Ranking

[David L Carlson]

I used your code but deleted sep="\t" since there were no tabs in your email 
and added the fill= argument I mentioned before.



David



 Original message 
From: Ashta 
Date: 11/14/2015 6:40 PM (GMT-06:00)
To: David L Carlson 
Cc: R help 
Subject: Re: [R] Ranking

Thank you David,

My intention was if I change the status column  to numeric
0= Lost and 1 Won, then I can use this numeric variables  to calculate
the  Percent game Won by each country.
how did you read the data first?
That was my problem.   The actual data is in a file have to be read or laded.

Thank you !






On Sat, Nov 14, 2015 at 6:10 PM, David L Carlson  wrote:
> It is always good to read the manual page for a function, but especially when 
> it is not working as you expected. In this case if you look at the arguments 
> for read.table(), you will find one called fill=TRUE that is useful in this 
> case.
>
> Based on your ifelse(), you seem to be assuming that a blank is not missing 
> data but a lost game. You may also discover that in your example wins are 
> coded as w and W.  Since character variables get converted to factors by 
> default, you could use something like:
>
>> levels(test$STATUS) <- c("L", "W", "W")
>> addmargins(xtabs(~Country+STATUS, test), 2)
>STATUS
> Country L W Sum
> FRA 2 3   5
> GER 1 3   4
> SPA 2 1   3
> UNK 1 2   3
> USA 1 2   3
>
> I'll let you figure out how to get the last column.
>
> David L. Carlson
> Department of Anthropology
> Texas A&M University
>
> -Original Message-
> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Ashta
> Sent: Saturday, November 14, 2015 4:28 PM
> To: R help 
> Subject: [R] Ranking
>
> Hi all,
>
> I have the following raw data some records  don't have the second variable.
>
> test <- read.table(textConnection(" Country  STATUS
> USA
> USAW
> USAW
> GER
> GERW
> GERw
> GERW
> UNKW
> UNK
> UNKW
> FRA
> FRA
> FRAW
> FRAW
> FRAW
> SPA
> SPAW
> SPA  "),header = TRUE,  sep= "\t")
> test
>
> It is not reading it correctly.
>
> Error in scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings,  :
>   line 17 did not have 2 elements
>
>
>
> After reading   I want change the status column  to numeric so that I
> can use the table function
>
> test$STATUS <- ifelse(is.na(test$STATUS), 0,  1)
>
> at the end I want the following table (Country, Won, Lost , Number of
> games played and % of score ) and pick the top 3 countries.
>
> COUNTRY   Won   Lost   NG%W
>  USA 21 3  (2/3)*100
>  GER 31 4  (3/4)*100
>  UNK 21 3  (2/3)*100
>  FRA 3 25  (3/5)*100
>  SPA 1 2 3  (1/3)*100
>
> Thank you in  advance
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

[[alternative HTML version deleted]]

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

Re: [R] Why does a custom function called is.numeric.factor break lattice?

[Jeff Newmiller]

You need to read about S3 classes, and either make your custom function behave 
the way that function needs to behave or use a different function name for your 
custom function.

I think this is an example of the old saying that if it hurts when you slam 
your head against the wall, then don't do that.
---
Jeff NewmillerThe .   .  Go Live...
DCN:Basics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

On November 15, 2015 6:22:42 PM PST, sbihorel 
 wrote:
>Hi,
>
>Pretty much everything is in the title of the post. An example is
>below.
>
>library(lattice)
>data <- 
>data.frame(x=rep(1:10,8),y=rnorm(80),trt=factor(rep(1:4,each=20)),groups=rep(1:8,each=10))
>xyplot <- xyplot(y~x|trt,data,groups=groups)
>
>is.numeric.factor <- function(){
>   print('hello world')
>}
>
>xyplot <- xyplot(y~x|trt,data,groups=groups)
>
>Thanks for shedding some light on this.
>
>__
>R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide
>http://www.R-project.org/posting-guide.html
>and provide commented, minimal, self-contained, reproducible code.

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

Re: [R] Why does a custom function called is.numeric.factor break lattice?

[Bert Gunter]

Think about it.

I shall assume that you are familiar with S3 methods.  What do you
think would happen when xyplot code calls is.numeric() on a factor
object expecting it to call the is.numeric primitive but, instead,
finding a factor method defined, calls that? Note that your factor
method has no arguments, but the is.numeric() primitive does. Hence
when the code calls the primitive on the factor object, the error you
saw is thrown.

I would say that this is a weakness of the informal S3 "class" system,
although you probably should not have been surprised that is.numeric
is called on factors as the "x" argument, so you were inviting trouble
by defining a factor method that overrides this behavior.
Nevertheless, I would argue that one cannot know in general when this
occurs for other S3 classes, and that therefore allowing methods for
is.numeric() is dangerous.

Of course, full qualification in the original xyplot code
(base::is.numeric() rather than is.numeric() ) would avoid such
things, but that's a drag.

Contrary opinions and corrections to any flawed understanding on my
part are welcome, of course.

Cheers,
Bert


Bert Gunter

"Data is not information. Information is not knowledge. And knowledge
is certainly not wisdom."
   -- Clifford Stoll


On Sun, Nov 15, 2015 at 6:22 PM, sbihorel
 wrote:
> Hi,
>
> Pretty much everything is in the title of the post. An example is below.
>
> library(lattice)
> data <-
> data.frame(x=rep(1:10,8),y=rnorm(80),trt=factor(rep(1:4,each=20)),groups=rep(1:8,each=10))
> xyplot <- xyplot(y~x|trt,data,groups=groups)
>
> is.numeric.factor <- function(){
>   print('hello world')
> }
>
> xyplot <- xyplot(y~x|trt,data,groups=groups)
>
> Thanks for shedding some light on this.
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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

[R] Why does a custom function called is.numeric.factor break lattice?

[sbihorel]

Hi,

Pretty much everything is in the title of the post. An example is below.

library(lattice)
data <- 
data.frame(x=rep(1:10,8),y=rnorm(80),trt=factor(rep(1:4,each=20)),groups=rep(1:8,each=10))

xyplot <- xyplot(y~x|trt,data,groups=groups)

is.numeric.factor <- function(){
  print('hello world')
}

xyplot <- xyplot(y~x|trt,data,groups=groups)

Thanks for shedding some light on this.

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

Re: [R] Cautioning optim() users about "Nelder-Mead" default - (originally) Optim instability

[John C Frain]

In econometrics it was common to start an optimization with Nelder-Mead and
then switch to one of the other algorithms to finish the optimization. As
John Nash states NM gets one close. switching then speeds the final
solution.

John

John C Frain
3 Aranleigh Park
Rathfarnham
Dublin 14
Ireland
www.tcd.ie/Economics/staff/frainj/home.html
mailto:fra...@tcd.ie
mailto:fra...@gmail.com

On 15 November 2015 at 20:05, Mark Leeds  wrote:

> and just to add to john's comments, since he's too modest, in my
> experience,  the algorithm in the rvmmin  package ( written by john ) shows
> great improvement compared to the L-BFGS-B  algorithm so I don't use
> L-BFGS-B anymore.  L-BFGS-B often has a dangerous convergence issue  in
> that it can claim to converge when it hasn't. which, to
> me is worse than not converging.  Most likely it has to do with the link
> below which was pointed out to me by john a while back.
>
> http://www.ece.northwestern.edu/~morales/PSfiles/acm-remark.pdf
>
>
> On Sun, Nov 15, 2015 at 2:41 PM, ProfJCNash  wrote:
>
> > Agreed on the default algorithm issue. That is important for users to
> > know, and I'm happy to underline it. Also that CG (which is based on one
> > of my codes) should be deprecated. BFGS (also based on one of my codes
> > from long ago) does much better than I would ever have expected.
> >
> > Over the years I've tried different Nelder-Mead implementations. Cannot
> > say I've found any that is always better than that in optim() (also
> > based on an old code of mine), though nmkb() from dfoptim package seems
> > to do better a lot of the time, and it has a transformation method for
> > bounds, which may be useful, but does have the awkwardness that one
> > cannot start on a bound. For testing a function, I don't think it makes
> > a lot of difference which variant of NM one uses if the trace is on to
> > catch never-ending runs. For production use, it is a really good idea to
> > try different methods on a sample of likely cases and choose a method
> > that does well. That is the motivation for the optimx package or the
> > opm() function of the newer optimz (on R-forge) that I'm still
> > polishing. optimz has a function optimr() that has the same call as
> > optim() but incorporates over a dozen optimizers via method = "(selected
> > method)".
> >
> > As a gradient-free choice, the Powell codes from minqa or other packages
> > (there are several implementations) can sometimes have spectacular
> > performance, but they also flub rather more regularly than Nelder-Mead
> > in my experience. That is, when they are good, they are very very good,
> > and when they are not they are horrid. (Plagiarism!)
> >
> > JN
> >
> > On 15-11-15 12:46 PM, Ravi Varadhan wrote:
> > > Hi John,
> > > My main point is not about Nelder-Mead per se.  It is *primarily* about
> > the Nelder-Mead implementation in optim().
> > >
> > > The users of optim() should be cautioned regarding the default
> algorithm
> > and that they should consider alternatives such as "BFGS" in optim(), or
> > other implementations of Nelder-Mead.
> > >
> > > Best regards,
> > > Ravi
> > > 
> > > From: ProfJCNash 
> > > Sent: Sunday, November 15, 2015 12:21 PM
> > > To: Ravi Varadhan; 'r-help@r-project.org'; lorenzo.ise...@gmail.com
> > > Cc: b...@xs4all.nl; Gabor Grothendieck
> > > Subject: Re: Cautioning optim() users about "Nelder-Mead" default -
> > (originally) Optim instability
> > >
> > > Not contradicting Ravi's message, but I wouldn't say Nelder-Mead is
> > > "bad" per se. It's issues are that it assumes the parameters are all on
> > > the same scale, and the termination (not convergence) test can't use
> > > gradients, so it tends to get "near" the optimum very quickly -- say
> > > only 10% of the computational effort -- then spends an awful amount of
> > > effort deciding it's got there. It often will do poorly when the
> > > function has nearly "flat" zones e.g., long valley with very low slope.
> > >
> > > So my message is still that Nelder-Mead is an unfortunate default -- it
> > > has been chosen I believe because it is generally robust and doesn't
> > > need gradients. BFGS really should use accurate gradients, preferably
> > > computed analytically, so it would only be a good default in that case
> > > or with very good approximate gradients (which are costly
> > > computationally).
> > >
> > > However, if you understand what NM is doing, and use it accordingly, it
> > > is a valuable tool. I generally use it as a first try BUT turn on the
> > > trace to watch what it is doing as a way to learn a bit about the
> > > function I am minimizing. Rarely would I use it as a production
> > minimizer.
> > >
> > > Best, JN
> > >
> > > On 15-11-15 12:02 PM, Ravi Varadhan wrote:
> > >> Hi,
> > >>
> > >>
> > >>
> > >> While I agree with the comments about paying attention to parameter
> > >> scaling, a major issue here is that the default optimization
> algorithm,
> > >> Nelder-

[R] OT, apropos of nothing ....

[Rolf Turner]

In respect of Bert Gunter's signature quote:

"Data is not information. Information is not knowledge. And knowledge
is certainly not wisdom."

The other day my wife saw a grocery truck with the following remark 
emblazoned on its side:


"Knowledge is being aware that a tomato is a fruit.  Wisdom is 
refraining from putting tomatoes in a fruit salad."


cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

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

Re: [R] Cautioning optim() users about "Nelder-Mead" default - (originally) Optim instability

[Mark Leeds]

and just to add to john's comments, since he's too modest, in my
experience,  the algorithm in the rvmmin  package ( written by john ) shows
great improvement compared to the L-BFGS-B  algorithm so I don't use
L-BFGS-B anymore.  L-BFGS-B often has a dangerous convergence issue  in
that it can claim to converge when it hasn't. which, to
me is worse than not converging.  Most likely it has to do with the link
below which was pointed out to me by john a while back.

http://www.ece.northwestern.edu/~morales/PSfiles/acm-remark.pdf


On Sun, Nov 15, 2015 at 2:41 PM, ProfJCNash  wrote:

> Agreed on the default algorithm issue. That is important for users to
> know, and I'm happy to underline it. Also that CG (which is based on one
> of my codes) should be deprecated. BFGS (also based on one of my codes
> from long ago) does much better than I would ever have expected.
>
> Over the years I've tried different Nelder-Mead implementations. Cannot
> say I've found any that is always better than that in optim() (also
> based on an old code of mine), though nmkb() from dfoptim package seems
> to do better a lot of the time, and it has a transformation method for
> bounds, which may be useful, but does have the awkwardness that one
> cannot start on a bound. For testing a function, I don't think it makes
> a lot of difference which variant of NM one uses if the trace is on to
> catch never-ending runs. For production use, it is a really good idea to
> try different methods on a sample of likely cases and choose a method
> that does well. That is the motivation for the optimx package or the
> opm() function of the newer optimz (on R-forge) that I'm still
> polishing. optimz has a function optimr() that has the same call as
> optim() but incorporates over a dozen optimizers via method = "(selected
> method)".
>
> As a gradient-free choice, the Powell codes from minqa or other packages
> (there are several implementations) can sometimes have spectacular
> performance, but they also flub rather more regularly than Nelder-Mead
> in my experience. That is, when they are good, they are very very good,
> and when they are not they are horrid. (Plagiarism!)
>
> JN
>
> On 15-11-15 12:46 PM, Ravi Varadhan wrote:
> > Hi John,
> > My main point is not about Nelder-Mead per se.  It is *primarily* about
> the Nelder-Mead implementation in optim().
> >
> > The users of optim() should be cautioned regarding the default algorithm
> and that they should consider alternatives such as "BFGS" in optim(), or
> other implementations of Nelder-Mead.
> >
> > Best regards,
> > Ravi
> > 
> > From: ProfJCNash 
> > Sent: Sunday, November 15, 2015 12:21 PM
> > To: Ravi Varadhan; 'r-help@r-project.org'; lorenzo.ise...@gmail.com
> > Cc: b...@xs4all.nl; Gabor Grothendieck
> > Subject: Re: Cautioning optim() users about "Nelder-Mead" default -
> (originally) Optim instability
> >
> > Not contradicting Ravi's message, but I wouldn't say Nelder-Mead is
> > "bad" per se. It's issues are that it assumes the parameters are all on
> > the same scale, and the termination (not convergence) test can't use
> > gradients, so it tends to get "near" the optimum very quickly -- say
> > only 10% of the computational effort -- then spends an awful amount of
> > effort deciding it's got there. It often will do poorly when the
> > function has nearly "flat" zones e.g., long valley with very low slope.
> >
> > So my message is still that Nelder-Mead is an unfortunate default -- it
> > has been chosen I believe because it is generally robust and doesn't
> > need gradients. BFGS really should use accurate gradients, preferably
> > computed analytically, so it would only be a good default in that case
> > or with very good approximate gradients (which are costly
> > computationally).
> >
> > However, if you understand what NM is doing, and use it accordingly, it
> > is a valuable tool. I generally use it as a first try BUT turn on the
> > trace to watch what it is doing as a way to learn a bit about the
> > function I am minimizing. Rarely would I use it as a production
> minimizer.
> >
> > Best, JN
> >
> > On 15-11-15 12:02 PM, Ravi Varadhan wrote:
> >> Hi,
> >>
> >>
> >>
> >> While I agree with the comments about paying attention to parameter
> >> scaling, a major issue here is that the default optimization algorithm,
> >> Nelder-Mead, is not very good.  It is unfortunate that the optim
> >> implementation chose this as the "default" algorithm.  I have several
> >> instances where people have come to me with poor results from using
> >> optim(), because they did not realize that the default algorithm is
> >> bad.  We (John Nash and I) have pointed this out before, but the R core
> >> has not addressed this issue due to backward compatibility reasons.
> >>
> >>
> >>
> >> There is a better implementation of Nelder-Mead in the "dfoptim"
> package.
> >>
> >>
> >>
> >> ​require(dfoptim)
> >>
> >> mm_def1 <- nmk(par = par_ini1, min.perc_erro

Re: [R] Cautioning optim() users about "Nelder-Mead" default - (originally) Optim instability

[ProfJCNash]

Agreed on the default algorithm issue. That is important for users to
know, and I'm happy to underline it. Also that CG (which is based on one
of my codes) should be deprecated. BFGS (also based on one of my codes
from long ago) does much better than I would ever have expected.

Over the years I've tried different Nelder-Mead implementations. Cannot
say I've found any that is always better than that in optim() (also
based on an old code of mine), though nmkb() from dfoptim package seems
to do better a lot of the time, and it has a transformation method for
bounds, which may be useful, but does have the awkwardness that one
cannot start on a bound. For testing a function, I don't think it makes
a lot of difference which variant of NM one uses if the trace is on to
catch never-ending runs. For production use, it is a really good idea to
try different methods on a sample of likely cases and choose a method
that does well. That is the motivation for the optimx package or the
opm() function of the newer optimz (on R-forge) that I'm still
polishing. optimz has a function optimr() that has the same call as
optim() but incorporates over a dozen optimizers via method = "(selected
method)".

As a gradient-free choice, the Powell codes from minqa or other packages
(there are several implementations) can sometimes have spectacular
performance, but they also flub rather more regularly than Nelder-Mead
in my experience. That is, when they are good, they are very very good,
and when they are not they are horrid. (Plagiarism!)

JN

On 15-11-15 12:46 PM, Ravi Varadhan wrote:
> Hi John,
> My main point is not about Nelder-Mead per se.  It is *primarily* about the 
> Nelder-Mead implementation in optim().  
> 
> The users of optim() should be cautioned regarding the default algorithm and 
> that they should consider alternatives such as "BFGS" in optim(), or other 
> implementations of Nelder-Mead.
> 
> Best regards,
> Ravi
> 
> From: ProfJCNash 
> Sent: Sunday, November 15, 2015 12:21 PM
> To: Ravi Varadhan; 'r-help@r-project.org'; lorenzo.ise...@gmail.com
> Cc: b...@xs4all.nl; Gabor Grothendieck
> Subject: Re: Cautioning optim() users about "Nelder-Mead" default - 
> (originally) Optim instability
> 
> Not contradicting Ravi's message, but I wouldn't say Nelder-Mead is
> "bad" per se. It's issues are that it assumes the parameters are all on
> the same scale, and the termination (not convergence) test can't use
> gradients, so it tends to get "near" the optimum very quickly -- say
> only 10% of the computational effort -- then spends an awful amount of
> effort deciding it's got there. It often will do poorly when the
> function has nearly "flat" zones e.g., long valley with very low slope.
> 
> So my message is still that Nelder-Mead is an unfortunate default -- it
> has been chosen I believe because it is generally robust and doesn't
> need gradients. BFGS really should use accurate gradients, preferably
> computed analytically, so it would only be a good default in that case
> or with very good approximate gradients (which are costly
> computationally).
> 
> However, if you understand what NM is doing, and use it accordingly, it
> is a valuable tool. I generally use it as a first try BUT turn on the
> trace to watch what it is doing as a way to learn a bit about the
> function I am minimizing. Rarely would I use it as a production minimizer.
> 
> Best, JN
> 
> On 15-11-15 12:02 PM, Ravi Varadhan wrote:
>> Hi,
>>
>>
>>
>> While I agree with the comments about paying attention to parameter
>> scaling, a major issue here is that the default optimization algorithm,
>> Nelder-Mead, is not very good.  It is unfortunate that the optim
>> implementation chose this as the "default" algorithm.  I have several
>> instances where people have come to me with poor results from using
>> optim(), because they did not realize that the default algorithm is
>> bad.  We (John Nash and I) have pointed this out before, but the R core
>> has not addressed this issue due to backward compatibility reasons.
>>
>>
>>
>> There is a better implementation of Nelder-Mead in the "dfoptim" package.
>>
>>
>>
>> ​require(dfoptim)
>>
>> mm_def1 <- nmk(par = par_ini1, min.perc_error, data = data)
>>
>> mm_def2 <- nmk(par = par_ini2, min.perc_error, data = data)
>>
>> mm_def3 <- nmk(par = par_ini3, min.perc_error, data = data)
>>
>> print(mm_def1$par)
>>
>> print(mm_def2$par)
>>
>> print(mm_def3$par)
>>
>>
>>
>> In general, better implementations of optimization algorithms are
>> available in packages such as "optimx", "nloptr".  It is unfortunate
>> that most naïve users of optimization in R do not recognize this.
>> Perhaps, there should be a "message" in the optim help file that points
>> this out to the users.
>>
>>
>>
>> Hope this is helpful,
>>
>> Ravi
>>
>>
>

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-

Re: [R] Cautioning optim() users about "Nelder-Mead" default - (originally) Optim instability

[Ravi Varadhan]

Hi John,
My main point is not about Nelder-Mead per se.  It is *primarily* about the 
Nelder-Mead implementation in optim().  

The users of optim() should be cautioned regarding the default algorithm and 
that they should consider alternatives such as "BFGS" in optim(), or other 
implementations of Nelder-Mead.

Best regards,
Ravi

From: ProfJCNash 
Sent: Sunday, November 15, 2015 12:21 PM
To: Ravi Varadhan; 'r-help@r-project.org'; lorenzo.ise...@gmail.com
Cc: b...@xs4all.nl; Gabor Grothendieck
Subject: Re: Cautioning optim() users about "Nelder-Mead" default - 
(originally) Optim instability

Not contradicting Ravi's message, but I wouldn't say Nelder-Mead is
"bad" per se. It's issues are that it assumes the parameters are all on
the same scale, and the termination (not convergence) test can't use
gradients, so it tends to get "near" the optimum very quickly -- say
only 10% of the computational effort -- then spends an awful amount of
effort deciding it's got there. It often will do poorly when the
function has nearly "flat" zones e.g., long valley with very low slope.

So my message is still that Nelder-Mead is an unfortunate default -- it
has been chosen I believe because it is generally robust and doesn't
need gradients. BFGS really should use accurate gradients, preferably
computed analytically, so it would only be a good default in that case
or with very good approximate gradients (which are costly
computationally).

However, if you understand what NM is doing, and use it accordingly, it
is a valuable tool. I generally use it as a first try BUT turn on the
trace to watch what it is doing as a way to learn a bit about the
function I am minimizing. Rarely would I use it as a production minimizer.

Best, JN

On 15-11-15 12:02 PM, Ravi Varadhan wrote:
> Hi,
>
>
>
> While I agree with the comments about paying attention to parameter
> scaling, a major issue here is that the default optimization algorithm,
> Nelder-Mead, is not very good.  It is unfortunate that the optim
> implementation chose this as the "default" algorithm.  I have several
> instances where people have come to me with poor results from using
> optim(), because they did not realize that the default algorithm is
> bad.  We (John Nash and I) have pointed this out before, but the R core
> has not addressed this issue due to backward compatibility reasons.
>
>
>
> There is a better implementation of Nelder-Mead in the "dfoptim" package.
>
>
>
> ​require(dfoptim)
>
> mm_def1 <- nmk(par = par_ini1, min.perc_error, data = data)
>
> mm_def2 <- nmk(par = par_ini2, min.perc_error, data = data)
>
> mm_def3 <- nmk(par = par_ini3, min.perc_error, data = data)
>
> print(mm_def1$par)
>
> print(mm_def2$par)
>
> print(mm_def3$par)
>
>
>
> In general, better implementations of optimization algorithms are
> available in packages such as "optimx", "nloptr".  It is unfortunate
> that most naïve users of optimization in R do not recognize this.
> Perhaps, there should be a "message" in the optim help file that points
> this out to the users.
>
>
>
> Hope this is helpful,
>
> Ravi
>
>

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

[R] Cautioning optim() users about "Nelder-Mead" default - (originally) Optim instability

[Ravi Varadhan]

Hi,



While I agree with the comments about paying attention to parameter scaling, a 
major issue here is that the default optimization algorithm, Nelder-Mead, is 
not very good.  It is unfortunate that the optim implementation chose this as 
the "default" algorithm.  I have several instances where people have come to me 
with poor results from using optim(), because they did not realize that the 
default algorithm is bad.  We (John Nash and I) have pointed this out before, 
but the R core has not addressed this issue due to backward compatibility 
reasons.



There is a better implementation of Nelder-Mead in the "dfoptim" package.



?require(dfoptim)

mm_def1 <- nmk(par = par_ini1, min.perc_error, data = data)

mm_def2 <- nmk(par = par_ini2, min.perc_error, data = data)

mm_def3 <- nmk(par = par_ini3, min.perc_error, data = data)

print(mm_def1$par)

print(mm_def2$par)

print(mm_def3$par)



In general, better implementations of optimization algorithms are available in 
packages such as "optimx", "nloptr".  It is unfortunate that most na�ve users 
of optimization in R do not recognize this.  Perhaps, there should be a 
"message" in the optim help file that points this out to the users.



Hope this is helpful,

Ravi


[[alternative HTML version deleted]]

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

Re: [R] Two Time Fixed Effects - LFE package

[Miluji Sb]

Dear Andrew,

Thank you for your reply. Its an R question. The weeks are coded as 1-53
for each year and I would like to control weeks and years as time fixed
effects.

Will this be an issue if I estimate this type of regression using the LFE
package?

felm(outcome ~ temperature + precipitation | city + year + week

Thanks again!

Sincerely,

MS

On Sun, Nov 15, 2015 at 12:13 AM, Andrew Crane-Droesch 
wrote:

> Is this an R question or an econometrics question?  I'll assume that it is
> an R question.  If your weeks are coded sequentially (i.e.: weeks since a
> particular date), then they'll be strictly determined by year.  If however
> you're interested in the effect of a particular week of the year (week 7,
> for example), then you'll need to recode your week variable as a factor
> with 52 levels.  For that you'd likely need the "%%" operator.  For example:
>
> 1> 1:10%%3
>  [1] 1 2 0 1 2 0 1 2 0 1
>
>
>
>
> On 11/14/2015 05:18 PM, Miluji Sb wrote:
>
>> I have weekly panel data for more than a hundred cities. The independent
>> variables are temperature and precipitation. The time dimensions are year
>> and week and likely have time invariant characteristics and are all
>> important for proper estimation.
>>
>> Could I use the LFE (or plm) package to estimate something like this by
>> including the location and two time fixed-effects?
>>
>> felm(outcome ~ temperature + precipitation | city + year + week
>>
>> Thanks!
>>
>> MS
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>

[[alternative HTML version deleted]]

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

Re: [R] Cautioning optim() users about "Nelder-Mead" default - (originally) Optim instability

[ProfJCNash]

Not contradicting Ravi's message, but I wouldn't say Nelder-Mead is
"bad" per se. It's issues are that it assumes the parameters are all on
the same scale, and the termination (not convergence) test can't use
gradients, so it tends to get "near" the optimum very quickly -- say
only 10% of the computational effort -- then spends an awful amount of
effort deciding it's got there. It often will do poorly when the
function has nearly "flat" zones e.g., long valley with very low slope.

So my message is still that Nelder-Mead is an unfortunate default -- it
has been chosen I believe because it is generally robust and doesn't
need gradients. BFGS really should use accurate gradients, preferably
computed analytically, so it would only be a good default in that case
or with very good approximate gradients (which are costly
computationally).

However, if you understand what NM is doing, and use it accordingly, it
is a valuable tool. I generally use it as a first try BUT turn on the
trace to watch what it is doing as a way to learn a bit about the
function I am minimizing. Rarely would I use it as a production minimizer.

Best, JN

On 15-11-15 12:02 PM, Ravi Varadhan wrote:
> Hi,
> 
>  
> 
> While I agree with the comments about paying attention to parameter
> scaling, a major issue here is that the default optimization algorithm,
> Nelder-Mead, is not very good.  It is unfortunate that the optim
> implementation chose this as the "default" algorithm.  I have several
> instances where people have come to me with poor results from using
> optim(), because they did not realize that the default algorithm is
> bad.  We (John Nash and I) have pointed this out before, but the R core
> has not addressed this issue due to backward compatibility reasons. 
> 
>  
> 
> There is a better implementation of Nelder-Mead in the "dfoptim" package.
> 
>  
> 
> ​require(dfoptim)
> 
> mm_def1 <- nmk(par = par_ini1, min.perc_error, data = data)
> 
> mm_def2 <- nmk(par = par_ini2, min.perc_error, data = data)
> 
> mm_def3 <- nmk(par = par_ini3, min.perc_error, data = data)
> 
> print(mm_def1$par)
> 
> print(mm_def2$par)
> 
> print(mm_def3$par)
> 
>  
> 
> In general, better implementations of optimization algorithms are
> available in packages such as "optimx", "nloptr".  It is unfortunate
> that most naïve users of optimization in R do not recognize this. 
> Perhaps, there should be a "message" in the optim help file that points
> this out to the users. 
> 
>  
> 
> Hope this is helpful,
> 
> Ravi
> 
>

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

Re: [R] Alternatives for explicit for() loops

[marammagdysalem]

Thanks a lot Jim and Boris for replying.

Sent from my iPhone

> On Nov 9, 2015, at 1:13 AM, jim holtman  wrote:
> 
> You need to take a close look at the function incomb that you are creating.  
> I see what appears to be a constant value 
> ("*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))") being computed that you might 
> only have to compute once before the function.  You are also referencing many 
> variables (m, LED, j, ...) that are not being passed in that are in the 
> global environment; it is probably better to pass them in to the function.  I 
> am not sure what 'pbapply' is doing for you since I see this is new to the 
> code that you first sent out.
> 
> I would be good it you told us what the function is trying to do; you are 
> showing us how you want to do it, not what you want to do.  Are there other 
> ways of doing it?  If speed is your problem, then consider the "Rcpp" package 
> and write the function is C++ which might be faster, but again, take a look 
> at what you are doing to see if there are other ways.  I don't have time to 
> dig into the code, since there is a lack of comments, to understand why you 
> are using, e.g., 'choose', 'prod', etc.).  There are probably a lot of ways 
> of speeding up the code, if could tell us what you want to accomplish.
> 
> 
> Jim Holtman
> Data Munger Guru
>  
> What is the problem that you are trying to solve?
> Tell me what you want to do, not how you want to do it.
> 
>> On Sun, Nov 8, 2015 at 4:48 AM, Maram SAlem  
>> wrote:
>> Thanks all for replying.
>> 
>> In fact I've used the the Rprof() function and found out that the incomb() 
>> function (in my code above)  takes about 80% of the time, but I didn't 
>> figure out which part of the function is causing the delay. So I thought 
>> that this may be due to the for() loops. 
>> I MUST run this code for rather large values of n and m, so is there any way 
>> that can help me do that without having to wait for more than three days to 
>> reach an output. N.B. I'll have to repeat these runs for may be 70 or 80 
>> times , and this means HUGE time
>> 
>> I'd appreciate any sort of help.
>> Thanks in advance.
>> 
>> Maram Salem
>> 
>>> On 6 November 2015 at 16:54, jim holtman  wrote:
>>> If you have code that is running for a long time, then take a small case 
>>> that only runs for 5-10 minutes and turn on the RProfiler so that you can 
>>> see where you are spending your time.  In most cases, it is probably not 
>>> the 'for' loops that are causing the problem, but some function/calculation 
>>> you are doing within the loop that is consuming the time, and until you 
>>> determine what section of code that is, is it hard to tell exactly what the 
>>> problem is, much less the solution.
>>> 
>>> 
>>> Jim Holtman
>>> Data Munger Guru
>>>  
>>> What is the problem that you are trying to solve?
>>> Tell me what you want to do, not how you want to do it.
>>> 
 On Wed, Nov 4, 2015 at 9:09 AM, Maram SAlem  
 wrote:
 Hi Jim,
 
 Thanks a lot for replying.
 
 In fact I'm trying to run a simulation study that enables me to calculate 
 the Bayes risk of a sampling plan selected from progressively type-II 
 censored Weibull model. One of the steps involves evaluating the expected 
 test time, which is a rather complicated formula that involves nested 
 multiple summations where the counters of the summation signs are 
 dependent, that's why I thought of I should create the incomb() function 
 inside the loop, or may be I didn't figure out how to relate its arguments 
 to the ones inside the loop had I created it outside it.  I'm trying to 
 create a matrix of all the possible combinations involved in the 
 summations and then use the apply() function on each row of that matrix. 
 The problem is that the code I wrote works perfectly well for rather small 
 values of the sample size,n, and the censoring number, m (for example, 
 n=8,m=4),but when n and m are increased (say, n=25,m=15) the code keeps on 
 running for days with no output. That's why I thoug!
 ht I should try to avoid explicit loops as much as possible, so I did my best 
in this regard but still the code takes too long to execute,(more than three 
days), thus, i believe there must be something wrong.
 
 Here's the full code:
 
 library(pbapply)
 f1 <- function(n, m) {
stopifnot(n > m)
r0 <- t(diff(combn(n-1, m-1)) - 1L)
r1 <- rep(seq(from=0, len=n-m+1), choose( seq(to=m-2, by=-1, 
 len=n-m+1), m-2))
cbind(r0[, ncol(r0):1, drop=FALSE], r1, deparse.level=0)
 }
 simpfun<- function (x,n,m,p,alpha,beta)
   {
   a<-factorial(n-m)/(prod((factorial(x)))*(factorial((n-m)- sum(x
   b <-  ((m-1):1)
   c<- a*((p)^(sum(x)))*((1-p)^(((m-1)*(n-m))- sum(x%*%(as.matrix(b)
 d <- n - cumsum(x) - (1:(m-1))
   e<- n*(prod(d))*c
 LD<-list()
for (i in 1:(m-1))  {
LD[[i]]<-seq(

Re: [R] Cautioning optim() users about "Nelder-Mead" default - (originally) Optim instability

[lorenzo.ise...@gmail.com]

Thanks a lot, Ravi.
Indeed you best understood the point of my email.
I am perfectly aware that most of the optimization algorithms find
local rather than global minima and therefore the choice of the
initial parameters plays (at least in principle) a role.
Nevertheless, my optimization problem is rather trivial and I did not
bother to look for anything beyond the most basic tool in R for
optimization.
What surprised me is that an algorithm different from the default one
in optim() is extremely robust to a partially deliberate bad choice
ofthe initial parameters, whereas the standard one is not.
You perfectly answered my question.
Regards

Lorenzo


On Sun, Nov 15, 2015 at 05:02:52PM +, Ravi Varadhan wrote:

Hi,



While I agree with the comments about paying attention to parameter scaling, a major 
issue here is that the default optimization algorithm, Nelder-Mead, is not very good.  It 
is unfortunate that the optim implementation chose this as the "default" 
algorithm.  I have several instances where people have come to me with poor results from 
using optim(), because they did not realize that the default algorithm is bad.  We (John 
Nash and I) have pointed this out before, but the R core has not addressed this issue due 
to backward compatibility reasons.



There is a better implementation of Nelder-Mead in the "dfoptim" package.



?require(dfoptim)

mm_def1 <- nmk(par = par_ini1, min.perc_error, data = data)

mm_def2 <- nmk(par = par_ini2, min.perc_error, data = data)

mm_def3 <- nmk(par = par_ini3, min.perc_error, data = data)

print(mm_def1$par)

print(mm_def2$par)

print(mm_def3$par)



In general, better implementations of optimization algorithms are available in packages such as 
"optimx", "nloptr".  It is unfortunate that most naïve users of optimization in R do not 
recognize this.  Perhaps, there should be a "message" in the optim help file that points this out 
to the users.



Hope this is helpful,

Ravi



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

Re: [R] Two Time Fixed Effects - LFE package

[Andrew Crane-Droesch]

I don't see why not, but I also don't see why you need to take my word 
for it when you can compare the output of felm against the output of lm, 
with dummy variables for all the factors.  If that many dummies is 
computationally tough, just work with a subset.

On 11/15/2015 08:37 AM, Miluji Sb wrote:
> Dear Andrew,
>
> Thank you for your reply. Its an R question. The weeks are coded as 
> 1-53 for each year and I would like to control weeks and years as time 
> fixed effects.
>
> Will this be an issue if I estimate this type of regression using the 
> LFE package?
>
> felm(outcome ~ temperature + precipitation | city + year + week
>
> Thanks again!
>
> Sincerely,
>
> MS
>
> On Sun, Nov 15, 2015 at 12:13 AM, Andrew Crane-Droesch 
> mailto:andre...@gmail.com>> wrote:
>
> Is this an R question or an econometrics question?  I'll assume
> that it is an R question.  If your weeks are coded sequentially
> (i.e.: weeks since a particular date), then they'll be strictly
> determined by year.  If however you're interested in the effect of
> a particular week of the year (week 7, for example), then you'll
> need to recode your week variable as a factor with 52 levels.  For
> that you'd likely need the "%%" operator.  For example:
>
> 1> 1:10%%3
>  [1] 1 2 0 1 2 0 1 2 0 1
>
>
>
>
> On 11/14/2015 05:18 PM, Miluji Sb wrote:
>
> I have weekly panel data for more than a hundred cities. The
> independent
> variables are temperature and precipitation. The time
> dimensions are year
> and week and likely have time invariant characteristics and
> are all
> important for proper estimation.
>
> Could I use the LFE (or plm) package to estimate something
> like this by
> including the location and two time fixed-effects?
>
> felm(outcome ~ temperature + precipitation | city + year + week
>
> Thanks!
>
> MS
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org  mailing
> list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
>


[[alternative HTML version deleted]]

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