[R] repository for ubuntu/linux mint for R 4.0.0

[ProfJCNash]

In updating (an older computer with) Linux Mint 18.3 I tried to add
the repository

deb https://cloud.r-project.org/bin/linux/ubuntu xenial-cran40/

as per the "Download R for Linux" instructions. This gave an error
that there was no Release.key file.

After some investigation, I found that

deb https://cran.r-project.org/bin/linux/ubuntu xenial-cran40/

i.e., CRAN not CLOUD. With this change, I could install R 4.0.

Is this a known glitch?

JN

__
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] I can't get seq to behave how I think it should

[ProfJCNash]

As one of the approximately 30 names on the 1985 IEEE 754 standard, I
should be first to comment about representations. However, a quite
large fraction of the computers I've owned or used were decimal beasts.
This doesn't remove all the issues, of course, but some of these
input-output conversions would be avoided. I rather doubt we'll ever
see an R version for decimal arithmetic, since there'd be a lot of
awkwardness with comparing with the usual version, and a lot of the test
comparisons would likely fail. On the other hand, it might serve as a
reminder that IEEE arithmetic, while a great step in computation, is not
the only possibility and does not eliminate all the difficulties with
finite precision arithmetic.

Best, JN




On 2019-01-17 11:52 a.m., peter dalgaard wrote:
> 
> 
>> On 17 Jan 2019, at 15:56 , POLWART, Calum (COUNTY DURHAM AND DARLINGTON NHS 
>> FOUNDATION TRUST) via R-help  wrote:
>>
>> Well I get the issue with finite precision. As in SQRT(2) * SQRT(2) is not 2.
> 
> As Jeff indicates, you also need to get that just like 3rds and 7ths cannot 
> be represented exactly in base 10, 5ths and 10ths cannot be represented 
> exactly in base 2 (only powers of 1/2 and their multiples can).
> 
> Specifically, 1.4 decimal is 1.0110011001100 binary
> 
> -pd
> 
>>
>> What surprised me was that seq(1.4, 2.1, by=0.001) starts at 
>> 1.3999 and not 1.4!
>>
>>
>> -Original Message-
>> From: PIKAL Petr [mailto:petr.pi...@precheza.cz]
>> Sent: 17 January 2019 14:30
>> To: POLWART, Calum (COUNTY DURHAM AND DARLINGTON NHS FOUNDATION TRUST); Ben 
>> Tupper
>> Cc: r-help@r-project.org
>> Subject: RE: [R] I can't get seq to behave how I think it should
>>
>> Hi
>>
>> It is not seq problem, it is floating point numbers representation in finit 
>> precision problem. Ben pointed to it and you could learn about it from FAQ 
>> 7.31.
>>
>> Cheers
>> Petr
>>
>>> -Original Message-
>>> From: POLWART, Calum (COUNTY DURHAM AND DARLINGTON NHS FOUNDATION
>>> TRUST) 
>>> Sent: Thursday, January 17, 2019 2:56 PM
>>> To: PIKAL Petr ; Ben Tupper
>>> 
>>> Cc: r-help@r-project.org
>>> Subject: RE: [R] I can't get seq to behave how I think it should
>>>
>>> Thanks guys.
>>>
>>> I've used Petr's method and its working for me.
>>>
>>> If the data had been from a calculation I'd have rounded it... just
>>> didn't expect seq to break it!
>>>
>>> C
>>>
>>> -Original Message-
>>> From: PIKAL Petr [mailto:petr.pi...@precheza.cz]
>>> Sent: 17 January 2019 13:53
>>> To: Ben Tupper; POLWART, Calum (COUNTY DURHAM AND DARLINGTON NHS
>>> FOUNDATION TRUST)
>>> Cc: r-help@r-project.org
>>> Subject: RE: [R] I can't get seq to behave how I think it should
>>>
>>> Hi
>>>
>>> Or you could use rounding.
>>> which(round(lut, 3)==1.8)
>>> [1] 401
>>>
>>> Cheers
>>> Petr
>>>
 -Original Message-
 From: R-help  On Behalf Of Ben Tupper
 Sent: Thursday, January 17, 2019 2:43 PM
 To: POLWART, Calum (COUNTY DURHAM AND DARLINGTON NHS
>>> FOUNDATION TRUST)
 
 Cc: r-help@r-project.org
 Subject: Re: [R] I can't get seq to behave how I think it should

 Hi,

 This looks like a floating point reality bump - see

 https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-thin
 k- these-numbers-are-equal_003f
 

 You can use other methods to finding your row - I would opt for
 findInterval()

> lut = seq(1.4, 2.1, by=0.001)
> findInterval(1.8, lut)
 [1] 401

 findInterval() uses a rapid search to find the index in the look up
 table (lut) that is just less than  or equal to the search value (in
 your example
>>> 1.8).

 Cheers,
 Ben

> On Jan 17, 2019, at 8:33 AM, POLWART, Calum (COUNTY DURHAM AND
 DARLINGTON NHS FOUNDATION TRUST) via R-help 
 wrote:
>
> I am using seq with the expression seq(1.4, 2.1, by=0.001) to
> create a sequence of references from 1.4 to 2.1 in 0.001
> increments.  They appear to be created correctly.  They have a
> related pair of data which for the purposes of this we will call
> val.  I'm interested in the content on the row with seq = 1.8. But
> I can't seem to get it returned.  I can get other values but not
> 1.8!  yet looking at row
> 401 there is nothing to indicate an issue
>
>> a = 1.4
>> b = 2.1
>> seq = seq(a, b, by=0.001)
>> val = ceiling(seq * 50)
>> s=data.frame(seq, val)
>> s$val[seq==1.799]
> [1] 90
>> s$val[s$seq==1.8]
> numeric(0)
>> s$val[seq==1.8]
> numeric(0)
>> s$val[s$seq==1.800]
> numeric(0)
>> s$val[s$seq==1.801]
> [1] 91
>> head(s[s$seq>1.798,])
> seq val
> 400 1.799  90
> 401 1.800  90
> 402 1.801  91
> 403 1.802  91
> 404 1.803  91
> 405 1.804  91
>
>
> Can anyone explain what's 

Re: [R] Error in nlsModel

[ProfJCNash]

I'm not the author of nlsModel, so would prefer not to tinker with it.

But "singular gradient" is a VERY common problem with nls() that is used
by nlsModel as I understand it. The issue is actually a singular
Jacobian matrix resulting from a rather weak approximation of the
derivatives (a simple forward approximation as far as I can determine,
and using a fairly large step (1e-7), though a choice I'd probably make
too for this approach.

Duncan Murdoch and I wrote nlsr to do analytic derivatives where
possible. If you can use that (i.e., extract the modeling part of
nlsModel and call nlxb() from nlsr directly), I suspect you will have
better luck. If you still get singularity, it likely means that you
really have parameters that are some combination of each other.

JN

On 2018-10-08 05:14 AM, Belinda Hum Bei Lin wrote:
> Hello,
> 
> It is my first time using R studio and I am facing the error of
> "Error in nlsModel(formula, mf, start, wts) :
>   singular gradient matrix at initial parameter estimates"
> when I try to run my script. From what I read online, I understand that the
> error might be due to the parameters. However, I do not know how to choose
> the right set of parameters. Is there anyone who could advice me on how to
> do this?
> 
> Below are my script details:
> rm(list=ls()) #remove ALL objects
> cat("\014") # clear console window prior to new run
> Sys.setenv(LANG = "en") #Let's keep stuff in English
> Sys.setlocale("LC_ALL","English")
> 
> ##
> #import necessary packages
> #
> 
> ##To install the packages use the function install.packages. Installing is
> done once.
> #install.packages("ggplot2")
> #install.packages("minpack.lm")
> #install.packages("nlstools")
> 
> ##Activate the packages. This needs to be done everytime before running the
> script.
> require(ggplot2)
> require(minpack.lm)
> require(nlstools)
> 
> 
> 
> #
> #define the Weibull function
> #
> Weibull<-function(tet1, tet2,x){
>   1-exp(-exp(tet1+tet2*log10(x)))
> }
> 
> #
> ##define the inverse of the Weibull function. put in effect and get
> concentration as output
> #
> iWeibull<-function(tet1,tet2,x){
>   10^((log(-log(1-x))-tet1)/tet2)
> }
> 
> 
> #
> #define the Logit function
> #
> Logit<-function(tet1, tet2,x){
>   1/(1+exp(-tet1-tet2*log10(x)))
> }
> 
> #
> ##define the inverse of the Logit function
> #
> iLogit<-function(tet1,tet2,x){
>   10^(-(log(1/x-1)+tet1)/tet2)
> }
> 
> #
> #define the Probit function
> #
> Probit<-function(tet1, tet2, x){
>   pnorm(tet1+tet2*(log10(x)))
> }
> 
> #
> ##define the inverse of the Probit function
> #
> iProbit<-function(tet1,tet2,x){
>   10^((qnorm(x)-tet1)/tet2)
> }
> 
> #
> # Establish data to fit
> # data given here are the data for Diuron from the example datasets
> #
> # Of course one could also import an external datafile via e.g.
> # read.table, read.csv functions
> 
> ### example to choose a file for import with the read.csv function, where
> "," is seperating the columns,
> # header=TURE tells R that the first row contains the titles of the
> columns, and
> # stringsAsFactors = FALSE specify that the characters should not be
> converted to factors. For more info run ?read.csv
> effectdata<-read.csv(file.choose(),sep=",",stringsAsFactors = FALSE,header
> = TRUE)
> ?read.csv
> ###
> 
> #
> conc<-c(0,
> 0,
> 0,
> 0,
> 0,
> 0,
> 0.000135696,
> 0.000135696,
> 0.000135696,
> 0.000152971,
> 0.000152971,
> 0.000152971,
> 0.000172445,
> 0.000172445,
> 0.000172445,
> 0.000194398,
> 0.000194398,
> 0.000194398,
> 0.000219146,
> 0.000219146,
> 0.000219146,
> 0.000247044,
> 0.000247044,
> 0.000247044
> )
> 
> effect<-c(5.342014355,
>   13.46249176,
>   -9.249022885,
>   -6.666486351,
>   1.00292152,
>   -3.891918402,
>   12.63136345,
>   -2.372582186,
>   8.601073479,
>   1.309926638,
>   0.772728968,
>   -7.01067202,
>   30.65306236,
>   28.10819667,
>   17.94875421,
>   73.00440617,
>   71.33593917,
>   62.23994217,
>   99.18897648,
>   99.05982514,
>   99.2325145,
>   100.2402872,
>   100.1276669,
>   100.1501468
> )
> 
> #build input dataframe
> effectdata<-data.frame(conc,effect)
> 
> #plot the data just to get a first glance of the data
> ggplot()+
>   geom_point(data=effectdata,aes(x=conc,y=effect), size = 5)+
>   scale_x_log10("conc")
> 
> 
> #delete controls
> effectdata_without_controls<-subset(effectdata,effectdata$conc>0)
> 
> 
> #save controls in a seperate dataframe called effectdata_control, which
> will be added to the ggplot in the end.
> #since you can't have 0 on a logscale we 

Re: [R] Maximum number of iterations (maxit) does not work in hydroPSO-modFit-optimr in r

[ProfJCNash]

As the author of optimx (there is a new version just up), I can say that
maxit only "works" if the underlying solver is set up to use it.

You seem to be mistaken, however, in this case, as the following example
shows.

> library(optimx)
> library(adagio)
> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 18.3

Matrix products: default
BLAS: /usr/lib/openblas-base/libblas.so.3
LAPACK: /usr/lib/libopenblasp-r0.2.18.so

locale:
 [1] LC_CTYPE=en_CA.UTF-8   LC_NUMERIC=C
LC_TIME=en_CA.UTF-8
 [4] LC_COLLATE=en_CA.UTF-8 LC_MONETARY=en_CA.UTF-8
LC_MESSAGES=en_CA.UTF-8
 [7] LC_PAPER=en_CA.UTF-8   LC_NAME=C  LC_ADDRESS=C

[10] LC_TELEPHONE=C LC_MEASUREMENT=en_CA.UTF-8
LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] adagio_0.7.1 optimx_2018-7.10

loaded via a namespace (and not attached):
[1] compiler_3.5.1tools_3.5.1   numDeriv_2016.8-1
> test <- optimr(c(-1.2, 1), fn=fnRosenbrock, method="L-BFGS-B",
control=list(maxit=2, trace=2))
Parameter scaling:[1] 1 1
N = 2, M = 5 machine precision = 2.22045e-16
This problem is unconstrained.
final  value 4.120516
stopped after 3 iterations

So it seems to have worked, albeit one iteration later than specified
(the test is not always immediate).

Note optimr() is just a wrapper. There may be improvements that could be
made to use maxit and similar controls. Note that the code is in plain R
and optimr() was written in a way that allows it to be maintained. The
most difficult problem -- which may be related to the OP's complaint in
some cases -- is translating maxit to whatever is used by the solver,
then copying back the relevant information about number of iterations.

John Nash


On 2018-10-03 01:30 PM, Ahmed Attia wrote:
> The argument maxit used in libraries optimr, hydroPSO, and FME to
> control the maximal number of iterations to be performed does not work
> at ALL!!
> 
> This needs to be fixed...
> 
> 
> 
> mysamp <- function(n, m, s, lwr, upr, nnorm) {
>   samp <- rnorm(nnorm, m, s)
>   samp <- samp[samp >= lwr & samp <= upr]
>   if (length(samp) >= n) {
> return(sample(samp, n))
>   }
>   stop(simpleError("Not enough values to sample from. Try increasing nnorm."))
> }
> 
> x <- mysamp(1000,0.8117648,0.1281839,0,1,100)
> y <- mysamp(1000,0.7530346,0.1865946,0,1,100)
> 
> n <- length(x)
> xgroup <- c('a','b')
> ygroup <- c('c','d')
> 
> for (i in c(2:n)){
>   xgroup[i] <-  c('b')
>   ygroup[i] <-  c('d')
> }
> 
> dataset <- data.frame(x = x, y = y, xgroup = as.factor(xgroup),
>   ygroup = as.factor(ygroup))
> dataset$ObsNo <- 1:n
> dataset <- dataset[order(dataset$x,decreasing=F),]
> 
> xopt <- c(dataset$x[1],0.95)#Initial critical x,y values
> 
> 
> my.function <- function(xopt){
>   for (i in c(1:n)){
> dataset$xgroup[i] <- if(dataset$x[i] < xopt[1]) 'a' else 'b'
> dataset$ygroup[i] <- if(dataset$y[i] < xopt[2]) 'c' else 'd'
> dataset$q.i[i] <- with(dataset, ifelse
>(dataset$xgroup[i]=='a' &
> dataset$ygroup[i]=='d', 1, 0))
> dataset$q.ii[i] <- with(dataset, ifelse
> (dataset$xgroup[i]=='b' &
> dataset$ygroup[i]=='d', 1, 0))
> dataset$q.iii[i] <- with(dataset, ifelse
>  (dataset$xgroup[i]=='b' &
> dataset$ygroup[i]=='c', 1, 0))
> dataset$q.iv[i] <- with(dataset, ifelse
> (dataset$xgroup[i]=='a' &
> dataset$ygroup[i]=='c', 1, 0))
> 
> dataset$q.err[i] <- sum(dataset$q.i[i] + dataset$q.iii[i])
>   }
>   min.qerr <- sum(dataset$q.err)
>   q.I <- sum(dataset$q.i)
>   q.II <- sum(dataset$q.ii)
>   q.III <- sum(dataset$q.iii)
>   q.IV <- sum(dataset$q.iv)
>   q.Iandq.III <- sum(dataset$q.err)
>   print(c(q.I,
>   q.II,
>   q.III,
>   q.IV,
>   q.Iandq.III))
>   return(min.qerr)
> }
> 
> my.function(xopt)
> 
> #---Algorithm
> 
> xmin=c(0,0.60)
> xmax=c(0.95,0.95)
> 
> res= hydroPSO(par=xopt,my.function,method="spso2011",
>   lower=xmin,upper=xmax,control=list(maxit=3))
> 
> 
> res= modFit(my.function,p=xopt,method="Pseudo",
> lower=xmin,upper=xmax,control=list(numiter=3))
> 
> res= optimr(par=xopt,my.function,method="L-BFGS-B",
> lower=xmin,upper=xmax,control=list(maxit=3,trace=T))
> 
> #Only the hjkb form dfoptim that has argument maxfeval and it works.
> 
> res=hjkb(par=xopt,my.function,
>  lower=xmin,upper=xmax,control=list(maxfeval=3,tol=2^-10,info=T))
> 
> 
> 
> 
> 
> 
> Ahmed Attia, Ph.D.
> Agronomist & Crop Modeler
> 
> __
> 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, 

Re: [R] Optim function returning always initial value for parameter to be optimized

[ProfJCNash]

Did you check the gradient? I don't think so. It's zero, so of course
you end up where you start.

Try

data.input= data.frame(state1 = (1:500), state2 = (201:700) )
err.th.scalar <- function(threshold, data){

state1 <- data$state1
state2 <- data$state2

op1l <- length(state1)
op2l <- length(state2)

op1.err <- sum(state1 <= threshold)/op1l
op2.err <- sum(state2 >= threshold)/op2l

total.err <- (op1.err + op2.err)

return(total.err)
}

soln <- optim(par = 300, fn=err.th.scalar, data = data.input, method =
"BFGS")
soln
require("numDeriv")
gtest <- grad(err.th.scalar, x=300, data = data.input)
gtest


On 2018-02-09 09:05 AM, BARLAS Marios 247554 wrote:
> data.input= data.frame(state1 = (1:500), state2 = (201:700) )
> 
> with data that partially overlap in terms of values. 
> 
> I want to minimize the assessment error of each state by using this function:
> 
> err.th.scalar <- function(threshold, data){
>   
>   state1 <- data$state1
>   state2 <- data$state2
>   
>   op1l <- length(state1)
>   op2l <- length(state2)
>   
>   op1.err <- sum(state1 <= threshold)/op1l
>   op2.err <- sum(state2 >= threshold)/op2l
>   
>   total.err <- (op1.err + op2.err)
> 
>   return(total.err)
> }
> 
> 
> SO I'm trying to minimize the total error. This Total Error should be a U 
> shape essentially.
> 
> 
> I'm using optim as follows: 
> 
> optim(par = 300, fn=err.th.scalar, data = data.input, method = "BFGS")


Maybe develop an analytic gradient if it is very small, as the numeric
approximation can then be zero even when the true gradient is not.

JN

__
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] MCMC Estimation for Four Parametric Logistic (4PL) Item Response Model

[ProfJCNash]

If you have the expression of the model, package nlsr should be able to
form the Jacobian analytically for nonlinear least squares. Likelihood
approaches allow for more sophisticated loss functions, but the
optimization is generally much less reliable because one is working
with essentially squared quantities and possibly multiple minima where
some are not the ones you want.

JN

On 2018-01-18 12:32 PM, Doran, Harold wrote:
> I know of no existing functions for estimating the parameters of this model 
> using MCMC or MML. Many years ago, I wrote code to estimate this model using 
> marginal maximum likelihood. I wrote this based on the using nlminb and 
> gauss-hermite quadrature points from statmod. 
> 
> I could not find that code to share with you, but I do have code for 
> estimating the 3PL in this way and you could modify the likelihood for the 
> upper asymptote yourself.
> 
> 
> 
> -Original Message-
> From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of ALYSSA FATMAH 
> MASTURA
> Sent: Thursday, January 18, 2018 10:15 AM
> To: r-help@r-project.org
> Subject: [R] MCMC Estimation for Four Parametric Logistic (4PL) Item Response 
> Model
> 
> Good day Sir/Ma'am! This is Alyssa Fatmah S. Mastura taking up Master of 
> Science in Statistics at Mindanao State University-Iligan Institute 
> Technology (MSU-IIT), Philippines. I am currently working on my master's 
> thesis titled "Comparing the Three Estimation Methods for the Four Parametric 
> Logistic (4PL) Item Response Model". While I am looking for a package about 
> Markov chain Monte Carlo (MCMC) method for 4PL model in this R forum, I found 
> an sirt package of an MCMC method but only for the three parametric normal 
> ogive (3PNO) item response model. However, my study focus on 4PL model. In 
> line with this, I would like to know if there exist a function of MCMC method 
> for 4PL model in R language. I am asking for your help to inform me if such 
> function exist. I highly appreciate your response on this matter. Thank you 
> so much. Have a great day ahead!
> 
> 
> 
> Alyssa
> 
> 
> Virus-free.
> www.avast.com
> 
> <#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
> 
> --
> ---
> *DISCLAIMER AND CONFIDENTIALITY NOTICE* The Mindanao Sta...{{dropped:9}}
> 
> __
> 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] Finding optim.R function

[ProfJCNash]

The codes were taken from the 2nd edition of my book Compact Numerical
Methods for Computers, where they are in Pascal. They were converted by
p2c to c, so are pretty opaque and likely difficult to modify. Moreover,
they are based on 1970s codes I wrote for the first edition. Why not
look at optimr (CRAN) or the more extensive optimrx (R-forge) where
there are calls to pure R versions with improvements in the codes as
well as bounds constraints on parameters for some. If you have
suggestions or queries about the newer codes, contact me off-list and
we'll see what can be done.

JN (who will be at UseR! next week)


On 2017-06-27 12:46 PM, Tauras Vilgalys wrote:
> Hello, could anybody direct me where to find code for optim.R? I was able to 
> find the C code at http://docs.rexamine.com/R-devel/optim_8c.html, but the R 
> version would be easier for me to work with and modify.
> 
> 
> Thank you!
> 
> 
>   [[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] installing rgl

[ProfJCNash]

I had this problem this week in Linux Mint (debian/ubuntu based) and 
needed to install some libgl* and libglu* packages. A search for
"rgl  glu.h" and "rgl gl.h" found the appropriate suggestions, though I 
probably installed a couple of unnecessary packages too.


JN

On 2017-01-11 03:14 PM, Duncan Murdoch wrote:

On 11/01/2017 3:03 PM, Weiner, Michael wrote:

-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com]
Sent: Wednesday, January 11, 2017 2:55 PM
To: Weiner, Michael ; r-help@r-project.org
Subject: Re: [R] installing rgl


On this page

http://forums.fedoraforum.org/showthread.php?t=294543

eventually it turned out that a similar problem was fixed by

yum install libpng-devel


Thank you for your response Duncan, unfortunately that didn't help,
though I do see in config.log:

configure:4429: checking for glEnd in -lGL
configure:4454: gcc -o conftest -g -O2  -DHAVE_PNG_H
-I/usr/include/libpng16  conftest.c -lGL   -L/usr/lib64 -lpng16 -lX11 >&5
/usr/bin/ld: skipping incompatible
/usr/lib/gcc/x86_64-redhat-linux/6.2.1/../../../libGL.so when
searching for -lGL
/usr/bin/ld: skipping incompatible /lib/libGL.so when searching for -lGL
/usr/bin/ld: skipping incompatible /usr/lib/libGL.so when searching
for -lGL
/usr/bin/ld: cannot find -lGL

So something else is up



I don't know Fedora at all so I don't know what you'd need to do this,
but I'd suggest asking to uninstall and reinstall mesa-libGL-devel and
mesa-libGLU-devel (and maybe libpng-devel).

Duncan Murdoch


Michael


===


 Please consider the environment before printing this e-mail

Cleveland Clinic is ranked as one of the top hospitals in America by
U.S.News & World Report (2015).
Visit us online at http://www.clevelandclinic.org for a complete
listing of our services, staff and locations.


Confidentiality Note:  This message is intended for use ...{{dropped:18}}

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

Re: [R] Linear optimization with quadratic constraints

[ProfJCNash]

Small example code to set up the problem?

JN

On 2017-01-07 06:26 AM, Preetam Pal wrote:
> Hi Guys,
> Any help with this,please?
> Regards,
> Preetam
> 
> On Thu, Jan 5, 2017 at 4:09 AM, Preetam Pal  wrote:
> 
>> Hello guys,
>>
>> The context is ordinary multivariate regression with k (>1) regressors,
>> i.e. *Y = XB + Error*, where
>> Y = n X 1 vector of predicted variable,
>> X = n X (k + 1) matrix of regressor variables(including ones in the first
>> column)
>> B = (k+1) vector of coefficients, including intercept.
>>
>> Say, I have already estimated B as B_hat = (X'X)^(-1) X'Y.
>>
>> I have to solve the following program:
>>
>> *minimize f(B) = LB*   ( L is a fixed vector 1 X (k+1)   )
>> such that:
>> *[(B-B_hat)' * X'X * (B-B_hat) ] / [ ( Y - XB_hat)' (Y - XB_hat) ] *  is
>> less than a given value *c*.
>>
>> Note that this is a linear optimization program *with respect to B* with
>> quadratic constraints.
>>
>> I don't understand how we can solve this optimization - I was going
>> through some online resources, each of which involve manually computing
>> gradients of the objective as well as constraint functions - which I want
>> to avoid (at least manually doing this).
>>
>>
>> Can you please help with solving this optimization problem? The inputs
>> would be:
>>
>>- X and Y
>>- B_hat
>>- L
>>- c
>>
>>
>> Please let me know if any further information is required - the set-up is
>> pretty general.
>>
>> Regards,
>> Preetam
>>
> 
> 
>

__
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] Dealing with -Inf in a maximisation problem.

[ProfJCNash]

Rolf, What optimizers did you try? There are a few in the optimrx package on 
R-forge that handle bounds, and it may be
useful to set bounds in this case. Transformations using log or exp can be 
helpful if done carefully, but as you note,
they can make the function more difficult to optimize.

Be cautious about using the default numerical gradient approximations. optimrx 
allows selection of the numDeriv grad()
function, which is quite good. Complex step would be better, but you need a 
function which can be computed with complex
arguments. Unfortunately, numerical gradients often step over the cliff edge of 
computability of the function. The
bounds are not checked for the step to compute things like (f(x+h) - f(x) / h.

Cheers, JN

On 16-11-06 07:07 PM, William Dunlap via R-help wrote:
> Have you tried reparameterizing, using logb (=log(b)) instead of b?
> 
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
> 
> On Sun, Nov 6, 2016 at 1:17 PM, Rolf Turner  wrote:
> 
>>
>> I am trying to deal with a maximisation problem in which it is possible
>> for the objective function to (quite legitimately) return the value -Inf,
>> which causes the numerical optimisers that I have tried to fall over.
>>
>> The -Inf values arise from expressions of the form "a * log(b)", with b =
>> 0.  Under the *starting* values of the parameters, a must equal equal 0
>> whenever b = 0, so we can legitimately say that a * log(b) = 0 in these
>> circumstances.  However as the maximisation algorithm searches over
>> parameters it is possible for b to take the value 0 for values of
>> a that are strictly positive.  (The values of "a" do not change during
>> this search, although they *do* change between "successive searches".)
>>
>> Clearly if one is *maximising* the objective then -Inf is not a value of
>> particular interest, and we should be able to "move away".  But the
>> optimising function just stops.
>>
>> It is also clear that "moving away" is not a simple task; you can't
>> estimate a gradient or Hessian at a point where the function value is -Inf.
>>
>> Can anyone suggest a way out of this dilemma, perhaps an optimiser that is
>> equipped to cope with -Inf values in some sneaky way?
>>
>> Various ad hoc kludges spring to mind, but they all seem to be fraught
>> with peril.
>>
>> I have tried changing the value returned by the objective function from
>> "v" to exp(v) --- which maps -Inf to 0, which is nice and finite. However
>> this seemed to flatten out the objective surface too much, and the search
>> stalled at the 0 value, which is the antithesis of optimal.
>>
>> The problem arises in a context of applying the EM algorithm where the
>> M-step cannot be carried out explicitly, whence numerical optimisation.
>> I can give more detail if anyone thinks that it could be relevant.
>>
>> I would appreciate advice from younger and wiser heads! :-)
>>
>> 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/posti
>> ng-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

Re: [R] nls.lm

[ProfJCNash]

Berend's point is well-taken. It's a lot of work to re-jig a code, especially 
one more than
30 years old.

On the other hand, nlmrt is all-R, and it does more or less work on 
underdetermined systems as I
illustrated in a small script. The changes needed to treat the problem as Mike 
suggests are
pretty easy, but they imply that some aspects of the modeling computations have 
to be omitted.
Mike: If you give that a try, I'd be interested in outcomes, and am willing to 
answer questions
on how the code works.

Indeed, I want to include the sort of "diagnostics" that allow underdetermined 
systems
in my packages, so that some results are replaced with NA or NULL and helpful 
warnings rather
than "Failed m 
>> On 21 Oct 2016, at 06:00, Mike meyer <1101...@gmx.net> wrote:
>>
>> Let's take a different view of the problem.
>> Given f=(f_1,...,f_m):R^n -> R^m we want to minimize ||f(x)||.
>>
>> What distinguishes this from a general minimization problem is that you know 
>> the structure of the
>> objective function F(x)=||f(x)||² and have the individual constituents f_j.
>> Make use of that information as appropriate.
>>
>> This is more general than trying to solve the system f(x)=0 or fitting a 
>> model to data.
>> In this more general setting notions such as underdetermined/overdetermined 
>> system do not apply.
>>
>> The restricted view of model fitting serves only to confuse the issue.
>> For that reason it is (in my view) a bad idea to force the user to set up 
>> his problem in 
>> R-model notation.
>>
> 
> I assume that you have been referring to the R package minpack.lm.
> 
> I've had a look at the underlying Fortran code (from Minpack and developed by 
> More et.al.  made in a distant past) as used by the package.
> That underlying code returns an error when the condition:  number of 
> functions (m)  >=  the number of independent variables (n)
> is not satisfied i.e. when m < n.
> 
> Making that more general would entail a lot of thinking and reworking of the 
> code. As far as I can see it is not possible to just remove the condition 
> m>=n from the underlying Fortran. More (possibly many) changes would be 
> required. Blaming R and/or the package author/maintainer is unfair.
> 
> If you require a more general version of the algorithm or if you want 
> something else you will have to roll your own package/code.
> If you don't feel that minpack.lm is appropriate for your application and you 
> want changes you'll have to discuss matters with Moré 
> (http://www.mcs.anl.gov/~more/) if I got the correct link.
> 
> Berend
>

__
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] Optimize selected variables but not all - nloptr

[ProfJCNash]

optimr/optimrx are for unconstrained and bounds constrained problems.

I think you are going to have to do the masking yourself. It's more messy and 
tedious than
difficult. The optimization is on a new set of parameters newpar that has
fewer components, so in your objective function and constraints you have to 
expand it to the
full set of parameters, putting in the fixed values. A bit of back and forth, 
but
it does work.

However, do test, test, test to make sure things are working before you try the
optimization. It is REALLY easy to make a silly mistake. I see someone who has
made lots every time I look in a mirror.

Best,

JN



On 16-10-20 11:59 AM, Narendra Modi wrote:
> Thanks Prof Nash.
> The reason I used nlopr() in my problem is due to non linear
> constraints. I wonder if optimrx/optim can model the below scenario. I
> will be elated if it can.
> The problem in hand goes like this:
> 
> There are 2 injectors and 2 producers. Consider these as some entity.
> 
> I have an actual dataset, Act.Matrix; assume row=50, col=2.
> The predicted dataset is calculated by using the attached equation,
> i.e for each row of the prediction, that equation is used.
> 
> The parameters to be evaluated are Tau, and a set of Fij's such that
> sum of Fij's <=1 .
> Fij represents connectivity between each inj to each prod, hence for
> this example there will be a total of 4 fij. and 2 tau (one tau for
> each producer)
> 
> The lb for each parameter that I provided is 0 and ub as Inf
> 
> TAU fijfij2
>>>  [1,] 14.3  0.010.01449572
>>>  [2,] 14.3  0.2  0.
> 
> The constraints being Sum of each column of Fij<=1.
> 
> With nloptr (algorithm NLOPT_LN_COYLA) , I am able to solve it to some
> extent. but when the dataset is scaled (consider having 50 producers
> and 25 injectors!! ), the number of variables to be solved is much
> higher!
> 
> Do you think OPTIMX/OPTIMRX can handle non linear constraints like
> that? if it can then I can definitely use the MASKED parameters.
> 
> NM
> 
> 
> On Wed, Oct 19, 2016 at 8:48 PM, ProfJCNash <profjcn...@gmail.com> wrote:
>> I refer to such parameters as "masked" in my 2014 book Nonlinear parameter 
>> optimization with R tools.
>> Recently I put package optimrx on R-forge (and optimr with fewer solvers on 
>> CRAN) that allows for masks with
>> all the parameters. The masks can be specified as you suggest with 
>> start=lower=upper. However, for reasons
>> I won't go into here, nloptr solvers are not yet included. However, I 
>> suspect either Rvmmin or Rcgmin will
>> work fine.
>>
>> I will guess that nloptr does NOT cater for masks.
>>
>> JN
>>
>> On 16-10-19 05:59 PM, Narendra Modi wrote:
>>> Hello All,
>>> I have a matrix with initial values as below and I need to optimize
>>> the variables that are greater than 0.
>>>
>>>   TAU   fij  fij2
>>>  [1,] 14.33375 0.000 0.01449572
>>>  [2,] 14.33375 0.000 0.
>>>  [3,] 14.33375 0.000 0.
>>>  [4,] 14.33375 0.000 0.02206446
>>>  [5,] 14.33375 0.000 0.
>>>  [6,] 14.33375 0.000 0.
>>>  [7,] 14.33375 0.000 0.
>>>  [8,] 14.33375 0.8279846 0.
>>>  [9,] 14.33375 0.000 0.03695833
>>> [10,] 14.33375 0.000 0.
>>>
>>> Or structure(c(14.3337481730129, 14.3337481730129, 14.3337481730129,
>>> 14.3337481730129, 14.3337481730129, 14.3337481730129, 14.3337481730129,
>>> 14.3337481730129, 14.3337481730129, 14.3337481730129, 0, 0, 0,
>>> 0, 0, 0, 0, 0.827984553120177, 0, 0, 0.0144957197835888, 0, 0,
>>> 0.0220644627842788, 0, 0, 0, 0, 0.0369583294835073, 0), .Dim = c(10L,
>>> 3L), .Dimnames = list(NULL, c("TAU", "fij", "fij2")))
>>>
>>> Is it possible to provide lowerbound and upperbound as 0 for variables
>>> (< 0 in the initial matrix) and nloptr will consider them "unchanged"
>>> during optimization?
>>>
>>> Rstudio crashes when I try to do that. Is this a bug or I should
>>> approach it differently?
>>>
>>> NM
>>>
>>> __
>>> 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-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] nls.lm

[ProfJCNash]

>From a statistician's point of view, "nonsense" may be OK, but there are other 
>applications of R where
(partial or non-unique) solutions may be needed.

Yesterday I raised the question of how nonlinear least squares could be adapted 
to underdetermined problems.
Many folk are unaware of such possibilities, even in the linear case. In the 
nonlinear case, the m>n condition is
far from adequate to provide appropriate warnings to users. What should we do 
to detect non-unique solutions?

My interest is in building better nonlinear least squares and optimization 
software, and I find software users
(including R users) have some belief that those of us building the tools they 
use are magic wizards or fairy
godmothers who have provided everything they need in bullet-proof packages. 
This is so far from reality. All of
the tools I have worked with have weaknesses. To close the loopholes we need 
those "small reproducible examples",
including cases like these, so that better diagnostics can be devised.

We also need discussion on what to present as diagnostics and how to do so for 
maximum benefit and least
"get in the way". A two-way communication with users can aid immensely.

Sorry for the rant, but I'm in the midst of trying to prepare a unification of 
nls/nlmrt/minpack.lm and some
of the effort is pretty messy, especially in the area of derivatives and 
diagnostics.

Below is a little script that tries Berend's problem. The nonlinear least 
squares "runs" but the output fails
except for str(). The script will stop on failure at some points, so you need 
to paste some statements. I welcome
similar scripts/examples to build the necessary tests for improved packages.

JN

Here's the script.

# try to solve undetermined system by nonlinear least squares
X1 <- c(1,1)
Y1 <- c(1,1)
Z1 <- c(1,1)
RHS1 <- c(3,4)
X2 <- c(1,2)
RHS2 <- c(3,6)
mydata <- data.frame(X1, X2, Y1, Z1, RHS1, RHS2)
require(nlmrt)
st1 <- c(px=1,py=1,pz=1)
st0 <- c(px=0,py=0,pz=0)
sol10 <- nlxb(RHS1 ~ px*X1 + py*Y1 + pz*Z1, data=mydata, trace=1, start=st0)
summary(sol10)
print(sol10)
str(sol10)
sol11 <- nlxb(RHS1 ~ px*X1 + py*Y1 + pz*Z1, data=mydata, trace=1, start=st1)
summary(sol11)
print(sol11)
str(sol11)
## try RHS2
sol20 <- nlxb(RHS2 ~ px*X1 + py*Y1 + pz*Z1, data=mydata, trace=1, start=st0)
summary(sol20)
print(sol20)
str(sol20)
sol21 <- nlxb(RHS2 ~ px*X1 + py*Y1 + pz*Z1, data=mydata, trace=1, start=st1)
summary(sol21)
print(sol21)
str(sol21)
# change first column -- then we get solutions
sol220 <- nlxb(RHS2 ~ px*X2 + py*Y1 + pz*Z1, data=mydata, trace=1, start=st0)
summary(sol220)
print(sol220)
str(sol220)
sol221 <- nlxb(RHS2 ~ px*X2 + py*Y1 + pz*Z1, data=mydata, trace=1, start=st1)
summary(sol221)
print(sol221)
str(sol221)



On 16-10-20 07:05 AM, S Ellison wrote:
>> How do you reply to a specific post on this board instead of the thread?
> 
> You can reply to the individual, as I just did.
> 
> But I strongly suggest that you don't. You would be much better advised to 
> discontinue debate and follow the essential advice given by nls.lm, which - 
> no matter whether couched in terms of count of residuals - is simply to make 
> sure that you have more independent data than variables when seeking a unique 
> numerical solution by non-linear least squares. If you don't you'll get 
> nonsense.
> 
> 
> S Ellison
> 
> 
> 
> 
> 
> ***
> This email and any attachments are confidential. Any use...{{dropped:8}}
> 
> __
> 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] Optimize selected variables but not all - nloptr

[ProfJCNash]

I refer to such parameters as "masked" in my 2014 book Nonlinear parameter 
optimization with R tools.
Recently I put package optimrx on R-forge (and optimr with fewer solvers on 
CRAN) that allows for masks with
all the parameters. The masks can be specified as you suggest with 
start=lower=upper. However, for reasons
I won't go into here, nloptr solvers are not yet included. However, I suspect 
either Rvmmin or Rcgmin will
work fine.

I will guess that nloptr does NOT cater for masks.

JN

On 16-10-19 05:59 PM, Narendra Modi wrote:
> Hello All,
> I have a matrix with initial values as below and I need to optimize
> the variables that are greater than 0.
> 
>   TAU   fij  fij2
>  [1,] 14.33375 0.000 0.01449572
>  [2,] 14.33375 0.000 0.
>  [3,] 14.33375 0.000 0.
>  [4,] 14.33375 0.000 0.02206446
>  [5,] 14.33375 0.000 0.
>  [6,] 14.33375 0.000 0.
>  [7,] 14.33375 0.000 0.
>  [8,] 14.33375 0.8279846 0.
>  [9,] 14.33375 0.000 0.03695833
> [10,] 14.33375 0.000 0.
> 
> Or structure(c(14.3337481730129, 14.3337481730129, 14.3337481730129,
> 14.3337481730129, 14.3337481730129, 14.3337481730129, 14.3337481730129,
> 14.3337481730129, 14.3337481730129, 14.3337481730129, 0, 0, 0,
> 0, 0, 0, 0, 0.827984553120177, 0, 0, 0.0144957197835888, 0, 0,
> 0.0220644627842788, 0, 0, 0, 0, 0.0369583294835073, 0), .Dim = c(10L,
> 3L), .Dimnames = list(NULL, c("TAU", "fij", "fij2")))
> 
> Is it possible to provide lowerbound and upperbound as 0 for variables
> (< 0 in the initial matrix) and nloptr will consider them "unchanged"
> during optimization?
> 
> Rstudio crashes when I try to do that. Is this a bug or I should
> approach it differently?
> 
> NM
> 
> __
> 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] nls.lm

[ProfJCNash]

I sometimes find it useful to use nonlinear least squares for fitting an 
approximation i.e., zero
residual model. That could be underdetermined.

Does adding the set of residuals that is the parameters force a minimum length 
solution? If the
equations are inconsistent, then the residuals apart from the parameters would, 
I believe, be
non-zero. Of course, if they turn out non-zero, it doesn't mean there is no 
solution -- the method
may have failed, but if they are zero, we have constructed a solution.

While I've used this approach, I'm thinking that it deserves closer examination 
and testing.

Best, JN


On 16-10-19 11:28 AM, Berend Hasselman wrote:
> 
>> On 19 Oct 2016, at 14:09, Mike meyer <1101...@gmx.net> wrote:
>>
>> @pd: you know that a System of equations with more variables than equations 
>> is always solvable
>> and if a unique solution is desired one of mimimal norm can be used.
>>
> 
> Not true.
> 
> Take the system with 3 variables and 2 equations
> 
> x+y+z = 3
> x+y+z = 4
> 
> This does not have a solution.
> See https://en.wikipedia.org/wiki/Consistent_and_inconsistent_equations
> 
> Berend
> 
>> According to "Methods for nonlinear least squares problems" by Madsen, 
>> Nielsen and Tingleff the LM-algorithm
>> solves Systems of the form 
>>[J(x)'J(x)+\mu*I]x=...
>> with \mu>0 so that the Matrix on the left is always positive definite, 
>> especially nonsingular.
>>
>> __
>> 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-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] nls.lm

[ProfJCNash]

Peter is right that the conditions may be embedded in the underlying code. (Ask 
Kate!)

My nlmrt package is all in R, so the conditions are visible. I'm currently in 
process of rejigging this
using some work Duncan Murdoch helped with a while ago (I've had some other 
things get in the way), so
I can change such conditions. In the present case, I think I'd issue a warning, 
but haven't checked
whether I do or not.

I'd be very happy to have the usual "minimal reproducible example" scripts of 
issues with any nonlinear
least squares problems to use as tests to keep knocking the rough edges off the 
codes I'm working on.
Note, in particular, that the new code (named nlsr) tries to use analytic 
derivatives that are computed
from the model expression. Ultimately would like to have automatic derivatives 
to apply to functions too,
and migrate that to optimrx/optimr which are recently on Rforge/CRAN (optimr 
has fewer solvers to avoid
the "your program doesn't work" when a solver is changed/removed).

Best, JN

On 16-10-19 08:09 AM, Mike meyer wrote:
> @pd: you know that a System of equations with more variables than equations 
> is always solvable
> and if a unique solution is desired one of mimimal norm can be used.
> 
> According to "Methods for nonlinear least squares problems" by Madsen, 
> Nielsen and Tingleff the LM-algorithm
> solves Systems of the form 
> [J(x)'J(x)+\mu*I]x=...
> with \mu>0 so that the Matrix on the left is always positive definite, 
> especially nonsingular.
> 
> __
> 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] Finding starting values for the parameters using nls() or nls2()

[ProfJCNash]

The key lines are

library(nlmrt)
test <- nlxb(expf2, start= c(b0=.1, b1=1, th=.1), trace=TRUE, data=cl)

Thus I started with .1 1 and .1. The "solution" from nlxb, which is using 
analytic derivatives
and a very aggressive Marquardt code to keep trying even in bad situations, was
as you included below. Note that the singular values of the Jacobian are given 
(they are
recorded on the same table as the parameters, but do NOT correspond to the 
parameters.
The placement was simply a tidy place to put these numbers.)

The ratio of these sv's is 1.735e+16/0.004635 or approx 4E+18, so the condition 
number
of the traditional Gauss Newton approach is about 1E+37. Not a nice problem!

You probably should reformulate.

JN




On 16-10-10 10:41 AM, Pinglei Gao wrote:
> Thanks very much for your kindness help. I run your script then came out
> lots of outputs and I also studied the solution you posted. Forgive my
> ignorance, I still can't find the suitable starting values. Did I
> misunderstand something?
> 
> Best,
> 
> Pinglei Gao
> 
> -邮件原件-
> 发件人: ProfJCNash [mailto:profjcn...@gmail.com] 
> 发送时间: 2016年10月10日 10:41
> 收件人: Gabor Grothendieck; Pinglei Gao
> 主题: Re: [R] Finding starting values for the parameters using nls() or
> nls2()
> 
> I forgot to post the "solution" found by nlmrt:
> 
> nlmrt class object: x
> residual sumsquares =  1086.8  on  15 observations
> after  5001Jacobian and  6991 function evaluations
>   namecoeff  SE   tstat  pval  gradient
> JSingval
> b05.3274e-14NA NA NA  -6.614e+13
> 1.735e+16
> b1   33.5574NA NA NA  -3.466
> 11518
> th   -0.00721203NA NA NA  -740.8
> 0.004635
> 
> 
> Note the singular values -- this is the worst SV(max)/SV(min) ratio I've
> observed!
> 
> JN
> 
> 
>

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

Re: [R] Finding starting values for the parameters using nls() or nls2()

[ProfJCNash]

Despite nlmrt "solving" the OP's problem, I think Gabor's suggestion likely 
gives a more sensible approach
to the underlying modelling problem.

It is, of course, sometimes important to fit a particular model, in which case 
nls2 and nlmrt are set up to grind away.
And hopefully the follow-up to nlmrt I'm working on will have enough capability 
in getting analytic derivatives
to work for a wider class of models. Note that functional approaches in nlmrt 
and minpack.lm allow users to
provide derivatives. Too many users think numerical approximations are a 
panacea, but my experience is that
most problems benefit from very accurate derivatives, of which analytic 
expressions are generally the best.

JN

On 16-10-09 09:24 PM, Gabor Grothendieck wrote:
> If you are not tied to that model the SSasymp() model in R could be
> considered and is easy to fit:
> 
> # to plot points in order
> o <- order(cl$Area)
> cl.o <- cl[o, ]
> 
> fm <- nls(Retention ~ SSasymp(Area, Asym, R0, lrc), cl.o)
> summary(fm)
> 
> plot(Retention ~ Area, cl.o)
> lines(fitted(fm) ~ Area, cl.o, col = "red")
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

Re: [R] Finding starting values for the parameters using nls() or nls2()

[ProfJCNash]

I didn't try very hard, but got a solution from .1, 1, .1 with nlxb() from 
nlmrt. It took a lot
of iterations and looks to be pretty ill-conditioned. Note nlmrt uses analytic 
derivatives if it
can, and a Marquardt method. It is designed to be a pit bull -- tenacious, not 
fast.

I'm working on a replacement for this and nls(), but it will be a while. 
However, I welcome
short scripts like this as tests. My script below

Best, JN

cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),

Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )

expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}
expf2 <- "Retention ~ exp(b0*exp(b1*Area^th))"

# grid.Disperse <- expand.grid(list(b0 = seq(0.01,4, by = 0.01), th =
#c(0.02),b1 = seq(0.01, 4, by = 0.01)))

#> Disperse.m2a <- nls2(Retention ~expFct(Area, b0, b1,th), data = cl, start
#= grid.Disperse, algorithm = "brute-force")

# Disperse.m2a

library(nlmrt)
test <- nlxb(expf2, start= c(b0=.1, b1=1, th=.1), trace=TRUE, data=cl)



On 16-10-09 07:40 PM, peter dalgaard wrote:
> 
>> On 10 Oct 2016, at 00:40 , Bert Gunter  wrote:
>>
>> Well... (inline -- and I hope this isn't homework!)
>>
> 
> Pretty much same as I thought. 
> 
> Fixing th=0.02 in the grid search looks wrong. Bert's plot is pretty linear, 
> so th=1 is a good guesstimate. There's a slight curvature but to reduce it, 
> you would increase th, not decrease it. Running the regression, as Bert 
> suggests, indicates that b0=5.16 and b1= -0.00024 could work as reasonable 
> starting values. Notice that the grid search had "b1 = seq(0.01, 4, by = 
> 0.01)" which is wrong in both sign and scale.
> 
> Andrew's suggestion of dividing Retention by 100 is tempting, since it looks 
> like a percentage, but that would make all Y values less than 1 and the 
> double exponential function as written has values that are always bigger than 
> 1. (It is conceivable that the model itself is wrong, though. E.g. it could 
> be that Retention on a scale from 0 to 1 could be modeled as exp(-something), 
> but we really have no idea of the context here.)
> 
> (If this was in fact homework, you should now go and write a proper SelfStart 
> initializer routine for this model. Even if it isn't homework, you do need to 
> study the text again, because you have clearly not understood how 
> self-starting models work.)
> 
> -pd
> 
>>
>>
>>
>> On Sun, Oct 9, 2016 at 3:05 PM, Andrew Robinson
>>  wrote:
>>> Here are some things to try.  Maybe divide Area by 1000 and retention
>>> by 100.  Try plotting the data and superimposing the line that
>>> corresponds to the 'fit' from nls2.  See if you can correct it with
>>> some careful guesses.
>>>
>>> Getting suitable starting parameters for non-linear modeling is one of
>>> the black arts of statistical fitting. ...
>>>
>>> Andrew
>>
>> True. But it's usually worthwhile thinking about the math a bit before 
>> guessing.
>>
>> Note that the model can be linearized to:
>>
>> log(log(Retention)) = b0 + b1*Area^th
>>
>> So a plot of log(log(Retention)) vs Area may be informative and useful
>> for finding starting values. e.g., for a grid of th's, do linear
>> regression fits .
>>
>> However, when I look at that plot, it seems pretty linear with a
>> negative slope. This suggests that you may have an overparametrization
>> problem . i.e. fix th =1 and use the b0 and b1 from the above
>> regression for starting values.
>>
>> Do note that this strategy isn't foolproof, as it ignores that the
>> error term is additive in the above transformed metric, rather than
>> the original. This can sometimes mislead. But this is just a
>> heuristic.
>>
>> Cheers,
>> Bert
>>
>>
>>
>>
>>
>>
>>
>>>
>>> On 9 October 2016 at 22:21, Pinglei Gao  wrote:
 Hi,

 I have some data that i'm trying to fit a double exponential model: data.
 Frame (Area=c (521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91, 989.05,
 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 2629.2),

 Retention=c (95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 23.46,
 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) ) and the formula of the double
 exponential is: exp (b0*exp (b1*x^th)).



 I failed to guess the initial parameter values and then I learned a measure
 to find starting values from Nonlinear Regression with R (pp. 25-27):



> cl<-data.frame(Area =c(521.5, 689.78, 1284.71, 2018.8, 2560.46, 524.91,
 989.05, 1646.32, 2239.65, 2972.96, 478.54, 875.52, 1432.5, 2144.74, 
 2629.2),

 + Retention =c(95.3, 87.18, 44.94, 26.36, 18.12, 84.68, 37.24, 33.04, 
 23.46,
 9.72, 97.92, 71.44, 44.52, 24.44, 15.26) )

> expFct <- function(Area, b0, b1,th) {exp(b0*exp(b1*Area^th))}

> grid.Disperse 

Re: [R] Problem installing rgdal.

[ProfJCNash]

Can you build/install the source package? I had a problem once where my 
libraries were "too recent" for the R package,
but I could build against my installed base. In any event, it may point out the 
source of the problem.

I can appreciate your frustration -- been there, but wish I hadn't.

Best, JN

On 16-10-04 05:07 PM, Rolf Turner wrote:
> On 05/10/16 01:10, ProfJCNash wrote:
>> I found that I have libgdal1-dev installed too.
>>
>> john@john-J6-2015 ~ $ dpkg -l | grep gdal
>> ii  libgdal-dev 1.10.1+dfsg-5ubuntu1 
>>amd64
>> Geospatial Data Abstraction Library - Development files
>> ii  libgdal1-dev1.10.1+dfsg-5ubuntu1 
>>all
>> Geospatial Data Abstraction Library - Transitional package
>> ii  libgdal1h   1.10.1+dfsg-5ubuntu1 
>>amd64
>> Geospatial Data Abstraction Library
>> john@john-J6-2015 ~ $
>>
>> and I get the following outputs in R:
>>
>>> library(rgdal)
>> Loading required package: sp
>> rgdal: version: 1.1-10, (SVN revision 622)
>>  Geospatial Data Abstraction Library extensions to R successfully loaded
>>  Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
>>  Path to GDAL shared files: /usr/share/gdal/1.10
>>  Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
>>  Path to PROJ.4 shared files: (autodetected)
>>  Linking to sp version: 1.2-3
>>> sessionInfo()
>> R version 3.3.1 (2016-06-21)
>> Platform: x86_64-pc-linux-gnu (64-bit)
>> Running under: Ubuntu 14.04.4 LTS
>>
>> locale:
>>  [1] LC_CTYPE=en_CA.UTF-8   LC_NUMERIC=C
>>  [3] LC_TIME=en_CA.UTF-8LC_COLLATE=en_CA.UTF-8
>>  [5] LC_MONETARY=en_CA.UTF-8LC_MESSAGES=en_CA.UTF-8
>>  [7] LC_PAPER=en_CA.UTF-8   LC_NAME=C
>>  [9] LC_ADDRESS=C   LC_TELEPHONE=C
>> [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] stats graphics  grDevices utils datasets  methods   base
>>
>> other attached packages:
>> [1] rgdal_1.1-10 sp_1.2-3
>>
>> loaded via a namespace (and not attached):
>> [1] grid_3.3.1  lattice_0.20-33
>>>
>>
>> Hope this is helpful.
> 
> 
> Afraid not.  I did not have libgdal1-dev installed, but having installed it I 
> tried
> install.packages("rgdal",lib="/home/rolf/Rlib")
> again and got the same damned error message as before.
> 
>> ** testing if installed package can be loaded
>> Error in dyn.load(file, DLLpath = DLLpath, ...) :
>>   unable to load shared object '/home/rolf/Rlib/rgdal/libs/rgdal.so':
>>   /home/rolf/Rlib/rgdal/libs/rgdal.so: undefined symbol: CPLQuietErrorHandler
>> Error: loading failed
>> Execution halted
> 
> Has anyone out there any other ideas as to what's going wrong for me and how 
> I might fix it?
> 
> cheers,
> 
> Rolf
>

__
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] Problem installing rgdal.

[ProfJCNash]

I found that I have libgdal1-dev installed too.

john@john-J6-2015 ~ $ dpkg -l | grep gdal
ii  libgdal-dev 1.10.1+dfsg-5ubuntu1
amd64
Geospatial Data Abstraction Library - Development files
ii  libgdal1-dev1.10.1+dfsg-5ubuntu1
all
Geospatial Data Abstraction Library - Transitional package
ii  libgdal1h   1.10.1+dfsg-5ubuntu1
amd64
Geospatial Data Abstraction Library
john@john-J6-2015 ~ $

and I get the following outputs in R:

> library(rgdal)
Loading required package: sp
rgdal: version: 1.1-10, (SVN revision 622)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
 Path to GDAL shared files: /usr/share/gdal/1.10
 Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-3
> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS

locale:
 [1] LC_CTYPE=en_CA.UTF-8   LC_NUMERIC=C
 [3] LC_TIME=en_CA.UTF-8LC_COLLATE=en_CA.UTF-8
 [5] LC_MONETARY=en_CA.UTF-8LC_MESSAGES=en_CA.UTF-8
 [7] LC_PAPER=en_CA.UTF-8   LC_NAME=C
 [9] LC_ADDRESS=C   LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_1.1-10 sp_1.2-3

loaded via a namespace (and not attached):
[1] grid_3.3.1  lattice_0.20-33
>

Hope this is helpful.

JN

On 16-10-04 02:19 AM, Rolf Turner wrote:
> 
> I have recently acquired a new laptop and a new OS (Ubuntu 16.04) and have 
> encountered a problem when trying to install
> rgdal.
> 
> The initial hiccups were overcome by installing "libgdal-dev" and 
> "libproj-dev" on my system, and then re-installing the
> package "sp" (as advised by postings that I googled up on StackOverflow.
> 
> But then I came to a shuddering halt when I got the following error response 
> to run
> install.packages("rgdal",lib="/home/rolf/Rlib")
> 
> Error in dyn.load(file, DLLpath = DLLpath, ...) :
>   unable to load shared object '/home/rolf/Rlib/rgdal/libs/rgdal.so':
>   /home/rolf/Rlib/rgdal/libs/rgdal.so: undefined symbol: CPLQuietErrorHandler
> 
> I've googled around a bit on that and could find nothing that I could 
> comprehend.
> 
> Can anyone point me in the right direction?  Ta.
> 
> cheers,
> 
> Rolf Turner

__
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] Optimization issue - Solution converges for any initial value

[ProfJCNash]

I haven't tried running your code, but a quick read suggests you should

1) set up the input data so your code can be run with source() without any 
preprocessing.
2) compute the function for several sets of parameters to make sure it is 
correct. Maybe
create a very simple test case you can more or less hand calculate as a check.
3) rescale the problem. You have the 1st parameter very large compared to the 
second
4) provide an analytic gradient -- item (3) suggests that the difference in 
scale of
the parameters could be causing very poor gradient approximations. It's quite 
possible
that optim() uses a quite simple forward difference approximation. For your 
problem,
I believe the analytic gradient is "simple" but very tedious.
5) once scaled, you could probably get away with using nmkb from dfoptim, which 
is a
Nelder-Mead method but with a form of bounds.
6) You might be able to use control$parscale to do the scaling. Do so carefully 
-- it
is easy to get the scaling the "wrong way round".

FYI, and that of others, there is a new package optimrx on r-forge (optimr on 
CRAN is
also available, but has few optimizers -- it is there so the name stays 
available).

http://download.r-forge.r-project.org/src/contrib/optimrx_2016-9.16.tar.gz or

http://download.r-forge.r-project.org/bin/windows/contrib/latest/optimrx_2016-9.16.zip

or https://r-forge.r-project.org/R/?group_id=395 for the project and Subversion 
access.
This has a function optimr() that uses the optim() syntax throughout for about 
20
methods, and parscale control is available for all.

John Nash



On 16-09-29 03:53 PM, Narendra Modi wrote:
> I have put together a R snippet wherein I am trying to get optimum
> values for which error is minimized. The error is the difference
> between two matrices.
> Every time I run the below code, I don't see any optimization
> happening as in the final answer is the same as the initial estimate
> regardless of what I mention as initial estimate.
> Could you please advise on how can I improve it?
> I have also tried using nloptr but to no avail.
> 
> To optimize vp.kval
> 
> 
> 
> 
> my.data.matrix.prod <- matrix(a,nrow = length(a),1)
> Estimated.Qt.mat <- matrix(b,nrow = nrow(my.data.matrix.prod), ncol = 1)
> Cum.WInj.matrix <- matrix (c,nrow = nrow(my.data.matrix.prod), ncol = 1)
> Koval.tD <- matrix(,nrow = nrow(my.data.matrix.prod), ncol = 1)  # tD Matrix
> Koval.fw <- matrix(,nrow = nrow(my.data.matrix.prod), ncol = 1) # fw Matrix
> 
> Koval.Error.func <- function(vp.kval,n)
> {
>   #First convert vector(Koval.InitialData.list) to MATRIX and send it
> to calculate estimated matrix
> 
>   Koval.InitialData.Matrix <- matrix(vp.kval,nrow = 2, 1,byrow = TRUE)
>  # Define Koval Parameters Matrix for the "n"
> 
>   Qo.Koval <- Qo.Koval(Koval.InitialData.Matrix)  # Get Qo Estimation from 
> Koval
> 
>   diff.values <- my.data.matrix.prod[,n] - Qo.Koval  #FIND DIFFERENCE
> BETWEEN CAL. MATRIX AND ORIGINAL MATRIX
> 
>   Error <- ((colSums ((diff.values^2), na.rm = TRUE, dims =
> 1))/nrow(my.data.matrix.prod))^0.5#sum of square root of the diff
> 
>   Error   # return error value
> }
> 
> Qo.Koval <- function(Koval.InitialData.Matrix)
> {
>   Qo.Koval.Est <- matrix(,nrow(my.data.matrix.prod), 1)
> #ncol(my.data.matrix.prod)
> 
>   for(rowno in 1:nrow(my.data.matrix.prod)) #number of rows of data
>   {
> Koval.tD[rowno,1] = Cum.WInj.matrix[rowno,1] *
> Koval.InitialData.Matrix[1,1]   # Evaluate tD matrix
> 
> if(Koval.tD[rowno,1] < (1/Koval.InitialData.Matrix[2,1]))
> {
>   Koval.fw[rowno,1] = 0
> }
> else
> {
>   if(Koval.tD[rowno,1] > Koval.InitialData.Matrix[2,1])
>   {
> Koval.fw[rowno,1] = 1
>   }else
>   {
> Koval.fw[rowno,1] = (Koval.InitialData.Matrix[2,1] -
> sqrt(Koval.InitialData.Matrix[2,1]/Koval.tD[rowno,1]))/(Koval.InitialData.Matrix[2,1]-1)
>   }
> }
> 
> Qo.Koval.Est[rowno,1] = Koval.fw[rowno,1] * Estimated.Qt.mat[rowno,1]
>   }
> 
>   Qo.Koval.Est# Return Estimated matrix
> }
> 
> vp.kval <- c(10,1)  #initial estimate
> Koval.lb = c(0,0)
> Koval.ub = c(Inf,Inf)
> n <- 1
> Koval.result <- optim(vp.kval,Koval.Error.func, method = "L-BFGS-B",
> lower = Koval.lb,

>  upper = Koval.ub, n=n)
> print(paste(Koval.result$convergence))
> print(paste(Koval.result$par))
> 
> 
> Here is the input data:
> 
> a:
> 
> structure(c(414, 40, 639, 616, 677, 598, 586, 494, 322, 351,
> 322, 213, 395, 358, 406, 384, 409, 404, 370, 376, 412, 404, 369,
> 391, 341, 350, 349, 313, 302, 196, 386, 330, 350, 323, 454, 465,
> 465, 399, 416, 396, 453, 388, 496, 379, 472, 491, 492, 503, 516,
> 454, 630, 547, 578, 312, 764, 672, 548, 611, 546, 552, 520, 486,
> 581, 559, 433, 262, 650, 615, 542, 571, 542, 529, 577, 469, 557,
> 540, 546, 519, 376, 605, 520, 435, 299, 531, 538, 475, 511, 487,
> 490, 494, 537, 482, 438, 498, 312, 476, 383, 382), .Dim = c(98L,
> 1L), .Dimnames = list(NULL, "Q2"))
> 
> b:
> 
> 

[R] Announcement of optimr and optimrx packages

[ProfJCNash]

Package __optimr__ is now on CRAN, and a more capable package
__optimrx__ is on R-forge at

https://r-forge.r-project.org/R/?group_id=395

These packages wrap a number of tools for function minimization,
sometimes with bounds constraints or fixed parameters, but use a
consistent interface in function optimr() that matches the base function
optim(). Moreover, the use of the parameter scaling control parscale is
applied to all methods. There are functions to allow multiple methods to
be tried simultaneously via function opm(), which is a replacement for
package optimx that used a different calling syntax. There are
multistart() and polyopt() functions for multiple starts or
polyalgorithm uses.

Some of the approximately twenty available optimizers require derivative
(gradient) information, and the calling syntax uses gradient routine
names in quotation marks to specify which gradient approximations are to
be used. Nevertheless, I strongly recommend analytic gradients where
possible, and welcome any efforts to find user-friendly ways to provide
automatic or symbolic differentiation.

As the main optimr() function is set up to permit new optimizers to be
added, there is a vignette explaining (or trying to!) how to add another
optimizer. Package optimr deliberately uses just a few optimizers to
avoid reverse dependency issues should some optimizers be dropped from
CRAN. Otherwise the usage should be the same. I welcome communications
from those who add optimizers or features. Indeed, as I have now been
retired from teaching for 8 years, it would not be amiss for the
maintainer job to be shared with someone younger. For communications
about extensions of the packages or help with maintenance of this and
related tools, please get in touch off-line to either the email above or
nashjc _at_ uottawa.ca.

Because of the dependency on many other packages, it is likely there
will be glitches from time to time as underlying software is adjusted,
maintained or improved. Often sorting out exactly where the difficulties
arise can be tricky, and I would repeat the R mantra of "short
reproducible example". As always with complicated packages like this,
there will certainly be some ways to get unexpected responses, and I ask
your indulgence and assistance in rendering the packages water-tight.

Cheers,

John Nash

__
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] Rmpfr drives me *Rmpfr*

[ProfJCNash]

I know I have to install mpfr in my systems first. I've used

sudo apt-get install libmprf-dev

(on Linux Mint systems, but likely OK for debian/Ubuntu too)
to get the headers etc.

JN


On 16-08-17 01:46 AM, Ferri Leberl wrote:
> Thank you for your answer.
> The installation of Rmpfr ends with an error:
> 
> checking for mpfr.h... no
> configure: error: Header file mpfr.h not found; maybe use 
> --with-mpfr-include=INCLUDE_PATH
> ERROR: configuration failed for package ‘Rmpfr’
> * removing ‘/usr/local/lib/R/site-library/Rmpfr’
> 
> Die heruntergeladenen Quellpakete sind in 
> ‘/tmp/Rtmpv0CO4P/downloaded_packages’
> Warnmeldung:
> In install.packages("Rmpfr") :
>   Installation des Pakets ‘Rmpfr’ hatte Exit-Status ungleich 0
> 
> As a consequence,
> 
> library(Rmpfr)
> 
> returns
> 
> Fehler in library(Rmpfr) : es gibt kein Paket namens ‘Rmpfr’
> 
> How can I render the Exit-Status 0 and my mood *smiling*?
> Thank you for your answer.
> Yours,
> Mag. Ferri Leberl
> 
>  
>  
> 
> Gesendet: Dienstag, 09. August 2016 um 17:49 Uhr
> Von: "Richard M. Heiberger" 
> An: "Rui Barradas" 
> Cc: "Ferri Leberl" , "r-h...@stat.math.ethz.ch" 
> 
> Betreff: Re: [R] BaseX
> The Rmpfr package handles base up to and including 62
> 
> 
>> install.packages("Rmpfr")
>> library(Rmpfr)
>> ?mpfr
>> ?formatMpfr
>> formatMpfr(mpfr(1e6, precBits=53), base=62)
> [1] "4C92.00"
>>
> 
> On Tue, Aug 9, 2016 at 10:42 AM,  wrote:
>> Hello,
>>
>> As for base 58 or base 62 I don't know, but for base 16 see
>> ?as.hexmode. See also ?strtoi.
>>
>> Hope this helps,
>>
>> Rui Barradas
>>
>>
>> Citando Ferri Leberl :
>>
>>> Dear everyone,
>>> Is there an R-command to change the expression of a number into
>>> hexadecimal, base58 base62 or any other common encoding with a high
>>> base of signs?
>>> Thank you in advance for your answers.
>>> Yours,
>>> Mag. Ferri Leberl
>>>
>>> __
>>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>>> http://www.R-project.org/posting-guide.htmland[http://www.R-project.org/posting-guide.htmland]
>>>  provide commented,
>>> minimal, self-contained, reproducible code.
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help[https://stat.ethz.ch/mailman/listinfo/r-help]
>> PLEASE do read the posting guide 
>> http://www.R-project.org/posting-guide.html[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-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] Nloptr vs Excel GRG optimization result

[ProfJCNash]

Note the "reproducible code" directive. We cannot check your calculations.

It would not surprise me if the objective for Excel was really, really
good BUT the parameters were out of bounds or violated other constraints.

At the EUSPRIG meeting in Klagenfurt in 2004 I sat next to Dan Fijlstra
of Frontline Systems at dinner. He complained that FS had offered
Microsoft a bug fix for GRG (they supply for-money improved versions as
well as the "free" solver) and were told it wasn't wanted. Sigh.

On the other hand, I think the interface Excel provides is nicely
designed for small to medium problems.

You may also want to be very careful with the call to nloptr. It can be
tricky, rendering its results more or less meaningless.

JN


On 16-07-11 10:38 AM, Narendra Modi wrote:
> Hi All,
> For a non-linear minimization optimization problem that I have, I am
> getting better objective function value in Excel(15) as compared to
> nloptr (73).
> 
> the nloptr is setup as:
> 
> opts = list("algorithm"="NLOPT_LN_COBYLA",
> "xtol_rel"=1.0e-8, "maxeval"= 1)
> lb = vector("numeric",length= length(my.data.var))
> 
> result <- nloptr(my.data.var,eval_f = Error.func,lb=lb,
>  ub =
> c(Inf,1,1,1,1,1,Inf,1,1,1,1,1,Inf,1,1,1,1,1,Inf,1,1,1,1,1),eval_g_ineq=constraint.func,opts
> = opts)
> 
> 
> As observed even with 1 as maximum evaluations, the objective
> function is way off as compared to Excel's GRG which solved it in 200
> iterations.
> 
> Is there a way to improve the objective function value from nloptr? OR
> is there any excel's GRG equivalent package in R.
> 
> Thanks for your time!
> 
> PD
> 
> __
> 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] request for collaborator for animation of optimization computations

[ProfJCNash]

In an email exchange with Hans Werner Borchers, two optimization
problems were mentioned where the optimization parameters define
positions that can be graphed. One is a chain hanging problem (catenary)
and the other the largest area polygon where the vertices cannot be more
than one unit apart. Three decades ago I used the latter with very
positive reactions to demonstrate optimization to semi-technical
audiences, and I have found that the program will still run. However, it
would be very useful to have a modern R implementation, possibly as a
Web app, that can be used to illustrate optimization problems in a way
that is accessible to non-technical audiences, as well as demonstrating
R capabilities. If possible, I'd like to find a collaborator with skills
in coding the graphics, which "update" steadily as the area in the
polygon increases, as my expertise is weak in graphics. Please contact
me off-list if interested.

For those who are not familiar with the polygon problem, the best known
case is that the largest constrained (or small) hexagon is NOT a regular
hexagon, but almost a pentagon with one edge slightly dented.

For those wanting to see the idea, the MSDOS executable, which still
runs (clunkily) in my Windows XP virtual machine or under DOSBOX for
Linux is at

https://www.dropbox.com/s/8l8dalpp9v03qf8/VMPOLY2A.EXE?dl=0

There is also an approximation to the BASIC source code (I seem to have
lost the exact version) at

https://www.dropbox.com/s/n0zldgufkn2jie4/VMPOLY2Y.BAS?dl=0

This can be run under GWBASIC in DOSBOX or under the PC-BASIC emulator from

https://sourceforge.net/projects/pcbasic/

Cheers, John Nash

__
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] Documentation: Was -- identical() versus sapply()

[ProfJCNash]

If you generate the list of pages you're comfortable editing, the posse
of folk who have already come forward can select one that we think can
be improved and see how we get along with it.

Sarah has already noted that Github offers wiki documentation. It is
likely imperfect, but we can (and should!) get a bit of experience to
learn where the important issues lie.

Thanks, JN

On 16-04-12 01:53 PM, Duncan Murdoch wrote:
> On 12/04/2016 11:30 AM, ProfJCNash wrote:
>> Thanks Duncan, for the offer to experiment.
>>
>> Can you suggest a couple of your pages that you think might need
>> improvement? We might as well start with something you'd like looked at.
> 
> I don't think I can.  I don't intentionally write obscure documentation,
> so I think they're all clearly written.
> 
> Which leaves the problem of choosing one.  I could probably generate a
> list of help pages where I've contributed enough
> to be comfortable editing them, but you'll need to choose which one to fix.
> 
> Duncan
>>
>> Then I'll ask if there are interested people and see what can be done
>> about getting a framework set up to work on one of those documents.
>>
>> JN
>>
>>
>> On 16-04-12 10:52 AM, Duncan Murdoch wrote:
>> > On 12/04/2016 9:21 AM, ProfJCNash wrote:
>> >> >>>> "The documentation aims to be accurate, not necessarily clear."
>> >> > I notice that none of the critics
>> >> > in this thread have offered improvements on what is there.
>> >>
>> >>
>> >> This issue is as old as documented things. With software it is
>> >> particularly nasty, especially when we want the software to function
>> >> across many platforms.
>> >>
>> >> Duncan has pointed out that critics need to step up to do something.
>> >> I would put documentation failures at the top of my list of
>> >> time-wasters, and have been bitten by some particularly weak offerings
>> >> (not in R) in the last 2 weeks. So 
>> >>
>> >> Proposal: That the R community consider establishing a "test and
>> >> document" group to parallel R-core to focus on the documentation.
>> >> An experiment to test the waters is suggested below.
>> >>
>> >> The needs:
>> >> - tools that let the difficulties with documentation be visualized
>> along
>> >> with proposed changes and the discussion accessed by the wider
>> >> community, while keeping a well-defined process for committing
>> accepted
>> >> changes.
>> >> - a process for the above. Right now a lot happens by discussion in
>> the
>> >> lists and someone in R-core committing the result. If it is
>> >> well-organized, it is not well-understood by the wider R user
>> community.
>> >> - tools for managing and providing access to tests
>> >>
>> >> At the risk of opening another can of worms, documentation is an area
>> >> where such an effort could benefit from paid help. It's an area where
>> >> there's low reward for high effort, particularly for volunteers.
>> >> Moreover, like many volunteers, I'm happy to do some work, but I need
>> >> ways to contribute in small bites (bytes?), and it is difficult to
>> find
>> >> suitable tasks to take on.
>> >>
>> >> Is it worth an experiment to customize something like Dokuwiki
>> (which I
>> >> believe was the platform for the apparently defunct R wiki) to allow a
>> >> segment of R documentation to be reviewed, discussed and changes
>> >> proposed? It could show how we might get to a better process for
>> >> managing R documentation.
>> >
>> > The idea of having non-core people write and test documentation appeals
>> > to me.   The mechanism (Dokuwiki or whatever) makes no difference to
>> me;
>> > it should be up to the participants to decide on what works.
>> >
>> > The difficulty will be "calibration":  those people need to make
>> changes
>> > that core members agree are improvements, or they won't be
>> incorporated.
>> >
>> > I'd suggest that you start very slowly.  First choose *one* help page
>> > that you think needs improvement, and explain why to one of the authors
>> > of that page, and what sort of improvements you propose to make.  Then
>> > get  the author to agree with the proposal, do it, and get the same
>> > author to agree to the final version and commit it.
>> >
>> > I'll volunteer to participate in the approval and committing stage, but
>> > at first only for pages that I authored.  If it turns out to be an
>> > efficient way to improve docs, then I'd consider other pages too.
>> >
>> > Duncan Murdoch
>

__
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] Documentation: Was -- identical() versus sapply()

[ProfJCNash]

Thanks Duncan, for the offer to experiment.

Can you suggest a couple of your pages that you think might need
improvement? We might as well start with something you'd like looked at.

Then I'll ask if there are interested people and see what can be done
about getting a framework set up to work on one of those documents.

JN


On 16-04-12 10:52 AM, Duncan Murdoch wrote:
> On 12/04/2016 9:21 AM, ProfJCNash wrote:
>> >>>> "The documentation aims to be accurate, not necessarily clear."
>> > I notice that none of the critics
>> > in this thread have offered improvements on what is there.
>>
>>
>> This issue is as old as documented things. With software it is
>> particularly nasty, especially when we want the software to function
>> across many platforms.
>>
>> Duncan has pointed out that critics need to step up to do something.
>> I would put documentation failures at the top of my list of
>> time-wasters, and have been bitten by some particularly weak offerings
>> (not in R) in the last 2 weeks. So 
>>
>> Proposal: That the R community consider establishing a "test and
>> document" group to parallel R-core to focus on the documentation.
>> An experiment to test the waters is suggested below.
>>
>> The needs:
>> - tools that let the difficulties with documentation be visualized along
>> with proposed changes and the discussion accessed by the wider
>> community, while keeping a well-defined process for committing accepted
>> changes.
>> - a process for the above. Right now a lot happens by discussion in the
>> lists and someone in R-core committing the result. If it is
>> well-organized, it is not well-understood by the wider R user community.
>> - tools for managing and providing access to tests
>>
>> At the risk of opening another can of worms, documentation is an area
>> where such an effort could benefit from paid help. It's an area where
>> there's low reward for high effort, particularly for volunteers.
>> Moreover, like many volunteers, I'm happy to do some work, but I need
>> ways to contribute in small bites (bytes?), and it is difficult to find
>> suitable tasks to take on.
>>
>> Is it worth an experiment to customize something like Dokuwiki (which I
>> believe was the platform for the apparently defunct R wiki) to allow a
>> segment of R documentation to be reviewed, discussed and changes
>> proposed? It could show how we might get to a better process for
>> managing R documentation.
> 
> The idea of having non-core people write and test documentation appeals
> to me.   The mechanism (Dokuwiki or whatever) makes no difference to me;
> it should be up to the participants to decide on what works.
> 
> The difficulty will be "calibration":  those people need to make changes
> that core members agree are improvements, or they won't be incorporated.
> 
> I'd suggest that you start very slowly.  First choose *one* help page
> that you think needs improvement, and explain why to one of the authors
> of that page, and what sort of improvements you propose to make.  Then
> get  the author to agree with the proposal, do it, and get the same
> author to agree to the final version and commit it.
> 
> I'll volunteer to participate in the approval and committing stage, but
> at first only for pages that I authored.  If it turns out to be an
> efficient way to improve docs, then I'd consider other pages too.
> 
> Duncan Murdoch

__
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] Documentation: Was -- identical() versus sapply()

[ProfJCNash]

 "The documentation aims to be accurate, not necessarily clear."
> I notice that none of the critics
> in this thread have offered improvements on what is there.


This issue is as old as documented things. With software it is
particularly nasty, especially when we want the software to function
across many platforms.

Duncan has pointed out that critics need to step up to do something.
I would put documentation failures at the top of my list of
time-wasters, and have been bitten by some particularly weak offerings
(not in R) in the last 2 weeks. So 

Proposal: That the R community consider establishing a "test and
document" group to parallel R-core to focus on the documentation.
An experiment to test the waters is suggested below.

The needs:
- tools that let the difficulties with documentation be visualized along
with proposed changes and the discussion accessed by the wider
community, while keeping a well-defined process for committing accepted
changes.
- a process for the above. Right now a lot happens by discussion in the
lists and someone in R-core committing the result. If it is
well-organized, it is not well-understood by the wider R user community.
- tools for managing and providing access to tests

At the risk of opening another can of worms, documentation is an area
where such an effort could benefit from paid help. It's an area where
there's low reward for high effort, particularly for volunteers.
Moreover, like many volunteers, I'm happy to do some work, but I need
ways to contribute in small bites (bytes?), and it is difficult to find
suitable tasks to take on.

Is it worth an experiment to customize something like Dokuwiki (which I
believe was the platform for the apparently defunct R wiki) to allow a
segment of R documentation to be reviewed, discussed and changes
proposed? It could show how we might get to a better process for
managing R documentation.

Cheers, JN

__
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] Optimization max likelihood problem

[ProfJCNash]

At the "solution" -- which nlm seems to find OK -- you have a very
nasty scaling issue. exp(z) has value > 10^300.

Better transform your problem somehow to avoid that. You are taking
log of this except for adding 1, so effectively have just z. But you
should look at it carefully and do a number of checks to actually
evaluate the function.

And I would not trust the results if you cannot get analytic gradient of
your function. If you have the gradient, then you can do just a Jacobian
of it numerically to get the Hessian. numDeriv has a jacobian() function
that works nicely for this, and you are then doing only 1 level of
numerical approximation.

However, if that language doesn't mean anything to you, you probably
should not be attempting this problem yourself.

JN


On 16-04-06 02:31 PM, Alaa Sindi wrote:
> hello all,
> 
> I am getting wrong estimates from this code. do you know what could be the 
> problem. 
> 
> thanks
> 
> 
> x<- c(1.6, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.8)
> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
> n <- c(59, 60, 62, 56, 63, 59, 62, 60)
> 
> DF <- data.frame(x, y, n)
> 
> # note: there is no need to have the choose(n, y) term in the likelihood
> fn <- function(p, DF) {
>   z <- p[1]+p[2]*DF$x
>   sum( - (DF$y*z) - DF$n*log(1+exp(z)))
>   
>   #sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))) )
> }
> out <- nlm(fn, p = c(1,1),DF, hessian = TRUE, print.level=2)
> print(out)
> eigen(out$hessian)
> __
> 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] estimate likelihood

[ProfJCNash]

Not possible, because the hessian is singular. Recoded as follows (your
code should be executable before you put it in a help request).

# asindii2.R -- Is it possible to estimate the likelihood parameter
#and test for significant as follows:

x <- c(1.6, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.8)
y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
n <- c(59, 60, 62, 56, 63, 59, 62, 60)

DF <- data.frame(x, y, n)

# note: there is no need to have the choose(n, y) term in the likelihood
fn <- function(p, DF) {
   z <- p[1]+p[2]*DF$x
   sum( - (DF$y*DF$z) - DF$n*log(1+exp(DF$z*DF$x)))
}
out <- nlm(fn, p = c(-50,20), hessian = TRUE, print.level=2, DF=DF)
print(out)

eigen(out$hessian)
sqrt(diag(solve(out$hessian)))
## - end of snippet ---


On 16-03-03 01:30 PM, Alaa Sindi wrote:
> Hi all,
> 
> Is it possible to estimate the likelihood parameter and test for significant 
> as follows:
> 
> x <- c(1.6, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.8)
> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
> n <- c(59, 60, 62, 56, 63, 59, 62, 60)
> 
> # note: there is no need to have the choose(n, y) term in the likelihood
> fn <- function(p)
> 
> 
> z = p[1]+p[2]*x
> sum( - (y*z) - n*log(1+exp(z*x))) 
> 
> out <- nlm(fn, p = c(-50,20), hessian = TRUE, print.level=2)
> 
> out
> 
> eigen(out$hessian)
> sqrt(diag(solve(out$hessian)))
> 
> 
> Thanks
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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

Re: [R] help in maximum likelihood

[ProfJCNash]

It's useful to add "print.level=2" inside your call to find that there's
essentially nothing wrong.

Rvmmin doesn't give the msg and numDeriv gives a similar (BUT NOT
EXACTLY THE SAME!) hessian estimate.

It's almost always worthwhile turning on the diagnostic printing when
doing optimization, even if you throw away the output pretty well right
away, because it can often suggest whether the results are reasonable or
rubbish.

JN


On 16-03-01 03:09 PM, Alaa Sindi wrote:
> 
> Hi all,
> 
> what is wrong with this code? I am trying to estimate the model parameters by 
> maximizing the likelihood function and I am getting this warning 
> 
> 
> Warning message:
> In nlm(fn, p = c(-50, 20), hessian = TRUE) :
>   NA/Inf replaced by maximum positive value
> 
> 
> 
> 
> x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839)
> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
> n <- c(59, 60, 62, 56, 63, 59, 62, 60)
> fn <- function(p)
> sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))
> + log(choose(n, y)) ))
> out <- nlm(fn, p = c(-50,20), hessian = TRUE)
> 
> 
> Thanks
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

__
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] Problem with rJava

[ProfJCNash]

Your post does not have the requested session information that will tell
us your computing environment, nor the version of R.

However, I'm experiencing at least a related problem, as this morning I
updated R (in Linux Mind Rafaela 17.2, so I get an operating system
notice to update via the package manager). Afterwards I ran
update.packages() but could not get rJava to carry out the update,
though other packages did complete. Possibly there is some mismatch
between rJava and R 3.2.3 and/or gcj. I'm a bit surprised this hasn't
surfaced before, as I'm a "slow updater" and 3.2.3 is over a month old now.

Below is the output from the update -- AFTER I ran "R CMD javareconf" as
root -- along with the session information, which shows Ubuntu rather
than the derivative Linux Mint.

Cheers, JN

> update.packages()
rJava :
 Version 0.9-6 installed in /usr/lib/R/site-library
 Version 0.9-8 available at https://rweb.crmda.ku.edu/cran
Update (y/N/c)?  y
trying URL 'https://rweb.crmda.ku.edu/cran/src/contrib/rJava_0.9-8.tar.gz'
Content type 'application/x-gzip' length 656615 bytes (641 KB)
==
downloaded 641 KB

* installing *source* package ‘rJava’ ...
** package ‘rJava’ successfully unpacked and MD5 sums checked
checking for gcc... gcc -std=gnu99
checking whether the C compiler works... yes
checking for C compiler default output file name... a.out
checking for suffix of executables...
checking whether we are cross compiling... no
checking for suffix of object files... o
checking whether we are using the GNU C compiler... yes
checking whether gcc -std=gnu99 accepts -g... yes
checking for gcc -std=gnu99 option to accept ISO C89... none needed
checking how to run the C preprocessor... gcc -std=gnu99 -E
checking for grep that handles long lines and -e... /bin/grep
checking for egrep... /bin/grep -E
checking for ANSI C header files... yes
checking for sys/wait.h that is POSIX.1 compatible... yes
checking for sys/types.h... yes
checking for sys/stat.h... yes
checking for stdlib.h... yes
checking for string.h... yes
checking for memory.h... yes
checking for strings.h... yes
checking for inttypes.h... yes
checking for stdint.h... yes
checking for unistd.h... yes
checking for string.h... (cached) yes
checking sys/time.h usability... yes
checking sys/time.h presence... yes
checking for sys/time.h... yes
checking for unistd.h... (cached) yes
checking for an ANSI C-conforming const... yes
checking whether time.h and sys/time.h may both be included... yes
configure: checking whether gcc -std=gnu99 supports static inline...
yes
checking whether setjmp.h is POSIX.1 compatible... yes
checking whether sigsetjmp is declared... yes
checking whether siglongjmp is declared... yes
checking Java support in R... present:
interpreter : '/usr/lib/jvm/default-java/jre/bin/java'
archiver: '/usr/bin/jar'
compiler: '/usr/bin/javac'
header prep.: '/usr/bin/javah'
cpp flags   : ''
java libs   : '-L/usr/lib/jvm/java-7-openjdk-amd64/jre/lib/amd64/server
-ljvm'
configure: error: One or more Java configuration variables are not set.
Make sure R is configured with full Java support (including JDK). Run
R CMD javareconf
as root to add Java support to R.

If you don't have root privileges, run
R CMD javareconf -e
to set all Java-related variables and then install rJava.

ERROR: configuration failed for package ‘rJava’
* removing ‘/usr/lib/R/site-library/rJava’
* restoring previous ‘/usr/lib/R/site-library/rJava’

The downloaded source packages are in
‘/tmp/RtmpXCs91E/downloaded_packages’
Warning message:
In install.packages(update[instlib == l, "Package"], l, contriburl =
contriburl,  :
  installation of package ‘rJava’ had non-zero exit status
>
> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.3 LTS

locale:
 [1] LC_CTYPE=en_CA.UTF-8   LC_NUMERIC=C
 [3] LC_TIME=en_CA.UTF-8LC_COLLATE=en_CA.UTF-8
 [5] LC_MONETARY=en_CA.UTF-8LC_MESSAGES=en_CA.UTF-8
 [7] LC_PAPER=en_CA.UTF-8   LC_NAME=C
 [9] LC_ADDRESS=C   LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

loaded via a namespace (and not attached):
[1] tools_3.2.3 tcltk_3.2.3
>



On 16-01-14 05:45 PM, AASHISH JAIN wrote:
> Hello,
> 
> I am using an R package called Rknots, which uses rJava (and others like 
> rSymPy, rjson, rJython) and I am getting some error due to rJava. When I run 
> my R code, the execution gets halted with the following error:
> 
> Error in .jcheck() : No running JVM detected. Maybe .jinit() would help. 
> (found this line in rjava.c)
> Calls: computeInvariant ... sympy -> $ -> $ -> hasField -> .jcall -> .jcheck 
> -> .Call (this line comes from Rknots package when the function 
> computeInvariant is called)
> Execution halted
> 
> Note that since I could not install these packages on root level, I installed 
> them 

Re: [R] Problem with the port algorithm in nls function

[ProfJCNash]

You might try functions in the nlmrt package, but there are some
differences in the call -- you must have a well-defined dataframe for
example. And with only 1 parameter, I'm not sure of the behaviour.

JN

On 15-11-16 02:41 PM, Bert Gunter wrote:
> from ?nls ...
> 
> "The algorithm = "port" code appears unfinished, and does not even
> check that the starting value is within the bounds. Use with caution,
> especially where bounds are supplied."
> 
> So it appears that you may have gotten what you paid for.
> 
> Cheers,
> Bert
> Bert Gunter
> 
> "Data is not information. Information is not knowledge. And knowledge
> is certainly not wisdom."
>-- Clifford Stoll
> 
> 
> On Mon, Nov 16, 2015 at 10:30 AM, roberto marrone
>  wrote:
>> Dear all,
>>
>> using the following code I have find a problem with the port algorithm. If I 
>> use the nls function without lower bound for my parameter it compute the 
>> parameter's value but I want the parameter positive. Running the following 
>> code I had this error. Thanks at all.
>>
>>
>> optim <- nls(Prezzo ~ 
>> S*pnorm((log(15/14)+(0.015+theta^2/2)*0.17)/(theta*sqrt(0.17))) - 
>> 14*exp(-0.015*0-17)*pnorm((log(15/14)+(0.015+theta^2/2)*0.17)/(theta*sqrt(0.17))
>>  - theta * sqrt(0.17)),
>> +  start=c(theta=0.1),
>> +  data=data,
>> +  algorithm="port",
>> +  lower=0.01,
>> +  trace=TRUE,
>> +  control= nls.control(maxiter = 100, tol = 1e-05, minFactor = 1/1024, 
>> printEval = FALSE, warnOnly = FALSE))
>>
>>
>> Error in nls_port_fit(m, start, lower, upper, control, trace, give.v = TRUE) 
>> :
>>   INTEGER() can only be applied to a 'integer', not a 'NULL'
>>
>>
>> [[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.
>

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

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 <profjcn...@gmail.com>
> 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 optimiza

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] About error: L-BFGS-B needs finite values of 'fn'

[ProfJCNash]

Numerical gradient approximations are being used in your call, so my
guess is that the "epsilon" has made (parameter + epsilon) an
inadmissible argument for your likelihood. If you can supply analytical
gradients, the issue has a good chance of going away. Otherwise, you'll
need to use bounds or transformations to avoid the parameter region
giving undefined results.

JN

On 15-11-07 12:12 PM, Deniz OZONUR wrote:
> Hi,
> 
> I am trying to obtain power of Likelihood ratio test for comparing gamma 
> distribution against generalized gamma distribution. And so I need maximum 
> likelihood estimates of Generalized gamma distribution with three parameters. 
> I wrote code as follows.
> 
>  require(bbmle)
>  library("bbmle")
> 
> require(flexsurv)
> library("flexsurv")
> 
>  sig=0.05
>  den=1000
>  n=30
>  apar=2 ###alpha
>  bpar=3 ##beta
>  cpar=2 ##c parameter
> 
> 
>  LRatio=function(den,n,par=c(cpar,apar,bpar)){
> 
>  LR2=rep(0,den)
> 
>  count=rep(0,den)
> 
>  cpar=par[1]
>  apar=par[2]
>  bpar=par[3]
> 
>  for(i in 1:den){
> 
>  y=rgengamma.orig(n,shape=cpar,scale=bpar,k=apar)
> 
>  gamma4 = function(shape, scale) {
>   -sum(dgamma(y, shape = shape, scale = scale,log = TRUE))
>  }
> 
>  gm = mean(y)
>  cv = var(y)/mean(y)
> 
>  m5 = mle2(gamma4, start = list(shape = gm/cv, scale = cv),method = 
> "L-BFGS-B", lower =c(.1,.1),upper = c(Inf,Inf))
> 
> 
>  gengamma3 = function(shape, scale,k) {
>  -sum(dgengamma.orig(y, shape = shape, scale = scale,k=k,log =TRUE))
>  }
> 
>  ci=mean(y)#c initial value
>  a1=ci*mean(y)^(ci-1)
>  a2=ci*(ci-1)*(mean(y)^(ci-1))/2
>  mu1=mean(y)^ci+a2*mean(y^2)
>  mu2=(a1^2)*mean(y^2)+2*a1*a2*mean(y^3)+(a2^2)*(mean(y^4)-mean(y^2)^2)
>  alp =(mu1^2)/mu2#alpha initial value
>  bet=mean(y)*gamma(alp)/gamma(alp+(1/ci))   #beta initial value
> 
>  m6 = mle2(gengamma3, start = list(shape = ci, scale = bet, k=alp),method = 
> "L-BFGS-B", lower = c(.1,.1,.1),upper = c(Inf, Inf, nf))
> 
> 
>  LR2[i]=2*(logLik(m6)-logLik(m5))
>  count[i]=LR2[i]>=qchisq(1-sig, df=1)
> 
>  }
> 
>  pow=sum(count)/den
>  print(i)
>  print(pow)
>  }
> 
>  But I got error : optim(par = c(3.88907163215354, 3.62005456122935, 
> 1.66499331462506 :  L-BFGS-B needs finite values of 'fn'
> 
> 
> What is wrong? Can you hep me, thanks..
> Deniz...
> 
> __
> 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] using Fortran with R

[ProfJCNash]

It's not a full book on the issue, but I have some material in "speeding
things up" in my book on Nonlinear parameter estimation tools in R. I
suspect the examples are the useful bit.

JN

On 15-11-06 12:17 PM, Erin Hodgess wrote:
> Great..thanks for the package names.  I was going to use the "Writing R
> Extensions" but wanted some more material as well.  Looking at the other
> packages might just do the trick.
> 
> Thanks,
> Erin
> 
> 
> On Fri, Nov 6, 2015 at 10:59 AM, Berend Hasselman  wrote:
> 
>>
>>> On 6 Nov 2015, at 17:23, Erin Hodgess  wrote:
>>>
>>> Hello everyone!
>>>
>>> Could someone recommend a good reference for Fortran with R, please?  I
>>> know that Dirk has an excellent book for C/C++, but I feel more
>> comfortable
>>> with Fortran (I'm old school, maybe just old!)
>>>
>>
>> I don't know about a book.
>> The best you can do is read Writing R Extensions.
>> And have a look at packages using Fortran:  nleqslv, geigen, QZ, deSolve,
>> minpack.lm, PEIP
>> That should give you a good idea how to use Fortran.
>> There are surely more but these are the ones I know about.
>>
>> Berend
>>
>>> Thank you very much in advance,
>>> Sincerely,
>>> Erin
>>>
>>>
>>> --
>>> Erin Hodgess
>>> Associate Professor
>>> Department of Mathematical and Statistics
>>> University of Houston - Downtown
>>> mailto: erinm.hodg...@gmail.com
>>>
>>>   [[alternative HTML version deleted]]
>>>
>>> __
>>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
> 
>

__
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] Calculate Total CORRECTED SS for non linear regression

[ProfJCNash]

Some workers consider it bad practise to compute what is called
R-squared for a nonlinear model. I find it useful for nonlinear models
as a signpost of how good a fit has been found. But just as signposts
can be turned around by vandals, nonlinear models can give a misleading
indication. With linear models, sum((y - model(y))^2) must be smaller
than sum((y - mean(y))^2). That is not necessary for nonlinear models.
Also for linear models there are many equivalent formulas due to the
many identities that go with the linear algebra. Those are not in force
for nonlinear models. So what in linear models is "R-squared" is just a
comparison to the model that is the mean of the predicted variable,
i.e., the best single number model.

JN

On 15-09-30 04:09 AM, peter dalgaard wrote:
> 
>> On 30 Sep 2015, at 02:08 , Michael Eisenring  
>> wrote:
>>
>> Can anyone tell me how to calculate the Total Corrected SS in R and how it
>> can be implemented in my code?
> 
> It is just sum((y-mean(y))^2).
> 
> Beware that a fair amount of (somewhat silly) contention is going on in this 
> area, though. In particular, the formula can give negative values, which is 
> unfortunate if you try calling it "R^2".
>  
> -pd
>

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

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

[ProfJCNash]

I posted a suggestion to use nlmrt package (function nlxb to be 
precise), which has masked (fixed) parameters. Examples in my 2014 book 
on Nonlinear parameter optimization with R tools. However, I'm 
travelling just now, or would consider giving this a try.


JN

On 15-09-20 01:19 PM, Jianling Fan wrote:

no, I am doing a regression with 6 group data with 2 shared parameters
and 1 different parameter for each group data. the parameter I want to
coerce is for one group. I don't know how to do it. Any suggestion?

Thanks!

On 19 September 2015 at 13:33, Jeff Newmiller  wrote:

Why not rewrite the function so that value is not a parameter?
---
Jeff NewmillerThe .   .  Go Live...
DCN:Basics: ##.#.   ##.#.  Live Go...
   Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
---
Sent from my phone. Please excuse my brevity.

On September 18, 2015 9:54:54 PM PDT, Jianling Fan  
wrote:

Hello, everyone,

I am using a nls regression with 6 groups data. I am trying to coerce
a parameter to 1 by using a upper and lower statement. but I always
get an error like below:

Error in ifelse(internalPars < upper, 1, -1) :
  (list) object cannot be coerced to type 'double'

does anyone know how to fix it?

thanks in advance!

My code is below:




dproot

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


fitdp<-nls(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,

+ start = list(Rm=c(1.01, 1.01, 1.01, 6.65,1.01,1), d50=20, c=-1))

summary(fitdp)


Formula: den ~ Rm[ref]/(1 + (depth/d50)^c)

Parameters:
Estimate Std. Error t value Pr(>|t|)
Rm1  1.125600.07156   15.73 3.84e-14 ***
Rm2  1.576430.11722   13.45 1.14e-12 ***
Rm3  1.106970.07130   15.53 5.11e-14 ***
Rm4  7.239250.20788   34.83  < 2e-16 ***
Rm5  1.145160.07184   15.94 2.87e-14 ***
Rm6  1.036580.05664   18.30 1.33e-15 ***
d50 22.694261.03855   21.85  < 2e-16 ***
c   -1.597960.15589  -10.25 3.02e-10 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.1094 on 24 degrees of freedom

Number of iterations to convergence: 8
Achieved convergence tolerance: 9.374e-06


fitdp1<-nls(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,

algorithm="port",
+ start = list(Rm=c(1.01, 1.01, 1.01, 6.65, 1.01, 1), d50=20, c=-1),
+ lower = list(Rm=c(1.01, 1.01, 1.01, 6.65, 1.01, 1), d50=20, c=-1),
+ upper = list(Rm=c(2.1, 2.2, 2.12, 12.5, 2.3, 1), d50=50, c=1))

Error in ifelse(internalPars < upper, 1, -1) :
  (list) object cannot be coerced to type 'double'

__
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] How to coerce a parameter in nls?

[ProfJCNash]

I've not used it for group data, and suspect that the code to generate 
derivatives cannot cope with the bracket syntax. If you can rewrite the 
equation without the brackets, you could get the derivatives and solve 
that way. This will probably mean having a "translation" routine to glue 
things together.


JN

On 15-09-21 12:22 PM, Jianling Fan wrote:

Thanks Prof. Nash,

Sorry for late reply. I am learning and trying to use your nlmrt
package since I got your email. It works good to mask a parameter in
regression but seems does work for my equation. I think the problem is
that the parameter I want to mask is a group-specific parameter and I
have a "[]" syntax in my equation. However, I don't have your 2014
book on hand and couldn't find it in our library. So I am wondering if
nlxb works for group data?
Thanks a lot!

following is my code and I got a error form it.


fitdp1<-nlxb(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,

 + start =c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
Rm5=1.01, Rm6=1, d50=20, c=-1),
 + masked=c("Rm6"))

Error in deriv.default(parse(text = resexp), names(start)) :
   Function '`[`' is not in the derivatives table


Best regards,

Jianling


On 20 September 2015 at 12:56, ProfJCNash <profjcn...@gmail.com> wrote:

I posted a suggestion to use nlmrt package (function nlxb to be precise),
which has masked (fixed) parameters. Examples in my 2014 book on Nonlinear
parameter optimization with R tools. However, I'm travelling just now, or
would consider giving this a try.

JN


On 15-09-20 01:19 PM, Jianling Fan wrote:


no, I am doing a regression with 6 group data with 2 shared parameters
and 1 different parameter for each group data. the parameter I want to
coerce is for one group. I don't know how to do it. Any suggestion?

Thanks!

On 19 September 2015 at 13:33, Jeff Newmiller <jdnew...@dcn.davis.ca.us>
wrote:


Why not rewrite the function so that value is not a parameter?

---
Jeff NewmillerThe .   .  Go
Live...
DCN:<jdnew...@dcn.davis.ca.us>Basics: ##.#.   ##.#.  Live
Go...
Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.
rocks...1k

---
Sent from my phone. Please excuse my brevity.

On September 18, 2015 9:54:54 PM PDT, Jianling Fan
<fanjianl...@gmail.com> wrote:


Hello, everyone,

I am using a nls regression with 6 groups data. I am trying to coerce
a parameter to 1 by using a upper and lower statement. but I always
get an error like below:

Error in ifelse(internalPars < upper, 1, -1) :
   (list) object cannot be coerced to type 'double'

does anyone know how to fix it?

thanks in advance!

My code is below:




dproot


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


fitdp<-nls(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,


+ start = list(Rm=c(1.01, 1.01, 1.01, 6.65,1.01,1), d50=20, c=-1))


summary(fitdp)



Formula: den ~ Rm[ref]/(1 + (depth/d50)^c)

Parameters:
 Estimate Std. Error t value Pr(>|t|)
Rm1  1.125600.07156   15.73 3.84e-14 ***
Rm2  1.576430.11722   13.45 1.14e-12 ***
Rm3  1.106970.07130   15.53 5.11e-14 ***
Rm4  7.239250.20788   34.83  < 2e-16 ***
Rm5  1.145160.07184   15.94 2.87e-14 ***
Rm6  1.036580.05664   18.30 1.33e-15 ***
d50 22.694261.03855   21.85  < 2e-16 ***
c   -1.597960.15589  -10.25 3.02e-10 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.1094 on 24 degrees of freedom

Number of iterations to convergence: 8
Achieved convergence tolerance: 9.374e-06


fitdp1<-nls(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,


algorithm="port",
+ start = list(Rm=c(1.01, 1.01, 1.01, 6.65, 1.01, 1), d50=20, c=-1),
+ lower = list(Rm=c(1.01, 1.01, 1.01, 6.65, 1.01, 1), d50=20, c=-1),
+ upper = list(Rm=c(2.1, 2.2, 2.12, 12.5, 2.3, 1), d50=50, c=1))

Error in ifelse(internalPars < upper, 1, -1) :
   (list) object cannot be coerc

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

[ProfJCNash]

Apologies for replying and overlapping Gabor's contribution, which 
actually did the work!


Best, JN

On 15-09-21 01:03 PM, Jianling Fan wrote:

Thanks Gabor,

  That works good to rewrite the express the formula!

Thanks a lot!

Regards,

Jianling

On 21 September 2015 at 10:43, Gabor Grothendieck
<ggrothendi...@gmail.com> wrote:

Express the formula in terms of simple operations like this:

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

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

where we used this input:

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

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



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


Thanks Prof. Nash,

Sorry for late reply. I am learning and trying to use your nlmrt
package since I got your email. It works good to mask a parameter in
regression but seems does work for my equation. I think the problem is
that the parameter I want to mask is a group-specific parameter and I
have a "[]" syntax in my equation. However, I don't have your 2014
book on hand and couldn't find it in our library. So I am wondering if
nlxb works for group data?
Thanks a lot!

following is my code and I got a error form it.


fitdp1<-nlxb(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,

 + start =c(Rm1=1.01, Rm2=1.01, Rm3=1.01, Rm4=6.65,
Rm5=1.01, Rm6=1, d50=20, c=-1),
 + masked=c("Rm6"))

Error in deriv.default(parse(text = resexp), names(start)) :
   Function '`[`' is not in the derivatives table


Best regards,

Jianling


On 20 September 2015 at 12:56, ProfJCNash <profjcn...@gmail.com> wrote:

I posted a suggestion to use nlmrt package (function nlxb to be
precise),
which has masked (fixed) parameters. Examples in my 2014 book on
Nonlinear
parameter optimization with R tools. However, I'm travelling just now,
or
would consider giving this a try.

JN


On 15-09-20 01:19 PM, Jianling Fan wrote:


no, I am doing a regression with 6 group data with 2 shared parameters
and 1 different parameter for each group data. the parameter I want to
coerce is for one group. I don't know how to do it. Any suggestion?

Thanks!

On 19 September 2015 at 13:33, Jeff Newmiller
<jdnew...@dcn.davis.ca.us>
wrote:


Why not rewrite the function so that value is not a parameter?


---
Jeff NewmillerThe .   .  Go
Live...
DCN:<jdnew...@dcn.davis.ca.us>Basics: ##.#.   ##.#.  Live
Go...
Live:   OO#.. Dead: OO#..
Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.
rocks...1k


---
Sent from my phone. Please excuse my brevity.

On September 18, 2015 9:54:54 PM PDT, Jianling Fan
<fanjianl...@gmail.com> wrote:


Hello, everyone,

I am using a nls regression with 6 groups data. I am trying to coerce
a parameter to 1 by using a upper and lower statement. but I always
get an error like below:

Error in ifelse(internalPars < upper, 1, -1) :
   (list) object cannot be coerced to type 'double'

does anyone know how to fix it?

thanks in advance!

My code is below:




dproot


depth   den ref
1 20 0.573   1
2 40 0.780   1
3 60 0.947   1
4 80 0.990   1
5100 1.000   1
6 10 0.600   2
7 20 0.820   2
8 30 0.930   2
9 40 1.000   2
1020 0.480   3
1140 0.734   3
1260 0.961   3
1380 0.998   3
14   100 1.000   3
1520 3.2083491   4
1640 4.9683383   4
1760 6.2381133   4
1880 6.5322348   4
19   100 6.5780660 

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

[ProfJCNash]

Besides this, using bounds to fix (also called "mask") parameters is 
generally a very bad idea. Some optimization methods allow this 
explicitly. For nonlinear least squares nlmrt package has it, but I'm 
not sure I fully documented the process. For optimization, Rvmmin and 
Rcgmin both allow masks, but again the documentation is not fully 
developed. There is a section in my 2014 book Nonlinear Parameter 
Optimization in R. I'm on vacation, and don't have the page refs. at 
this moment, however.


I know there are some other packages that include the possibility of 
fixed parameters, and perhaps others could give pointers, as I don't 
have the references with me.


JN

On 15-09-19 05:10 AM, Sarah Goslee wrote:

You need to read the help more closely. start should be a list, as you've
done, but upper and lower should be vectors instead. Which is exactly what
your error message is telling you.

lower, upper

vectors of lower and upper bounds, replicated to be as long as start. If
unspecified, all parameters are assumed to be unconstrained. Bounds can
only be used with the "port" algorithm. They are ignored, with a warning,
if given for other algorithms.
Sarah

On Saturday, September 19, 2015, Jianling Fan  wrote:


Hello, everyone,

I am using a nls regression with 6 groups data. I am trying to coerce
a parameter to 1 by using a upper and lower statement. but I always
get an error like below:

Error in ifelse(internalPars < upper, 1, -1) :
   (list) object cannot be coerced to type 'double'

does anyone know how to fix it?

thanks in advance!

My code is below:




dproot

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


fitdp<-nls(den~Rm[ref]/(1+(depth/d50)^c),data=dproot,

+ start = list(Rm=c(1.01, 1.01, 1.01, 6.65,1.01,1), d50=20, c=-1))

summary(fitdp)


Formula: den ~ Rm[ref]/(1 + (depth/d50)^c)

Parameters:
 Estimate Std. Error t value Pr(>|t|)
Rm1  1.125600.07156   15.73 3.84e-14 ***
Rm2  1.576430.11722   13.45 1.14e-12 ***
Rm3  1.106970.07130   15.53 5.11e-14 ***
Rm4  7.239250.20788   34.83  < 2e-16 ***
Rm5  1.145160.07184   15.94 2.87e-14 ***
Rm6  1.036580.05664   18.30 1.33e-15 ***
d50 22.694261.03855   21.85  < 2e-16 ***
c   -1.597960.15589  -10.25 3.02e-10 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.1094 on 24 degrees of freedom

Number of iterations to convergence: 8
Achieved convergence tolerance: 9.374e-06


fitdp1<-nls(den~Rm[ref]/(1+(depth/d50)^c),data=dproot, algorithm="port",

+ start = list(Rm=c(1.01, 1.01, 1.01, 6.65, 1.01, 1), d50=20, c=-1),
+ lower = list(Rm=c(1.01, 1.01, 1.01, 6.65, 1.01, 1), d50=20, c=-1),
+ upper = list(Rm=c(2.1, 2.2, 2.12, 12.5, 2.3, 1), d50=50, c=1))

Error in ifelse(internalPars < upper, 1, -1) :
   (list) object cannot be coerced to type 'double'

__
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] Optimization Grid Search Slow

[ProfJCNash]

optimx does nothing to speed up optim or the other component optimizers. 
In fact, it does a lot of checking and extra work to improve reliability 
and add KKT tests that actually slow things down. The purpose of optimx 
is to allow comparison of methods and discovery of improved approaches 
to a problem. Is your function computing correctly?


Assuming you've got a correct function, then spending some time to speed 
up the function (I've found FORTRAN speediest) is likely your best hope.


JN



On 15-09-17 01:55 PM, Patzelt, Edward wrote:

R Help -

I am trying to use a grid search for a 2 free parameter reinforcement
learning model and the grid search is incredibly slow. I've used optimx but
can't seem to get reasonable answers. Is there a way to speed up this grid
search dramatically?


dat <- structure(list(choice = c(0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1,
  1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1,
0, 1, 0, 1, 0, 1, 0,
  0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1,
0, 0, 1, 0, 0, 1, 1,
  1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0,
0, 1, 0, 0, 0, 0, 1,
  1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1,
1, 0, 0, 0, 0, 0, 0,
  1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1,
1, 0, 0, 0, 0, 0, 1,
  1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1,
1, 0, 0, 1, 1, 0, 0,
  0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1,
0, 0, 1, 0, 0, 0, 0,
  1, 0, 1, 1, 1, 0), reward = c(0L, 0L, 0L,
0L, 1L, 1L, 0L, 0L,
1L, 0L, 0L,
0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L,
1L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 1L,
1L, 0L, 1L,
0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L,
0L, 0L, 1L,
1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L,
1L, 1L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 0L, 0L,
0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L,
1L, 0L, 0L,
1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L,
0L, 1L, 0L,
0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 1L, 0L,
1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L,
0L, 0L, 1L,
0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L), RepNum = c(1L,

1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,

1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,

1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,

1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,

2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,

2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,

2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,

2L, 2L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,

3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,

3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,

3L, 3L, 3L, 3L)), .Names
= c("choice", "reward", "RepNum"), row.names = c(NA,


  165L), class =
"data.frame")


CNTRACSID <- 0; subjectFit <- 0;
pLlist <- 0; pRlist <- 0; logLikelihood <- 0; trialProb <- 0;

hmmFunc <- function(delta, temperature){

   pLlist = 1
   pRlist = 1
   block = 0
   for (i in 1:length(dat$choice))
   {
 if (dat$RepNum[i] != block)
 {
   pL = 0.5
   pR = 0.5
   block = dat$RepNum[i]
 }
 # Markov Transitions
 pL <- pL*(1-delta) + pR*delta
 pR <- 1-pL
 # Apply feedback
 #denom <- p(F|L,C) * p(L) + p(F|R,C) * p(R)

 pflc <- ifelse(dat$choice[i] == dat$reward[i], .8, .2)
 pfrc <- 1 - pflc
 denom <- pflc * pL + pfrc * pR

 # What's the new belief given observation
 posteriorL <- pflc * pL/denom
 

Re: [R] Non linear regression - Von Bertalanffy Growth Function - singular gradient matrix at initial parameter estimates

[ProfJCNash]

Packages nlmrt or minpack.lm use a Marquardt method. minpack.lm won't 
proceed if the Jacobian singularity is at the starting point as far as 
I'm aware, but nlxb in nlmrt can sometimes get going. It has a policy 
that is aggressive in trying to improve the sum of squares, so will use 
more effort than nls when both work.


JN

On 15-08-18 12:08 PM, Xochitl CORMON wrote:

Dear all,

I am trying to estimate VBGF parameters K and Linf using non linear
regression and nls(). First I used a classic approach where I estimate
both parameters together as below with alkdyr being a subset per year
of my age-length-key database and running in a loop.

vbgf.par - nls(Lgtcm ~  Linf *(1 - exp(-K * (Age - tzero))), start =
c(K= 0.07, Linf = 177.1), data=alkdyr)

I obtain an estimation of both parameters that are strongly correlated.
Indeed after plotting Linf ~ K and fitting a linear regression I obtain
a function (Linf = a + b*K) with R2= 0.8 and a = 215, b = -763.

In this context, to take into account explicitly correlation between
parameters, I decided to fit a new non linear regression derivate from
VBGF but where Linf is expressed depending on K (I am most interested in
K). To do so, I tried this model:
vbgf.par - nls(Lgtcm ~  (a + (b*k)) *(1 - exp(-k * (Age - tzero))),
start = c(k= 0.07, a= 215, b=-763), data=alkdyr)

Unfortunately at this point I cannot go further as I get the error
message singular gradient matrix at initial parameter estimates.

I tried to use alg= plinear (which I am not sure I understand properly
yet). If I give a starting value for a and b only, I have an error
message stating step factor below minFactor (even when minFactor is
set to 1000).

Any help will be more than welcome as this is quite urgent

Best,

Xochitl C.






__
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] Error: Gradient function might be wrong - check it! OPTMX

[ProfJCNash]

There are tolerances in the checks, and sometimes scaling of the problem
leads to false positives on the checks.

Try control = list(starttest=FALSE, etc.) to suppress the test.

JN

On 15-08-13 01:20 PM, Olu Ola via R-help wrote:
 Hello,
 I am trying to estimate a non-linear GMM in R using Optimx. However, when I 
 do the start test to compare my analytical gradient with the numerical 
 gradient, I get the following error message
 
 Error: Gradient function might be wrong - check it!
 
 Below are the print out of my analytical and numerical gradient. They are 
 close to an extent and I want to use my analytical gradient. The first column 
 is the analytical gradient while the second column is the numerical gradient
 
  numgrad
  [1,]   -9862.118   -9862.103
  [2,]  -20855.001  -20854.985
  [3,]  -12277.764  -12277.767
  [4,]3332.8643332.865
  [5,]3562.8083562.804
  [6,]  -14971.913  -14971.915
  [7,]   27394.457   27394.451
  [8,]  -15634.531  -15634.532
  [9,]  -16893.399  -16893.402
 [10,]   -1330.469   -1330.468
 [11,]   16877.994   16877.991
 [12,]   -9406.643   -9406.643
 [13,]  -40451.881  -40451.881
 [14,]  -48052.941  -48052.944
 [15,]  -29997.114  -29997.116
 [16,]   34546.011   34546.002
 [17,]  -65540.191  -65540.200
 [18,]  -24400.970  -24400.972
 [19,]  -26655.760  -26655.760
 [20,]  -59128.576  -59128.576
 [21,]  -40026.926  -40026.931
 [22,]  -34780.662  -34780.655
 [23,] -119284.474 -119284.473
 [24,]   22336.013   22336.025
 [25,]  -40665.139  -40665.139
 [26,]  -21596.959  -21596.959
 [27,]  -33859.274  -33860.747
 [28,]  -30536.938  -30539.807
 [29,]  -95411.022  -95416.182
 [30,]  -33016.450  -33018.690
 [31,]  -90052.772  -90059.612
 [32,]8008.6848009.342
 [33,]  -45084.823  -45086.720
 [34,]  -34071.584  -34074.802
 [35,] -259960.860 -259990.135
 [36,]  -69986.472  -69990.537
 [37,]   24870.530   24874.012
 
 The code for my optimization is as follows:
 
 
 nlgmm = optimx(par=b0, fn=obj,gr=gra, method =BFGS, itnmax=1, 
 control=list(follow.on = TRUE,usenumDeriv=FALSE,kkttol=10,starttests=TRUE, 
 save.failures=TRUE, trace=0))
 
 Is there a way to turn off the gradient check? this is because if I do not 
 run the start test, it will still eventually display the error message. Any 
 help to cross this bridge will be greatly appreciated.
 
 Thank you
 
 __
 R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


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

Re: [R] nls in r

[ProfJCNash]

With package nlmrt, I get a solution, but the Jacobian is essentially
singular, so the model may not be appropriate. You'll need to read the
documentation to learn how to interpret the Jacobian singular values. Or
Chapter 6 of my book Nonlinear parameter optimization with R tools.

Here's the script. Note nlxb is a little different from nls() and
requires a data frame for the data.

library(stats)
x=c(30:110)
y=c(0.000760289, 0.000800320, 0.000830345, 0.000840353, 0.000860370,
0.000910414, 0.000990490, 0.001090594, 0.001200721, 0.001350912,
0.001531172, 0.001751533, 0.001961923, 0.002192402, 0.002463031,
0.002793899, 0.003185067, 0.003636604, 0.004148594, 0.004721127,
0.005394524, 0.006168989, 0.007014544, 0.007870894, 0.008758242,
0.009656474, 0.010565620, 0.011485709, 0.012396520, 0.013308162,
0.014271353, 0.015326859, 0.016525802, 0.017889059, 0.019447890,
0.021223636, 0.023145810, 0.025174229, 0.027391752, 0.029645106,
0.032337259, 0.035347424, 0.039375125, 0.043575783, 0.048003973,
0.052926206, 0.058307309, 0.064581189, 0.071947231, 0.080494476,
0.089891885, 0.100671526, 0.111971207, 0.124237571, 0.137975539,
0.153629149, 0.171239194, 0.190712664, 0.212079979, 0.235026373,
0.259457493, 0.282867017, 0.307830359, 0.334773680, 0.364001719,
0.395742526, 0.425596389, 0.458391314, 0.494558651, 0.534657357,
0.579354317, 0.616075034, 0.656680256, 0.701804548, 0.752133146,
0.808558032, 0.872226001, 0.944664487, 1.027837007, 1.124484096,
1.238426232)

vdata- data.frame(x=x, y=y)

a=0.0189
b=0.14328
delta=0.0005

# fit = nls(y ~
a*(exp(b*(x+0.5)))*((delta*b)/((delta*b)+(a*(exp(b*(x+0.5))-1^(0.5),
# start=list(a=a,b=b, delta=delta), trace=TRUE)
# predict(fit)
# plot(x,y,col=red, xlab=Usia,ylab=expression(paste(mu)))
# lines(x,predict(fit), col=blue)
# legend(topleft,
# c(expression(paste(mu)),Fit),col=c(red,blue),lty=1:1)

library(nlmrt)

 vformula - y ~
a*(exp(b*(x+0.5)))*((delta*b)/((delta*b)+a*(exp(b*(x+0.5))-1)))^(0.5)

fitjn = nlxb(formula= vformula, start=list(a=a,b=b,delta=delta),
trace=TRUE, data=vdata)
fitjn

Best, JN


On 15-08-12 06:37 AM, vidya wrote:
 I get this error 
 Error in numericDeriv(form[[3L]], names(ind), env) : 
   Missing value or an infinity produced when evaluating the model
 I was replace the starting value but still get error.
 
 Here is my code:
 library(stats)
 x=c(30:110)
 y=c(0.000760289, 0.000800320, 0.000830345, 0.000840353, 0.000860370,  
 0.000910414, 0.000990490, 0.001090594, 0.001200721, 0.001350912, 
 0.001531172, 0.001751533, 0.001961923, 0.002192402, 0.002463031, 
 0.002793899, 0.003185067, 0.003636604, 0.004148594, 0.004721127, 
 0.005394524, 0.006168989, 0.007014544, 0.007870894, 0.008758242, 
 0.009656474, 0.010565620, 0.011485709, 0.012396520, 0.013308162, 
 0.014271353, 0.015326859, 0.016525802, 0.017889059, 0.019447890, 
 0.021223636, 0.023145810, 0.025174229, 0.027391752, 0.029645106, 
 0.032337259, 0.035347424, 0.039375125, 0.043575783, 0.048003973, 
 0.052926206, 0.058307309, 0.064581189, 0.071947231, 0.080494476, 
 0.089891885, 0.100671526, 0.111971207, 0.124237571, 0.137975539, 
 0.153629149, 0.171239194, 0.190712664, 0.212079979, 0.235026373,
 0.259457493, 0.282867017, 0.307830359, 0.334773680, 0.364001719, 
 0.395742526, 0.425596389, 0.458391314, 0.494558651, 0.534657357, 
 0.579354317, 0.616075034, 0.656680256, 0.701804548, 0.752133146, 
 0.808558032, 0.872226001, 0.944664487, 1.027837007, 1.124484096, 
 1.238426232)
 
 a=0.0189
 b=0.14328
 delta=0.0005
 
 fit = nls(y ~
 a*(exp(b*(x+0.5)))*((delta*b)/((delta*b)+(a*(exp(b*(x+0.5))-1^(0.5), 
 start=list(a=a,b=b, delta=delta))
 predict(fit)
 plot(x,y,col=red, xlab=Usia,ylab=expression(paste(mu)))
 lines(x,predict(fit), col=blue)
 legend(topleft,
 c(expression(paste(mu)),Fit),col=c(red,blue),lty=1:1)
 
 
 I really appreciate for the helps. Thank you.
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/nls-in-r-tp4711012.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


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

Re: [R] how to install and use Rcplex package

[ProfJCNash]

It's important to look at the CRAN documentation, where it is quite
clear there are NO Windows binaries for this package.

My suggestion -- set up the (freeware) VirtualBox or a similar Virtual
Machine environment, install a Linux OS virtually, and install there. If
there is no Windows binary on CRAN, there is likely a very good reason.

You may actually find that something like Lubuntu or Linux Mint are more
convenient for working with R.

JN

On 15-07-24 09:54 AM, Yasin Gocgun wrote:
 Hi,
 
 Thanks for your email. When I use that line, I received the following
 messages from R:
 
 package ‘Rcplex’ is available as a source package but not as a binary
 
 package ‘Rcplex’ is not available (for R version 3.1.1)
 
 So, do you see what the problem is?
 
 Thank you,
 
 Yasin
 
 On Fri, Jul 24, 2015 at 4:04 PM, Rui Barradas ruipbarra...@sapo.pt wrote:
 Hello,

 Try

 install.packages(Rcplex)

 That's it.

 Hope this helps,

 Rui Barradas

 Em 24-07-2015 12:20, Yasin Gocgun escreveu:

 Hi,

 I need to know how to completely install Rcplex to be able to use
 cplex in R. According to the instructions given in the following
 webpage, I edited the Makevars.win file:

 https://cran.r-project.org/web/packages/Rcplex/INSTALL

 The revised Makevars.win file includes only the following lines:

 PKG_CPPFLAGS=-I/C:/ilog/CPLEX_Studio122/include

 PKG_LIBS=-L/C:/ilog/CPLEX_Studio122/cplex/lib/x64_windows_vs2008/stat_mda-lcplex122
 -lm

 First, I am not sure whether the above lines are completely correct.
 Second, I don'y know what is the next step in installing and using
 Rcplex. So can someone help me about this issue?

 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.

Re: [R] nls singular gradient matrix - fit parameters in integral's upper limits

[ProfJCNash]

The list rejects almost all attachments.

You could dput the data and put it in your posting.

You may also want to try a Marquardt solver. In R from my nlmrt or
compiled in Kate Mullen's minpack.lm. They are slightly different in
flavour and the call is a bit different from nls.

JN

On 15-07-16 08:21 AM, Laura Teresa Corredor Bohórquez wrote:
 ---The las post rejected two files I had attached, so I modified
 it.---
 
 
 Hi. I am trying to make a nls fit for a little bit complicated expression that
 includes two integrals with two of the fit parameters in their upper limits.
 
 I got the error Error in nlsModel(formula, mf, start, wts) :
 singular gradient
 matrix at initial parameter estimates. First of all, I have searched
 already in the previous answers, but didn´t help. The parameters 
 initialization
 seems to be ok, I have tried to change the parameters but no one works. If
 my function has just one integral everything works very nicely, but when 
 adding
 a second integral term just got the error. I don´t believe the function is
 over-parametrized, as I have performed other fits with much more parameters
 and they worked. I did try to enclose the data but the attachment was
 rejected.
 
 The minimal example is the following:
 
 # read the data from a csv file
 dados = read.csv(file.csv, header=FALSE, stringsAsFactors=FALSE)
 x = 0*(1:97)
 y = 0*(1:97)
 for(i in 1:97){
   x[i] = dados[i,1]
   y[i] = dados[i,2]
 }
 integrand - function(X) {
   return(X^4/(2*sinh(X/2))^2)
 }
 fitting = function(T1, T2, N, D, x){
   int1 = integrate(integrand, lower=0, upper = T1)$value
   int2 = integrate(integrand, lower=0, upper = T2)$value
   return(N*(D/x)^2*(exp(D/x)/(1+exp(D/x))^2
 )+(448.956*(x/T1)^3*int1)+(299.304*(x/T2)^3*int2))
 }
 fit = nls(y ~ fitting(T1, T2, N, D, x),
 start=list(T1=400,T2=200,N=0.01,D=2))
 
 --For reference, the fit that worked is the following:
 
 # read the data from a csv file
 dados = read.csv(file.csv, header=FALSE, stringsAsFactors=FALSE)
 x = 0*(1:97)
 y = 0*(1:97)
 for(i in 1:97){
   x[i] = dados[i,1]
   y[i] = dados[i,2]
 }
 integrand - function(X) {
   return(X^4/(2*sinh(X/2))^2)
 }
 fitting = function(T1, N, D, x){
   int = integrate(integrand, lower=0, upper = T1)$value
   return(N*(D/x)^2*(exp(D/x)/(1+exp(D/x))^2 )+(748.26)*(x/T1)^3*int)
 }
 fit = nls(y ~ fitting(T1 , N, D, x), start=list(T1=400,N=0.01,D=2))
 
 
 I cannot figure out what happen. I need to perform this fit for three
 integral components, but even for two I have this problem. I appreciate so
 much your help. Thank you.
 
   [[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] Invalid URL for R documentation

[ProfJCNash]

I'm also getting a 404. Tried https just in case.

But

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

works and I can copy the link for R-intro and it works there.

http://cran.utstat.utoronto.ca/doc/manuals/r-release/R-intro.html

Very odd.

JN


On 15-07-06 11:24 AM, Paul wrote:
 http://cran.r-project.org/doc/manuals/r-release/R-intro.html

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

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

[ProfJCNash]

n163 - mpfr(163, 500)

is how I set up the number.

JN


On 15-07-04 05:10 PM, Ravi Varadhan wrote:
 What about numeric constants, like `163'?  

 eval(Pre(exp(sqrt(163)*pi), 120))does not work.

 Thanks,
 Ravi
 
 From: David Winsemius dwinsem...@comcast.net
 Sent: Saturday, July 4, 2015 1:12 PM
 To: Duncan Murdoch
 Cc: r-help; John Nash; Ravi Varadhan
 Subject: Re: [R] : Ramanujan and the accuracy of floating point computations 
 - using Rmpfr in R

 On Jul 4, 2015, at 12:20 AM, Duncan Murdoch murdoch.dun...@gmail.com 
 wrote:

 On 04/07/2015 8:21 AM, David Winsemius wrote:
 On Jul 3, 2015, at 11:05 PM, Duncan Murdoch murdoch.dun...@gmail.com 
 wrote:

 On 04/07/2015 3:45 AM, David Winsemius wrote:
 On Jul 3, 2015, at 5:08 PM, David Winsemius dwinsem...@comcast.net 
 wrote:



 It doesn’t appear to me that mpfr was ever designed to handle 
 expressions as the first argument.
 This could be a start. Obviously one would wnat to include code to do 
 other substitutions probably using the all.vars function to pull out the 
 other “constants” and ’numeric’ values to make them of equivalent 
 precision. I’m guessing you want to follow the parse-tree and then 
 testing the numbers for integer-ness and then replacing by paste0( 
 “mpfr(“, val, “L, “, prec,”)” )

 Pre - function(expr, prec){ sexpr - deparse(substitute(expr) )
 Why deparse?  That's almost never a good idea.  I can't try your code (I
 don't have mpfr available), but it would be much better to modify the
 expression than the text representation of it.  For example, I think
 your code would modify strings containing pi, or variables with those
 letters in them, etc.  If you used substitute(expr) without the
 deparse(), you could replace the symbol pi with the call to the Const
 function, and be more robust.

 Really? I did try. I was  fairly sure that someone could do better but I 
 don’t see an open path along the lines you suggest. I’m pretty sure I 
 tried `substitute(expr, list(pi= pi))` when `expr` had been the formal 
 argument and got disappointed because there is no `pi` in the expression 
 `expr`. I _thought_ the problem was that `substitute` does not evaluate 
 its first argument, but I do admit to be pretty much of a klutz with this 
 sort of programming. I don’t think you need to have mpfr installed in 
 order to demonstrate this.
 The substitute() function really does two different things.
 substitute(expr) (with no second argument) grabs the underlying
 expression out of a promise.  substitute(expr, list(pi = pi)) tries to
 make the substitution in the expression expr, so it doesn't see pi”.
 Thank you. That was really helpful. I hope it “sticks” to  sufficiently 
 durable set of neurons.
 This should work:

 do.call(substitute, list(expr = substitute(expr), env=list(pi =
 Const(pi, 120

 (but I can't evaluate the Const function to test it).
 The expression `pi` as the argument to Const only needed to be 
 character()-ized and then evaluation on that result succeeded:

 library(mpfr)
 #  Pre - function(expr, prec){ do.call(substitute,
 list(expr = substitute(expr),
  env=list(pi = Const(pi, 
 prec }

 Pre(exp(pi),120)
 Error in Const(pi, prec) :
   'name' must be one of 'pi', 'gamma', 'catalan', 'log2'


  Pre - function(expr, prec){ do.call(substitute,
  list(expr = substitute(expr),
   env=list(pi = Const('pi', prec 
 }

 Pre(exp(pi),120)
 exp(list(S4 object of class mpfr1))
 eval( Pre(exp(pi),120) )
 1 'mpfr' number of precision  120   bits
 [1] 23.140692632779269005729086367948547394

 So the next step for delivering Ravi’s request would be to add the rest of 
 the ‘Const’ options:

 Pre - function(expr, prec){ do.call(substitute,
  list(expr = substitute(expr),
  env=list(pi = Const('pi', prec),
  catalan=Const('catalan', prec),
  log2=Const('log2', prec),
  gamma=Const('gamma', prec }


 eval(Pre(exp(gamma)))
 Error in stopifnot(is.numeric(prec)) :
   argument prec is missing, with no default
 eval(Pre(exp(gamma), 120))
 1 'mpfr' number of precision  120   bits
 [1] 1.781072417990197985236504103107179549

 This works:

 do.call(substitute, list(expr = substitute(expr), env = list(pi =
 quote(Const(pi, 120)

 because it never evaluates the Const function, it just returns some code
 that you could modify more if you liked.

 Duncan Murdoch
 David Winsemius, MD
 Alameda, CA, USA


__
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