[R] Optimization under an absolute value constraint
I need to optimize a multivariate function f(w, x, y, z, ...) under an absolute value constraint. For instance: min { (2x+y) (w-z) } under the constraint: |w| + |x| + |y| + |z| = 1.0 . Is there any R function that does this? Thank you for your help! Phil Xiang __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization under an absolute value constraint
this should be possible in the lasso2 package. url:www.econ.uiuc.edu/~rogerRoger Koenker email[EMAIL PROTECTED]Department of Economics vox: 217-333-4558University of Illinois fax: 217-244-6678Champaign, IL 61820 On Sep 7, 2007, at 1:17 PM, Phil Xiang wrote: I need to optimize a multivariate function f(w, x, y, z, ...) under an absolute value constraint. For instance: min { (2x+y) (w-z) } under the constraint: |w| + |x| + |y| + |z| = 1.0 . Is there any R function that does this? Thank you for your help! Phil Xiang __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting- guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization under an absolute value constraint
On 9/7/07, Phil Xiang [EMAIL PROTECTED] wrote: I need to optimize a multivariate function f(w, x, y, z, ...) under an absolute value constraint. For instance: min { (2x+y) (w-z) } under the constraint: |w| + |x| + |y| + |z| = 1.0 . Is there any R function that does this? Thank you for your help! I think that the minimum value of the function f(x) := -2*x*(1-x), with 0 = x = 1 is also the minimum value of the objective function of your problem (but correct me if I am wrong). Thus, x y w z -0.50 0 -0.5 -0.50 0.1 -0.4 -0.50 0.3 -0.2 0.5 0 -0.50 -0.50 0.5 0 0.5 0 -0.40.1 0.5 0 -0.20.3 0.5 0 0 0.5 are all solutions for your problem. Paul __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization problem
Try this. 1. following Ben remove the Randalstown point and reset the levels of the Location factor. 2. then replace solve with ginv so it uses the generalized inverse to calculate the hessian: alan2 - subset(alan, subset = Location != Randalstown) alan2$Location - factor(as.character(alan2$Location)) library(MASS) solve - ginv zinb.zc - zicounts(resp=Scars~.,x =~Location + Lar + Mass + Lar:Mass + Location:Mass,z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data = alan2) rm(solve) On 8/21/07, Ben Bolker [EMAIL PROTECTED] wrote: (Hope this gets threaded properly. Sorry if it doesn't.) Gabor: Lac and Lacfac being the same is irrelevant, wouldn't produce NAs (but would produce something like a singular Hessian and maybe other problems) -- but they're not even specified in this model. The bottom line is that you have a location with a single observation, so the GLM that zicounts runs to get the initial parameter values has an unestimable location:mass interaction for one location, so it gives an NA, so optim complains. In gruesome detail: ## set up data scardat = read.table(scars.dat,header=TRUE) library(zicounts) ## try to run model zinb.zc - zicounts(resp=Scars~., x =~Location + Lar + Mass + Lar:Mass + Location:Mass, z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data=scardat) ## tried to debug this by dumping zicounts.R to a file, modifying ## it to put a trace argument in that would print out the parameters ## and log-likelihood for every call to the log-likelihood function. dump(zicounts,file=zicounts.R) source(zicounts.R) zinb.zc - zicounts(resp=Scars~., x =~Location + Lar + Mass + Lar:Mass + Location:Mass, z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data=scardat,trace=TRUE) ## this actually didn't do any good because the negative log-likelihood ## function never gets called -- as it turns out optim() barfs when it ## gets its initial values, before it ever gets to evaluating the log-likelihood ## check the glm -- this is the equivalent of what zicounts does to ## get the initial values of the x parameters p1 - glm(Scars~Location + Lar + Mass + Lar:Mass + Location:Mass, data=scardat,family=poisson) which(is.na(coef(p1))) ## find out what the deal is table(scardat$Location) scar2 = subset(scardat,Location!=Randalstown) ## first step to removing the bad point from the data set -- but ... table(scar2$Location) ## it leaves the Location factor with the same levels, so ## now we have ZERO counts for one location: ## redefine the factor to drop unused levels scar2$Location - factor(scar2$Location) ## OK, looks fine now table(scar2$Location) zinb.zc - zicounts(resp=Scars~., x =~Location + Lar + Mass + Lar:Mass + Location:Mass, z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data=scar2) ## now we get another error (system is computationally singular when ## trying to compute Hessian -- overparameterized?) Not in any ## trivial way that I can see. It would be nice to get into the guts ## of zicounts and stop it from trying to invert the Hessian, which is ## I think where this happens. In the meanwhile, I have some other ideas about this analysis (sorry, but you started it ...) Looking at the data in a few different ways: library(lattice) xyplot(Scars~Mass,groups=Location,data=scar2,jitter=TRUE, auto.key=list(columns=3)) xyplot(Scars~Mass|Location,data=scar2,jitter=TRUE) xyplot(Scars~Lar,groups=Location,data=scar2, auto.key=list(columns=3)) xyplot(Scars~Mass|Lar,data=scar2) xyplot(Scars~Lar|Location,data=scar2) Some thoughts: (1) I'm not at all sure that zero-inflation is necessary (see Warton 2005, Environmentrics). This is a fairly small, noisy data set without huge numbers of zeros -- a plain old negative binomial might be fine. I don't actually see a lot of signal here, period (although there may be some) ... there's not a huge range in Lar (whatever it is -- the rest of the covariates I think I can interpret). It would be tempting to try to fit location as a random effect, because fitting all those extra degrees of freedom is going to kill you. On the other hand, GLMMs are a bit hairy. cheers Ben __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization problem
Hello Folks, Very new to R so bear with me, running 5.2 on XP. Trying to do a zero-inflated negative binomial regression on placental scar data as dependent. Lactation, location, number of tick larvae present and mass of mouse are independents. Dataframe and attributes below: Location Lac Scars Lar Mass Lacfac 1 Tullychurry 0 0 15 13.87 0 2 Somerset 0 0 0 15.60 0 3 Tollymore 0 0 3 16.43 0 4 Tollymore 0 0 0 16.55 0 5 Caledon 0 0 0 17.47 0 6 Hillsborough 1 5 0 18.18 1 7 Caledon 0 0 1 19.06 0 8 Portglenone 0 4 0 19.10 0 9 Portglenone 0 5 0 19.13 0 10Tollymore 0 5 3 19.50 0 11 Hillsborough 1 5 0 19.58 1 12 Portglenone 0 4 0 19.76 0 13 Caledon 0 8 0 19.97 0 14 Hillsborough 1 4 0 20.02 1 15 Tullychurry 0 3 3 20.13 0 16 Hillsborough 1 5 0 20.18 1 17 LoughNavar 1 5 0 20.20 1 18Tollymore 0 0 1 20.24 0 19 Hillsborough 1 5 0 20.48 1 20 Caledon 0 4 1 20.56 0 21 Caledon 0 3 2 20.58 0 22Tollymore 0 4 3 20.58 0 23Tollymore 0 0 2 20.88 0 24 Hillsborough 1 0 0 21.01 1 25 Portglenone 0 5 0 21.08 0 26 Tullychurry 0 2 5 21.28 0 27 Ballysallagh 1 4 0 21.59 1 28 Caledon 0 0 1 21.68 0 29 Hillsborough 1 5 0 22.09 1 30 Tullychurry 0 5 5 22.28 0 31 Tullychurry 1 6 75 22.43 1 32 Ballysallagh 1 5 0 22.57 1 33 Ballysallagh 1 4 0 22.67 1 34 LoughNavar 1 5 3 22.71 1 35 Hillsborough 1 4 0 23.01 1 36 Caledon 0 0 3 23.08 0 37 LoughNavar 1 5 0 23.53 1 38 Ballysallagh 1 4 0 23.55 1 39 Portglenone 1 6 0 23.61 1 40 Mt.Stewart 0 3 0 23.70 0 41 Somerset 0 5 0 23.83 0 42 Ballysallagh 1 5 0 23.93 1 43 Ballysallagh 1 5 0 24.01 1 44 Caledon 0 0 3 24.14 0 45 LoughNavar 0 6 0 24.30 0 46 LoughNavar 1 5 0 24.34 1 47 Hillsborough 1 4 0 24.45 1 48 Caledon 0 3 2 24.55 0 49 Tullychurry 0 5 44 24.83 0 50 Hillsborough 1 5 0 24.86 1 51 Ballysallagh 1 5 0 25.02 1 52 Tullychurry 0 0 9 25.27 0 53 Mt.Stewart 0 5 0 25.31 0 54 LoughNavar 1 4 8 25.43 1 55 Somerset 1 0 0 25.58 1 56 Hillsborough 1 5 0 25.82 1 57 Portglenone 1 2 0 26.02 1 58 Ballysallagh 1 5 0 26.19 1 59 Mt.Stewart 1 0 0 26.66 1 60 Randalstown 1 0 1 26.70 1 61 Somerset 0 4 0 27.01 0 62 Mt.Stewart 0 4 0 27.05 0 63 Somerset 0 3 0 27.10 0 64 Somerset 0 6 0 27.34 0 65 Somerset 0 0 0 27.87 0 66 LoughNavar 1 5 1 28.01 1 67 Tullychurry 1 6 42 28.55 1 68 Hillsborough 1 5 0 28.84 1 69 Portglenone 1 4 0 29.00 1 70 Somerset 1 4 0 31.87 1 71 Ballysallagh 1 5 0 33.06 1 72 LoughNavar 1 4 0 33.24 1 73 Somerset 1 4 0 33.36 1 alan : 'data.frame':73 obs. of 6 variables: $ Location: Factor w/ 10 levels Ballysallagh,..: 10 8 9 9 2 3 2 6 6 9 ... $ Lac : int 0 0 0 0 0 1 0 0 0 0 ... $ Scars : int 0 0 0 0 0 5 0 4 5 5 ... $ Lar : int 15 0 3 0 0 0 1 0 0 3 ... $ Mass: num 13.9 15.6 16.4 16.6 17.5 ... $ Lacfac : Factor w/ 2 levels 0,1: 1 1 1 1 1 2 1 1 1 1 ... The syntax I used to create the model is: zinb.zc - zicounts(resp=Scars~.,x =~Location + Lar + Mass + Lar:Mass + Location:Mass,z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data=alan) The error given is: Error in optim(par = parm, fn = neg.like, gr = neg.grad, hessian = TRUE, : non-finite value supplied by optim In addition: Warning message: fitted probabilities numerically 0 or 1 occurred in: glm.fit(zz, 1 - pmin(y, 1), family = binomial()) I understand this is a problem with the model I specified, could anyone help out?? Many thanks Alan Harrison Quercus Queen's University Belfast MBC, 97 Lisburn Road Belfast BT9 7BL T: 02890 972219 M: 07798615682 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization problem
Lac and Lacfac are the same. On 8/21/07, Alan Harrison [EMAIL PROTECTED] wrote: Hello Folks, Very new to R so bear with me, running 5.2 on XP. Trying to do a zero-inflated negative binomial regression on placental scar data as dependent. Lactation, location, number of tick larvae present and mass of mouse are independents. Dataframe and attributes below: Location Lac Scars Lar Mass Lacfac 1 Tullychurry 0 0 15 13.87 0 2 Somerset 0 0 0 15.60 0 3 Tollymore 0 0 3 16.43 0 4 Tollymore 0 0 0 16.55 0 5 Caledon 0 0 0 17.47 0 6 Hillsborough 1 5 0 18.18 1 7 Caledon 0 0 1 19.06 0 8 Portglenone 0 4 0 19.10 0 9 Portglenone 0 5 0 19.13 0 10Tollymore 0 5 3 19.50 0 11 Hillsborough 1 5 0 19.58 1 12 Portglenone 0 4 0 19.76 0 13 Caledon 0 8 0 19.97 0 14 Hillsborough 1 4 0 20.02 1 15 Tullychurry 0 3 3 20.13 0 16 Hillsborough 1 5 0 20.18 1 17 LoughNavar 1 5 0 20.20 1 18Tollymore 0 0 1 20.24 0 19 Hillsborough 1 5 0 20.48 1 20 Caledon 0 4 1 20.56 0 21 Caledon 0 3 2 20.58 0 22Tollymore 0 4 3 20.58 0 23Tollymore 0 0 2 20.88 0 24 Hillsborough 1 0 0 21.01 1 25 Portglenone 0 5 0 21.08 0 26 Tullychurry 0 2 5 21.28 0 27 Ballysallagh 1 4 0 21.59 1 28 Caledon 0 0 1 21.68 0 29 Hillsborough 1 5 0 22.09 1 30 Tullychurry 0 5 5 22.28 0 31 Tullychurry 1 6 75 22.43 1 32 Ballysallagh 1 5 0 22.57 1 33 Ballysallagh 1 4 0 22.67 1 34 LoughNavar 1 5 3 22.71 1 35 Hillsborough 1 4 0 23.01 1 36 Caledon 0 0 3 23.08 0 37 LoughNavar 1 5 0 23.53 1 38 Ballysallagh 1 4 0 23.55 1 39 Portglenone 1 6 0 23.61 1 40 Mt.Stewart 0 3 0 23.70 0 41 Somerset 0 5 0 23.83 0 42 Ballysallagh 1 5 0 23.93 1 43 Ballysallagh 1 5 0 24.01 1 44 Caledon 0 0 3 24.14 0 45 LoughNavar 0 6 0 24.30 0 46 LoughNavar 1 5 0 24.34 1 47 Hillsborough 1 4 0 24.45 1 48 Caledon 0 3 2 24.55 0 49 Tullychurry 0 5 44 24.83 0 50 Hillsborough 1 5 0 24.86 1 51 Ballysallagh 1 5 0 25.02 1 52 Tullychurry 0 0 9 25.27 0 53 Mt.Stewart 0 5 0 25.31 0 54 LoughNavar 1 4 8 25.43 1 55 Somerset 1 0 0 25.58 1 56 Hillsborough 1 5 0 25.82 1 57 Portglenone 1 2 0 26.02 1 58 Ballysallagh 1 5 0 26.19 1 59 Mt.Stewart 1 0 0 26.66 1 60 Randalstown 1 0 1 26.70 1 61 Somerset 0 4 0 27.01 0 62 Mt.Stewart 0 4 0 27.05 0 63 Somerset 0 3 0 27.10 0 64 Somerset 0 6 0 27.34 0 65 Somerset 0 0 0 27.87 0 66 LoughNavar 1 5 1 28.01 1 67 Tullychurry 1 6 42 28.55 1 68 Hillsborough 1 5 0 28.84 1 69 Portglenone 1 4 0 29.00 1 70 Somerset 1 4 0 31.87 1 71 Ballysallagh 1 5 0 33.06 1 72 LoughNavar 1 4 0 33.24 1 73 Somerset 1 4 0 33.36 1 alan : 'data.frame':73 obs. of 6 variables: $ Location: Factor w/ 10 levels Ballysallagh,..: 10 8 9 9 2 3 2 6 6 9 ... $ Lac : int 0 0 0 0 0 1 0 0 0 0 ... $ Scars : int 0 0 0 0 0 5 0 4 5 5 ... $ Lar : int 15 0 3 0 0 0 1 0 0 3 ... $ Mass: num 13.9 15.6 16.4 16.6 17.5 ... $ Lacfac : Factor w/ 2 levels 0,1: 1 1 1 1 1 2 1 1 1 1 ... The syntax I used to create the model is: zinb.zc - zicounts(resp=Scars~.,x =~Location + Lar + Mass + Lar:Mass + Location:Mass,z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data=alan) The error given is: Error in optim(par = parm, fn = neg.like, gr = neg.grad, hessian = TRUE, : non-finite value supplied by optim In addition: Warning message: fitted probabilities numerically 0 or 1 occurred in: glm.fit(zz, 1 - pmin(y, 1), family = binomial()) I understand this is a problem with the model I specified, could anyone help out?? Many thanks Alan Harrison Quercus Queen's University Belfast MBC, 97 Lisburn Road Belfast BT9 7BL T: 02890 972219 M: 07798615682 [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal,
Re: [R] Optimization problem
(Hope this gets threaded properly. Sorry if it doesn't.) Gabor: Lac and Lacfac being the same is irrelevant, wouldn't produce NAs (but would produce something like a singular Hessian and maybe other problems) -- but they're not even specified in this model. The bottom line is that you have a location with a single observation, so the GLM that zicounts runs to get the initial parameter values has an unestimable location:mass interaction for one location, so it gives an NA, so optim complains. In gruesome detail: ## set up data scardat = read.table(scars.dat,header=TRUE) library(zicounts) ## try to run model zinb.zc - zicounts(resp=Scars~., x =~Location + Lar + Mass + Lar:Mass + Location:Mass, z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data=scardat) ## tried to debug this by dumping zicounts.R to a file, modifying ## it to put a trace argument in that would print out the parameters ## and log-likelihood for every call to the log-likelihood function. dump(zicounts,file=zicounts.R) source(zicounts.R) zinb.zc - zicounts(resp=Scars~., x =~Location + Lar + Mass + Lar:Mass + Location:Mass, z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data=scardat,trace=TRUE) ## this actually didn't do any good because the negative log-likelihood ## function never gets called -- as it turns out optim() barfs when it ## gets its initial values, before it ever gets to evaluating the log-likelihood ## check the glm -- this is the equivalent of what zicounts does to ## get the initial values of the x parameters p1 - glm(Scars~Location + Lar + Mass + Lar:Mass + Location:Mass, data=scardat,family=poisson) which(is.na(coef(p1))) ## find out what the deal is table(scardat$Location) scar2 = subset(scardat,Location!=Randalstown) ## first step to removing the bad point from the data set -- but ... table(scar2$Location) ## it leaves the Location factor with the same levels, so ## now we have ZERO counts for one location: ## redefine the factor to drop unused levels scar2$Location - factor(scar2$Location) ## OK, looks fine now table(scar2$Location) zinb.zc - zicounts(resp=Scars~., x =~Location + Lar + Mass + Lar:Mass + Location:Mass, z =~Location + Lar + Mass + Lar:Mass + Location:Mass, data=scar2) ## now we get another error (system is computationally singular when ## trying to compute Hessian -- overparameterized?) Not in any ## trivial way that I can see. It would be nice to get into the guts ## of zicounts and stop it from trying to invert the Hessian, which is ## I think where this happens. In the meanwhile, I have some other ideas about this analysis (sorry, but you started it ...) Looking at the data in a few different ways: library(lattice) xyplot(Scars~Mass,groups=Location,data=scar2,jitter=TRUE, auto.key=list(columns=3)) xyplot(Scars~Mass|Location,data=scar2,jitter=TRUE) xyplot(Scars~Lar,groups=Location,data=scar2, auto.key=list(columns=3)) xyplot(Scars~Mass|Lar,data=scar2) xyplot(Scars~Lar|Location,data=scar2) Some thoughts: (1) I'm not at all sure that zero-inflation is necessary (see Warton 2005, Environmentrics). This is a fairly small, noisy data set without huge numbers of zeros -- a plain old negative binomial might be fine. I don't actually see a lot of signal here, period (although there may be some) ... there's not a huge range in Lar (whatever it is -- the rest of the covariates I think I can interpret). It would be tempting to try to fit location as a random effect, because fitting all those extra degrees of freedom is going to kill you. On the other hand, GLMMs are a bit hairy. cheers Ben __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization question
Dear R users, Imagine please an optimization problem: minimize sum S1+S2 Subject to : y - x = a + S1 x - y = a + S2 and we want to add two more constraints: y - x = b - S3 x - y = b - S4 where a is a small constant value and b is a large constant value, S1 and S2 are surplus and S3 and S4 are slack variables. S3 and S4 have to be maximized, not minimized in objective function. But how to write this? Is this correct? : minimize sum S1+ S2 - S3 -S4 where actually we want to minimize S1 and S2; and maximize S3 and S4. If it is not correct, how to formulate this? what to do ? Thank you for any guide. Tobias - Pinpoint customers who are looking for what you sell. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
On 7/16/07, massimiliano.talarico [EMAIL PROTECTED] wrote: I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? What do you mean by 'radq', Massimiliano? Paul __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
I'm sorry the function is sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; Have you any suggests. Thanks, Massimiliano What is radq? --- massimiliano.talarico [EMAIL PROTECTED] wrote: Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
G'day Massimiliano, On Mon, 16 Jul 2007 22:49:32 +0200 massimiliano.talarico [EMAIL PROTECTED] wrote: Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Note that given the constraint x1+x2+x3=1 and the positivity constraints on x1, x2 and x3, the conditions that all of them should be =1 are redundant. I would go about as follows: x1+x2+x3=1 means x1=1-x2-x3 plug this into the next constraint to obtain: a*(1-x2-x3)^2+b*x2^2+c*x3^2=0.04^2 for suitable a, b, c. Note that for any value of x3, you can solve this equation for x2. Hence, take a grid of x3 values between 0 and 1 and calculate the correpsonding x2 values. From x3 and x2 calculate x1. Discard all tuples that do not fulfill the positivity constraints and calculate the function values for the remaining ones. If the grid of x3 values is fine enough, that should give you an idea what the solution should be. Alternatively, you could write a function that takes values of x3 and does all the computations outline above (return -Inf if the value x3 is not feasible) and then pass the function to optimise() for numerical optimisation... Hope this helps. Cheers, Berwin === Full address = Berwin A TurlachTel.: +65 6515 4416 (secr) Dept of Statistics and Applied Probability+65 6515 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED] Singapore 117546http://www.stat.nus.edu.sg/~statba __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
My apologies, didn't see the boundary constraints. Try this one... f - function(x) (sqrt((x[1]*0.114434)^2+(x[2]*0.043966)^2+(x[3]*0.100031)^2)-0.04)^2 optim(par=rep(0,3),f,lower=rep(0,3),upper=rep(1,3),method=L-BFGS-B) and check ?optim --- massimiliano.talarico [EMAIL PROTECTED] wrote: I'm sorry the function is sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; Have you any suggests. Thanks, Massimiliano What is radq? --- massimiliano.talarico [EMAIL PROTECTED] wrote: Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
f - function(x) (sqrt((x[1]*0.114434)^2+(x[2]*0.043966)^2+(x[3]*0.100031)^2)-0.04)^2 optim(c(0,0,0),f) see ?optim for details on arguments, options, etc. --- massimiliano.talarico [EMAIL PROTECTED] wrote: I'm sorry the function is sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; Have you any suggests. Thanks, Massimiliano What is radq? --- massimiliano.talarico [EMAIL PROTECTED] wrote: Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization (MAX) with R
Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano Questo messaggio di posta elettronica contiene informazioni di carattere confidenziale rivolte esclusivamente al destinatario sopra indicato. E' vietato l'uso, la diffusione, distribuzione o riproduzione da parte di ogni altra persona. Nel caso aveste ricevuto questo messaggio di posta elettronica per errore, siete pregati di segnalarlo immediatamente al mittente e distruggere quanto ricevuto (compresi i file allegati) senza farne copia. Qualsivoglia utilizzo non autorizzato del contenuto di questo messaggio costituisce violazione dell'obbligo di non prendere cognizione della corrispondenza tra altri soggetti, salvo piu grave illecito, ed espone il responsabile alle relative conseguenze. This e-mail is confidential and may also contain privileged information. If you are not the intended recipient you are not authorised to read, print, save, process or disclose this message. If you have received this message by mistake, please inform the sender immediately and delete this e-mail, its attachments and any copies. Any use, distribution, reproduction or disclosure by any person other than the intended recipient is strictly prohibited and the person responsible may incur penalties. Thank you! [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
Thanks for your suggests, but I need to obtain the MAX of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x2=0; x3=0; Thanks and again Thanks, Massimiliano My apologies, didn't see the boundary constraints. Try this one... f - function(x) (sqrt((x[1]*0.114434)^2+(x[2]*0.043966)^2+(x[3]*0.100031) ^2)-0.04)^2 optim(par=rep(0,3),f,lower=rep(0,3),upper=rep (1,3),method=L-BFGS-B) and check ?optim --- massimiliano.talarico [EMAIL PROTECTED] wrote: I'm sorry the function is sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2) =0.04; Have you any suggests. Thanks, Massimiliano What is radq? --- massimiliano.talarico [EMAIL PROTECTED] wrote: Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2) =0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Fussy? Opinionated? Impossible to please? Perfect. Join Yahoo!'s user panel and lay it on us. http://surveylink.yahoo.com/gmrs/yahoo_panel_invite.asp?a=7 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
My apologies, I read the post over too quickly (even the second time). It's been a while since I've played around with anything other than box constraints, but this one is conducive to a brute-force approach (employing Berwin suggestions). The pseudo-code would look something like this: delta - 1e-3 # grid space of x3, the smaller the better oldvalue - -Inf # some initial value for objective function for( x3 in seq(0,1,by=delta) ) { ## calculate x1,x2 as per Berwin's response ## if all constraints are met, feasible - TRUE ## else feasible - FALSE if( !feasible ) next # if not feasible, go to next x3 value ## newvalue - value of objective function with x1,x2,x3 if( newvalue oldvalue ) { oldvalue - newvalue max.x1 - x1; max.x2 - x2; max.x3 - x3 } } You should end up with the desired values of max.x1, max.x2, max.x3. Hope this helps, ST --- massimiliano.talarico [EMAIL PROTECTED] wrote: Thanks for your suggests, but I need to obtain the MAX of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x2=0; x3=0; Thanks and again Thanks, Massimiliano My apologies, didn't see the boundary constraints. Try this one... f - function(x) (sqrt((x[1]*0.114434)^2+(x[2]*0.043966)^2+(x[3]*0.100031) ^2)-0.04)^2 optim(par=rep(0,3),f,lower=rep(0,3),upper=rep (1,3),method=L-BFGS-B) and check ?optim --- massimiliano.talarico [EMAIL PROTECTED] wrote: I'm sorry the function is sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2) =0.04; Have you any suggests. Thanks, Massimiliano What is radq? --- massimiliano.talarico [EMAIL PROTECTED] wrote: Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2) =0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Fussy? Opinionated? Impossible to please? Perfect. Join Yahoo!'s user panel and lay it on us. http://surveylink.yahoo.com/gmrs/yahoo_panel_invite.asp?a=7 that gives answers, not web links. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
Hi, Your problem can be solved analytically. Eliminate one of the variables, say x3, from the problem by using the equality x1 + x2 + x3 = 1. Then solve for the intersection of the circle (in x1 and x2) defined by the radical constraint, with the straight line defined by the objective function. There will be, at most, two intersection points. The extremum has to be one of these two points, provided they also satisfy the other inequalities (To me, this sounds an awful lot like a homework problem). Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of massimiliano.talarico Sent: Monday, July 16, 2007 4:50 PM To: r-help Subject: [R] Optimization Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
This is partially true since both the function to be maximized and the constraint are non-linear. One may substitute 1-x1-x2 for x3 and use (let say) Lagrange multipliers to get two non-linear equations with 2 unknowns for which there should be a function solving them. Then you must find the points where the constraint function intersects with the triangle {x1=0,x2=0,x1+x2=1}, which is easier (for each of the 3 edges you get a non-linear equation in one variable). So even though an (almost) analytical solution can be found it would be much more convenient to use an optimization function which (hopefully) does all this for you. Moshe. --- Ravi Varadhan [EMAIL PROTECTED] wrote: Hi, Your problem can be solved analytically. Eliminate one of the variables, say x3, from the problem by using the equality x1 + x2 + x3 = 1. Then solve for the intersection of the circle (in x1 and x2) defined by the radical constraint, with the straight line defined by the objective function. There will be, at most, two intersection points. The extremum has to be one of these two points, provided they also satisfy the other inequalities (To me, this sounds an awful lot like a homework problem). Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of massimiliano.talarico Sent: Monday, July 16, 2007 4:50 PM To: r-help Subject: [R] Optimization Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
G'day Moshe, On Tue, 17 Jul 2007 17:32:52 -0700 (PDT) Moshe Olshansky [EMAIL PROTECTED] wrote: This is partially true since both the function to be maximized and the constraint are non-linear. I am not sure what your definition of non-linear is, but in my book, and I believe by most mathematical/statistical definitions, the objective function is linear. The only non-linearity comes in through the second constraint. One may substitute 1-x1-x2 for x3 and use (let say) Lagrange multipliers to get two non-linear equations with 2 unknowns for which there should be a function solving them. Why would you want to use Lagrange multipliers? Isn't that a bit of an overkill? Once you substitute 1-x1-x2 for x3 in the second constraint, you have a quadratic equations in x1 and x2. So for any given value of x1 you can solve for x2 (or for any given value of x2 you can solve for x1). They still teach how to solve quadratic equations at school, don't they? ;-) Then you must find the points where the constraint function intersects with the triangle {x1=0,x2=0,x1+x2=1}, which is easier (for each of the 3 edges you get a non-linear equation in one variable). Even easier. Take an x1 between 0 and 1. If for that x1 the quadratic equation in x2 has no real solution, then x1 is not feasible. Otherwise find the values of x2 that solve the equation. Use each of these values together with x1 to calculate corresponding values of x3. Then check these tuples for feasibility. If they are feasible, evaluate the objective function and return the tuple with the larger function value. All the calculations outlined in the paragraph above are easily implemented in R, e.g. the function polyroot() returns the roots of a polynomial. Cheers, Berwin === Full address = Berwin A TurlachTel.: +65 6515 4416 (secr) Dept of Statistics and Applied Probability+65 6515 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED] Singapore 117546http://www.stat.nus.edu.sg/~statba __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
You are right!!! For some strange reason I substituted ^ (exponentiation) for *, so the problem became Max x1^0.021986+x2^0.000964+x3^0.02913 with these conditions: x1+x2+x3=1; sqrt((x1^0.114434)^2+(x2^0.043966)^2+(x3^0.100031)^2)=0.04; which is clearly non-linear. --- Berwin A Turlach [EMAIL PROTECTED] wrote: G'day Moshe, On Tue, 17 Jul 2007 17:32:52 -0700 (PDT) Moshe Olshansky [EMAIL PROTECTED] wrote: This is partially true since both the function to be maximized and the constraint are non-linear. I am not sure what your definition of non-linear is, but in my book, and I believe by most mathematical/statistical definitions, the objective function is linear. The only non-linearity comes in through the second constraint. One may substitute 1-x1-x2 for x3 and use (let say) Lagrange multipliers to get two non-linear equations with 2 unknowns for which there should be a function solving them. Why would you want to use Lagrange multipliers? Isn't that a bit of an overkill? Once you substitute 1-x1-x2 for x3 in the second constraint, you have a quadratic equations in x1 and x2. So for any given value of x1 you can solve for x2 (or for any given value of x2 you can solve for x1). They still teach how to solve quadratic equations at school, don't they? ;-) Then you must find the points where the constraint function intersects with the triangle {x1=0,x2=0,x1+x2=1}, which is easier (for each of the 3 edges you get a non-linear equation in one variable). Even easier. Take an x1 between 0 and 1. If for that x1 the quadratic equation in x2 has no real solution, then x1 is not feasible. Otherwise find the values of x2 that solve the equation. Use each of these values together with x1 to calculate corresponding values of x3. Then check these tuples for feasibility. If they are feasible, evaluate the objective function and return the tuple with the larger function value. All the calculations outlined in the paragraph above are easily implemented in R, e.g. the function polyroot() returns the roots of a polynomial. Cheers, Berwin === Full address = Berwin A TurlachTel.: +65 6515 4416 (secr) Dept of Statistics and Applied Probability +65 6515 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: [EMAIL PROTECTED] Singapore 117546 http://www.stat.nus.edu.sg/~statba __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization (MAX) with R
Talarico Massimiliano (UniCredit Xelion Banca) wrote on 07/17/2007 06:00 PM: Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; sqrt((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Use Lagrange multipliers and do it analytically? Ted. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization
Dear all, I need a suggest to obtain the max of this function: Max x1*0.021986+x2*0.000964+x3*0.02913 with these conditions: x1+x2+x3=1; radq((x1*0.114434)^2+(x2*0.043966)^2+(x3*0.100031)^2)=0.04; x1=0; x1=1; x2=0; x2=1; x3=0; x3=1; Any suggests ? Thanks in advanced, Massimiliano __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
It is of great help for your advice. Thanks a lot to you all. livia wrote: Hi, I would like to minimize the value of x1-x2, x2 is a fixed value of 0.01, x1 is the quantile of normal distribution (0.0032,x) with probability of 0.7, and the changing value should be x. Initial value for x is 0.0207. I am using the following codes, but it does not work. fr - function(x) { x1-qnorm(0.7,0.0032,x) x2=0.01 x1-x2 } xsd - optim(0.0207, fr, NULL,method=BFGS) It is the first time I am trying to use optimization. Could anyone give me some advice? -- View this message in context: http://www.nabble.com/Optimization-tf3941212.html#a11196890 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization
Hi, I would like to minimize the value of x1-x2, x2 is a fixed value of 0.01, x1 is the quantile of normal distribution (0.0032,x) with probability of 0.7, and the changing value should be x. Initial value for x is 0.0207. I am using the following codes, but it does not work. fr - function(x) { x1-qnorm(0.7,0.0032,x) x2=0.01 x1-x2 } xsd - optim(0.0207, fr, NULL,method=BFGS) It is the first time I am trying to use optimization. Could anyone give me some advice? -- View this message in context: http://www.nabble.com/Optimization-tf3941212.html#a11178663 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
You don't need optimization for the solution to your problem. You just need an understanding of the meaning of qnorm() and some simple algebra. Try: x- (0.01-0.0032)/qnorm(0.7,0,1) At 12:01 PM 6/18/2007, you wrote: Hi, I would like to minimize the value of x1-x2, x2 is a fixed value of 0.01, x1 is the quantile of normal distribution (0.0032,x) with probability of 0.7, and the changing value should be x. Initial value for x is 0.0207. I am using the following codes, but it does not work. fr - function(x) { x1-qnorm(0.7,0.0032,x) x2=0.01 x1-x2 } xsd - optim(0.0207, fr, NULL,method=BFGS) It is the first time I am trying to use optimization. Could anyone give me some advice? -- View this message in context: http://www.nabble.com/Optimization-tf3941212.html#a11178663 Sent from the R help mailing list archive at Nabble.com. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: [EMAIL PROTECTED] Least Cost Formulations, Ltd.URL: http://lcfltd.com/ 824 Timberlake Drive Tel: 757-467-0954 Virginia Beach, VA 23464-3239Fax: 757-467-2947 Vere scire est per causas scire __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
livia wrote: Hi, I would like to minimize the value of x1-x2, x2 is a fixed value of 0.01, x1 is the quantile of normal distribution (0.0032,x) with probability of 0.7, and the changing value should be x. Initial value for x is 0.0207. I am using the following codes, but it does not work. fr - function(x) { x1-qnorm(0.7,0.0032,x) x2=0.01 x1-x2 } xsd - optim(0.0207, fr, NULL,method=BFGS) I guess you want to use optimize() and change the last line of fr to (x1-x2)^2 as in: fr - function(x) { x1 - qnorm(0.7, 0.0032, x) x2 - 0.01 (x1-x2)^2 } optimize(fr, c(-5, 5)) Uwe Ligges It is the first time I am trying to use optimization. Could anyone give me some advice? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
On 18-Jun-07 16:01:03, livia wrote: Hi, I would like to minimize the value of x1-x2, x2 is a fixed value of 0.01, x1 is the quantile of normal distribution (0.0032,x) with probability of 0.7, and the changing value should be x. Initial value for x is 0.0207. I'm a bit puzzled by the question. If I understand it right, we can ignore x2 (since it is a fixed value) and simply consider minimising x1 (instead of x1-x2). Then, denoting by P(u) the cumulative normal distribution function for mean=0 and variance=1 (i.e. in R: pnorm(u,0,1)), and by Q(p) its inverse, corresponding to qnorm(p,0,1), we have (again if I have understood right): P((x1 - 0.0032)/x) = 0.7 so x1 = 0.0032 + x*Q(0.7) and therefore, since Q(0.7) 0 and x must be positive, the value of x1 can be made as close to 0.032 as you please (but greater than 0.032) by taking x small enough. Hence there is no strictly minimising value of x, but the greatest lower bound of all possible values of x1 is 0.032. Then you can subtract x2. The fact that there is no positive value of x which gives this bound as the value probably explains the failure of your optim() attempt. Best wishes, Ted. E-Mail: (Ted Harding) [EMAIL PROTECTED] Fax-to-email: +44 (0)870 094 0861 Date: 18-Jun-07 Time: 17:46:01 -- XFMail -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
From the help page: Note: 'optim' will work with one-dimensional 'par's, but the default method does not work well (and will warn). Use 'optimize' instead. Next, there is a constraint of x=0 that you are not imposing. Finally, it is easy to see that qnorm(0.7, 0.0032, x) is monotome in x, so the solution is x=0. In fact, x1 = 0.0032 + sqrt(x) * qnorm(0.7). optim(0.0207, fr) does a good enough job, as does optimize(fr, low=0, up=0.05) Advice: numerical optimization is not a black box, and has to be used with some analysis of the problem to hand. See e.g. MASS4, chapter 16. On Mon, 18 Jun 2007, livia wrote: Hi, I would like to minimize the value of x1-x2, x2 is a fixed value of 0.01, x1 is the quantile of normal distribution (0.0032,x) with probability of 0.7, and the changing value should be x. Initial value for x is 0.0207. I am using the following codes, but it does not work. fr - function(x) { x1-qnorm(0.7,0.0032,x) x2=0.01 x1-x2 } xsd - optim(0.0207, fr, NULL,method=BFGS) It is the first time I am trying to use optimization. Could anyone give me some advice? -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization
Hi, my first guess is that the algorithm returns a negative value in some step - recall that you start from 0.0207!! This negative value is then passed as standard error to qnorm and that cannot work... My guess is based on a small experiment where I tried a different starting point (.02 is so close to 0 that one cannot see anything): xsd - optim(20, fr, NULL,method=BFGS,control=list(trace=6)) The warnings which you didn't include also tell you about NaNs in qnorm() - another strong indication of wrong arguments to qnorm(). Try constrained optimization to resctrict to positive values. See ?constrOptim or use optim() with a method allowing for box constraints - see ?optim, arguments lower, upper. Petr livia napsal(a): Hi, I would like to minimize the value of x1-x2, x2 is a fixed value of 0.01, x1 is the quantile of normal distribution (0.0032,x) with probability of 0.7, and the changing value should be x. Initial value for x is 0.0207. I am using the following codes, but it does not work. fr - function(x) { x1-qnorm(0.7,0.0032,x) x2=0.01 x1-x2 } xsd - optim(0.0207, fr, NULL,method=BFGS) It is the first time I am trying to use optimization. Could anyone give me some advice? -- Petr Klasterecky Dept. of Probability and Statistics Charles University in Prague Czech Republic __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization and simulation
Dear all, I would need to maximize a self-defined 'target' function(see below) with respect to theta, where v follows a log-normal distribution with mean 'mu(x)' and a constant variance. For each v drawn from its distribution, one maximized value and optimal theta are produced. I'd like to do this repeatedly and store the maximized value and corresponding theta. I wrote the following code that can produce a result. But the problem is that the result doesn't seem to be the optimized one because if I arbitrarily choose theta=c(4,27), I will get much bigger value than those simulated ones. I am not sure where is wrong. Could anyone help me with this? Thank you in advance! Here is the code: rdt-matrix(0,10,3) for (i in 1:10){ #Define the target function target-function(theta){ d0=theta[1] T0=theta[2] mu-function(x,d=d0,ta=23.86,ti=11.067,ppt=1.321,a=-12.31, S=5.62, L=3.83, b=c(0.338,-0.0055,0.113,-0.00466,-0.008,-0.205,-0.044,0.266,1.719,-0.169)){ return(a+sum(b*c(x,x^2,d,d^2,S,L,ta,ti,ppt,ti*ppt)))} v-function(x,d=d0){ return(rlnorm(1,mean=mu(x),sd=0.835^0.5))} p-function(x,d=d0){ if(8.179+0.00023*(5.29+0.16*x-0.21*d)^545){return(45)} return(8.179+0.00023*(5.29+0.16*x-0.21*d)^5) } Y3-function(x,d=d0){ return(p(x,d)*v(x,d)) } r-Y3(T0) return(r) } result-optim(c(2,10),target,lower=c(0.1,0.5),upper=c(66,50), method=L-BFGS-B,control=list(maxit=10,fnscale=-1)) rdt[i,1]-result$value rdt[i,2:3]-result$par } - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization and garch
Good day, Here I was trying to write a code for Garch(1,1) . As garch problem is more or less an optimization problem I also tried to get the algorithm for nlminb function. What I saw that if use this function 'nlminb I can easyly get the estimate of parameters. But any other function is not working. I tried to write my own code for optimization using Quasi-Newton Methods etc but although it is working for ordinary non-linear function, it fails in garch case. Therefore I am trying to get a step by step documentation for nlminb function. I already gone though its help page got a look on http://netlib.bell-labs.com/cm/cs/cstr/153.pdf. But it did not solve my problem. In this regards, can anyone give me any step-by-step approach or theory behind the calculation that 'nlminb uses? Any help will be highly appreciable. Thanks and regards, __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Optimization and garch
Optimizing GARCH likelihoods is notoriously difficult. I suspect that you will find 'nlminb' to be less than perfect, though it is relatively good. In particular you are likely to see different behavior depending on whether or not the data are in percent. A reference is Winker and Maringer (2006) The convergence of optimization based estimators: Theory and application to a GARCH-model. Available online. Patrick Burns [EMAIL PROTECTED] +44 (0)20 8525 0696 http://www.burns-stat.com (home of S Poetry and A Guide for the Unwilling S User) stat stat wrote: Good day, Here I was trying to write a code for Garch(1,1) . As garch problem is more or less an optimization problem I also tried to get the algorithm for nlminb function. What I saw that if use this function 'nlminb I can easyly get the estimate of parameters. But any other function is not working. I tried to write my own code for optimization using Quasi-Newton Methods etc but although it is working for ordinary non-linear function, it fails in garch case. Therefore I am trying to get a step by step documentation for nlminb function. I already gone though its help page got a look on http://netlib.bell-labs.com/cm/cs/cstr/153.pdf. But it did not solve my problem. In this regards, can anyone give me any step-by-step approach or theory behind the calculation that 'nlminb uses? Any help will be highly appreciable. Thanks and regards, __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] R: Optimization
Residual standard error: 0.5766 on 104 degrees of freedom Multiple R-Squared: 0.523, Adjusted R-squared: 0.5001 F-statistic: 22.8 on 5 and 104 DF, p-value: 2.179e-15 summary(fit.1) Call: lm(formula = log((Yield + 1)/(Sedi + 1)) ~ log((Sal + 1)/(Sedi + 1)) + log((Bathy + 1)/(Sedi + 1)) + log((Chl + 1)/(Sedi + 1)) + log((Hydro + 1)/(Sedi + 1)) + log((Oxy + 1)/(Sedi + 1)), data = subset(data.df)) Residuals: Min 1Q Median 3Q Max -1.05552 -0.23441 0.01263 0.18355 1.24473 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 1.849160.15474 11.950 2e-16 *** log((Sal + 1)/(Sedi + 1))5.038631.00977 4.990 2.46e-06 *** log((Bathy + 1)/(Sedi + 1)) 0.052790.20767 0.254 0.7999 log((Chl + 1)/(Sedi + 1)) -10.966822.27270 -4.825 4.86e-06 *** log((Hydro + 1)/(Sedi + 1)) 1.746310.28980 6.026 2.64e-08 *** log((Oxy + 1)/(Sedi + 1))3.959381.98206 1.998 0.0484 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3751 on 103 degrees of freedom Multiple R-Squared: 0.5606, Adjusted R-squared: 0.5393 F-statistic: 26.28 on 5 and 103 DF, p-value: 2.2e-16 anova(fit1,fit0) Analysis of Variance Table Res.Df RSS Df Sum of Sq F Pr(F) 1102 14.2409 2103 14.3650 -1 -0.1241 0.8891 0.3479 Thanks for the help Simone Vincenzi _ Simone Vincenzi, PhD Student Department of Environmental Sciences University of Parma Parco Area delle Scienze, 33/A, 43100 Parma, Italy Phone: +39 0521 905696 Fax: +39 0521 906611 e.mail: [EMAIL PROTECTED] -Messaggio originale- Da: Spencer Graves [mailto:[EMAIL PROTECTED] Inviato: lunedì 4 settembre 2006 21.31 A: Simone Vincenzi Cc: r-help@stat.math.ethz.ch Oggetto: Re: [R] Optimization Have you considered talking logarithms of the expression you mentioned: log(Yield) = a1*log(A)+b1*log(B)+c2*log(C)+... where a1 = a/(a+b+...), etc. This model has two constraints not present in ordinary least squares: First, the intercept is assumed to be zero. Second, the coefficients in this log formulation must sum to 1. If I were you, I might use something like lm to test them both. To explain how, I'll modify the notation, replacing A by X1, B by X2, ..., up to Xkm1 (= X[k-1]) and Xk for k different environmental variables. Then I might try something like the following: fit0 - lm(log(Yield) ~ log(X1) + ... + log(Xk)-1 ) fit1 - lm(log(Yield) ~ log(X1) + ... + log(Xk) ) fit.1 - lm(log(Yield/Xk) ~ log(X1/Xk) + ... + log(Xkm1/Xk) ) fit.0 - lm(log(Yield/Xk) ~ log(X1/Xk) + ... + log(Xkm1/Xk)-1 ) anova(fit1, fit0) would test the no-constant model, and if I haven't made a mistake in this, anova(fit0, fit.0) and anova(fit1, fit.1) would test the constraint that all the coefficients should sum to 1. If you would like further help from this listserve, please provide commented, minimal, self-contained, reproducible code to help potential respondents understand your question and concerns (as suggested in the posting guide www.R-project.org/posting-guide.html). Hope this helps. Spencer Graves Simone Vincenzi wrote: Dear R-list, I'm trying to estimate the relative importance of 6 environmental variables in determining clam yield. To estimate clam yield a previous work used the function Yield = (A^a*B^b*C^c...)^1/(a+b+c+...) where A,B,C... are the values of the environmental variables and the weights a,b,c... have not been calibrated on data but taken from literature. Now I'd like to estimate the weights a,b,c... by using a dataset with 110 observations of yield and values of the environmental variables. I'm wondering if it is feasible or if the number of observation is too low, if some data transformation is needed and which R function is the most appropriate to try to estimate the weights. Any help would be greatly appreciated. Simone Vincenzi _ Simone Vincenzi, PhD Student Department of Environmental Sciences University of Parma Parco Area delle Scienze, 33/A, 43100 Parma, Italy Phone: +39 0521 905696 Fax: +39 0521 906611 e.mail: [EMAIL PROTECTED] -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- No virus found in this incoming message. -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R: Optimization
= data.df) Residuals: Min 1Q Median 3Q Max -1.05485 -0.23759 0.01331 0.18692 1.23803 Coefficients: Estimate Std. Error t value Pr(|t|) log(Sal + 1) 5.071931.00584 5.042 1.98e-06 *** log(Bathy + 1) 0.058040.20684 0.281 0.779561 log(Chl + 1) -8.539412.11720 -4.033 0.000106 *** log(Hydro + 1) 1.738350.28815 6.033 2.56e-08 *** log(Oxy + 1) 4.199511.98459 2.116 0.036750 * log(Sedi + 1) 1.159530.14807 7.831 4.51e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3735 on 103 degrees of freedom Multiple R-Squared: 0.9148, Adjusted R-squared: 0.9098 F-statistic: 184.2 on 6 and 103 DF, p-value: 2.2e-16 summary(fit1) Call: lm(formula = log(Yield + 1) ~ log(Sal + 1) + log(Bathy + 1) + log(Chl + 1) + log(Hydro + 1) + log(Oxy + 1) + log(Sedi + 1), data = (data.df)) Residuals: Min1QMedian3Q Max -1.057766 -0.227720 0.006146 0.192200 1.222543 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) -4.494004.76603 -0.943 0.3479 log(Sal + 1) 5.134991.00861 5.091 1.63e-06 *** log(Bathy + 1) 0.066850.20716 0.323 0.7476 log(Chl + 1) -2.489096.75718 -0.368 0.7134 log(Hydro + 1) 1.706740.29025 5.880 5.22e-08 *** log(Oxy + 1) 4.637582.03928 2.274 0.0251 * log(Sedi + 1) 1.140050.14959 7.621 1.34e-11 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3737 on 102 degrees of freedom Multiple R-Squared: 0.7589, Adjusted R-squared: 0.7447 F-statistic: 53.51 on 6 and 102 DF, p-value: 2.2e-16 summary(fit.0) Call: lm(formula = log((Yield + 1)/(Sedi + 1)) ~ log((Sal + 1)/(Sedi + 1)) + log((Bathy + 1)/(Sedi + 1)) + log((Chl + 1)/(Sedi + 1)) + log((Hydro + 1)/(Sedi + 1)) + log((Oxy + 1)/(Sedi + 1)) - 1, data = subset(data.df)) Residuals: Min 1Q Median 3Q Max -1.5762 -0.2591 0.1065 0.5023 1.4533 Coefficients: Estimate Std. Error t value Pr(|t|) log((Sal + 1)/(Sedi + 1))3.0170 1.5304 1.971 0.05134 . log((Bathy + 1)/(Sedi + 1)) -0.3982 0.3139 -1.268 0.20748 log((Chl + 1)/(Sedi + 1))8.8137 2.3941 3.681 0.00037 *** log((Hydro + 1)/(Sedi + 1)) 0.3309 0.4066 0.814 0.41760 log((Oxy + 1)/(Sedi + 1)) -12.5242 2.1881 -5.724 1.01e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.5766 on 104 degrees of freedom Multiple R-Squared: 0.523, Adjusted R-squared: 0.5001 F-statistic: 22.8 on 5 and 104 DF, p-value: 2.179e-15 summary(fit.1) Call: lm(formula = log((Yield + 1)/(Sedi + 1)) ~ log((Sal + 1)/(Sedi + 1)) + log((Bathy + 1)/(Sedi + 1)) + log((Chl + 1)/(Sedi + 1)) + log((Hydro + 1)/(Sedi + 1)) + log((Oxy + 1)/(Sedi + 1)), data = subset(data.df)) Residuals: Min 1Q Median 3Q Max -1.05552 -0.23441 0.01263 0.18355 1.24473 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 1.849160.15474 11.950 2e-16 *** log((Sal + 1)/(Sedi + 1))5.038631.00977 4.990 2.46e-06 *** log((Bathy + 1)/(Sedi + 1)) 0.052790.20767 0.254 0.7999 log((Chl + 1)/(Sedi + 1)) -10.966822.27270 -4.825 4.86e-06 *** log((Hydro + 1)/(Sedi + 1)) 1.746310.28980 6.026 2.64e-08 *** log((Oxy + 1)/(Sedi + 1))3.959381.98206 1.998 0.0484 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3751 on 103 degrees of freedom Multiple R-Squared: 0.5606, Adjusted R-squared: 0.5393 F-statistic: 26.28 on 5 and 103 DF, p-value: 2.2e-16 anova(fit1,fit0) Analysis of Variance Table Res.Df RSS Df Sum of Sq F Pr(F) 1102 14.2409 2103 14.3650 -1 -0.1241 0.8891 0.3479 Thanks for the help Simone Vincenzi _ Simone Vincenzi, PhD Student Department of Environmental Sciences University of Parma Parco Area delle Scienze, 33/A, 43100 Parma, Italy Phone: +39 0521 905696 Fax: +39 0521 906611 e.mail: [EMAIL PROTECTED] -Messaggio originale- Da: Spencer Graves [mailto:[EMAIL PROTECTED] Inviato: lunedì 4 settembre 2006 21.31 A: Simone Vincenzi Cc: r-help@stat.math.ethz.ch Oggetto: Re: [R] Optimization Have you considered talking logarithms of the expression you mentioned: log(Yield) = a1*log(A)+b1*log(B)+c2*log(C)+... where a1 = a/(a+b+...), etc. This model has two constraints not present in ordinary least squares: First, the intercept is assumed to be zero. Second, the coefficients in this log formulation must
Re: [R] Optimization
Have you considered talking logarithms of the expression you mentioned: log(Yield) = a1*log(A)+b1*log(B)+c2*log(C)+... where a1 = a/(a+b+...), etc. This model has two constraints not present in ordinary least squares: First, the intercept is assumed to be zero. Second, the coefficients in this log formulation must sum to 1. If I were you, I might use something like lm to test them both. To explain how, I'll modify the notation, replacing A by X1, B by X2, ..., up to Xkm1 (= X[k-1]) and Xk for k different environmental variables. Then I might try something like the following: fit0 - lm(log(Yield) ~ log(X1) + ... + log(Xk)-1 ) fit1 - lm(log(Yield) ~ log(X1) + ... + log(Xk) ) fit.1 - lm(log(Yield/Xk) ~ log(X1/Xk) + ... + log(Xkm1/Xk) ) fit.0 - lm(log(Yield/Xk) ~ log(X1/Xk) + ... + log(Xkm1/Xk)-1 ) anova(fit1, fit0) would test the no-constant model, and if I haven't made a mistake in this, anova(fit0, fit.0) and anova(fit1, fit.1) would test the constraint that all the coefficients should sum to 1. If you would like further help from this listserve, please provide commented, minimal, self-contained, reproducible code to help potential respondents understand your question and concerns (as suggested in the posting guide www.R-project.org/posting-guide.html). Hope this helps. Spencer Graves Simone Vincenzi wrote: Dear R-list, I'm trying to estimate the relative importance of 6 environmental variables in determining clam yield. To estimate clam yield a previous work used the function Yield = (A^a*B^b*C^c...)^1/(a+b+c+...) where A,B,C... are the values of the environmental variables and the weights a,b,c... have not been calibrated on data but taken from literature. Now I'd like to estimate the weights a,b,c... by using a dataset with 110 observations of yield and values of the environmental variables. I'm wondering if it is feasible or if the number of observation is too low, if some data transformation is needed and which R function is the most appropriate to try to estimate the weights. Any help would be greatly appreciated. Simone Vincenzi _ Simone Vincenzi, PhD Student Department of Environmental Sciences University of Parma Parco Area delle Scienze, 33/A, 43100 Parma, Italy Phone: +39 0521 905696 Fax: +39 0521 906611 e.mail: [EMAIL PROTECTED] -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization
Dear R-list, I'm trying to estimate the relative importance of 6 environmental variables in determining clam yield. To estimate clam yield a previous work used the function Yield = (A^a*B^b*C^c...)^1/(a+b+c+...) where A,B,C... are the values of the environmental variables and the weights a,b,c... have not been calibrated on data but taken from literature. Now I'd like to estimate the weights a,b,c... by using a dataset with 110 observations of yield and values of the environmental variables. I'm wondering if it is feasible or if the number of observation is too low, if some data transformation is needed and which R function is the most appropriate to try to estimate the weights. Any help would be greatly appreciated. Simone Vincenzi _ Simone Vincenzi, PhD Student Department of Environmental Sciences University of Parma Parco Area delle Scienze, 33/A, 43100 Parma, Italy Phone: +39 0521 905696 Fax: +39 0521 906611 e.mail: [EMAIL PROTECTED] -- __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Optimization problem: selecting independent rows to maximizethe mean
Does R have packages for such multi-objectives optimization problems ? The rgenoud (R-GENetic Optimization Using Derivatives) package allows for multiple object optimization problems. See the lexical option which searches for the Pareto front. The package is written for NP-hard problems (but they are...well...difficult). See CRAN or: http://sekhon.berkeley.edu/rgenoud/ Cheers, Jas. === Jasjeet S. Sekhon Associate Professor Survey Research Center UC Berkeley http://sekhon.berkeley.edu/ V: 510-642-9974 F: 617-507-5524 === nojhan wrote: Le Wed, 01 Mar 2006 13:07:07 -0800, Berton Gunter a ?crit : 2) That the mean and sd can be simultaneously optimized as you describe-- what if the subset with maximum mean also has bigger than minimal sd? Then you have two choices : 1) balance the two objectives with weights, according to the importance you give to each one 2) get a list of non-dominated solutions (a Pareto front) Does R have packages for such multi-objectives optimization problems ? Moreover, does it have a package for difficult (i.e. NP-hard) problems ? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Optimization problem: selecting independent rows to maximizethe mean
Regarding multi-object optimization, I just got 0 hits from RSiteSearch(multi-objective optimization) and RSiteSearch(multiobjective optimization). However, it shouldn't be too difficult to write a wrapper function to blend other functions however you would like, then use optim or nlminb or one of the other optimizers in R. I don't feel qualified to even comment on 'difficult (i.e. NP-hard) problems'. hope this helps, spencer graves nojhan wrote: Le Wed, 01 Mar 2006 13:07:07 -0800, Berton Gunter a écrit : 2) That the mean and sd can be simultaneously optimized as you describe-- what if the subset with maximum mean also has bigger than minimal sd? Then you have two choices : 1) balance the two objectives with weights, according to the importance you give to each one 2) get a list of non-dominated solutions (a Pareto front) Does R have packages for such multi-objectives optimization problems ? Moreover, does it have a package for difficult (i.e. NP-hard) problems ? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Optimization problem: selecting independent rows to maximizethe mean
Le Wed, 01 Mar 2006 13:07:07 -0800, Berton Gunter a écrit : 2) That the mean and sd can be simultaneously optimized as you describe-- what if the subset with maximum mean also has bigger than minimal sd? Then you have two choices : 1) balance the two objectives with weights, according to the importance you give to each one 2) get a list of non-dominated solutions (a Pareto front) Does R have packages for such multi-objectives optimization problems ? Moreover, does it have a package for difficult (i.e. NP-hard) problems ? -- nojhan __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Optimization problem: selecting independent rows to maximize the mean
Dear R community, I have a dataframe with 500,000 rows and 102 columns. The rows represent spatial polygons, some of which overlap others (i.e., not all rows are independent of each other). Given a particular row, the first column contains a unique RowID. The second column contains the Variable of interest. The remaining 100 columns (Overlap1 ... Overlap100) each contain a row ID that overlaps this row (but if this row overlaps fewer than 100 other rows then the remainder of the columns OL1...OL100 contain NA). Here's the problem: I need to select the subset of 500 independent rows that maximizes the mean and minimizes the stdev of Variable. Clearly this requires iterative selection and comparison of rows, because each newly-selected row must be compared to rows already selected to ensure it does not overlap them. At each step, a row already selected might be removed from the subset if it can be replaced with another that increases the mean and/or reduces the stdev. The above description is a simplification of my problem, but it's a start. As I am new to R (and programming in general) I'm not sure how to start thinking about this, or even where to look. I'd appreciate any ideas that might help. Thank you, Mark __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Optimization problem: selecting independent rows to maximizethe mean
This sounds either easy via a greedy algorithm or NP-hard. Moreover, it is not clear to me that 1) A subset of 500 indpendent rows exists, where I presume independent means pairwise nonoverlapping; 2) That the mean and sd can be simultaneously optimized as you describe-- what if the subset with maximum mean also has bigger than minimal sd? -- Bert Gunter Genentech Non-Clinical Statistics South San Francisco, CA The business of the statistician is to catalyze the scientific learning process. - George E. P. Box -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Mark Sent: Wednesday, March 01, 2006 12:40 PM To: r-help@stat.math.ethz.ch Subject: [R] Optimization problem: selecting independent rows to maximizethe mean Dear R community, I have a dataframe with 500,000 rows and 102 columns. The rows represent spatial polygons, some of which overlap others (i.e., not all rows are independent of each other). Given a particular row, the first column contains a unique RowID. The second column contains the Variable of interest. The remaining 100 columns (Overlap1 ... Overlap100) each contain a row ID that overlaps this row (but if this row overlaps fewer than 100 other rows then the remainder of the columns OL1...OL100 contain NA). Here's the problem: I need to select the subset of 500 independent rows that maximizes the mean and minimizes the stdev of Variable. Clearly this requires iterative selection and comparison of rows, because each newly-selected row must be compared to rows already selected to ensure it does not overlap them. At each step, a row already selected might be removed from the subset if it can be replaced with another that increases the mean and/or reduces the stdev. The above description is a simplification of my problem, but it's a start. As I am new to R (and programming in general) I'm not sure how to start thinking about this, or even where to look. I'd appreciate any ideas that might help. Thank you, Mark __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Optimization problem: selecting independent rows to maximize the mean
Package lpSolve might help. On 3/1/06, Mark [EMAIL PROTECTED] wrote: Dear R community, I have a dataframe with 500,000 rows and 102 columns. The rows represent spatial polygons, some of which overlap others (i.e., not all rows are independent of each other). Given a particular row, the first column contains a unique RowID. The second column contains the Variable of interest. The remaining 100 columns (Overlap1 ... Overlap100) each contain a row ID that overlaps this row (but if this row overlaps fewer than 100 other rows then the remainder of the columns OL1...OL100 contain NA). Here's the problem: I need to select the subset of 500 independent rows that maximizes the mean and minimizes the stdev of Variable. Clearly this requires iterative selection and comparison of rows, because each newly-selected row must be compared to rows already selected to ensure it does not overlap them. At each step, a row already selected might be removed from the subset if it can be replaced with another that increases the mean and/or reduces the stdev. The above description is a simplification of my problem, but it's a start. As I am new to R (and programming in general) I'm not sure how to start thinking about this, or even where to look. I'd appreciate any ideas that might help. Thank you, Mark __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] optimization with inequalities
I have to estimate the following model for several group of observations : y(1-y) = p[1]*(x^2-y) + p[2]*y*(x-1) + p[3]*(x-y) with constraints : p[1]+p[3] = 1 p[1]+p[2]+p[3]+1 = 0 p[3] = 0 I use the following code : func - sum((y(1-y) - p[1]*(x^2-y) + p[2]*y*(x-1) + p[3]*(x-y))^2) estim - optim( c(1,0,0),func, method=L-BFGS-B , lower=c(1-p[3], -p[1]-p[3]-1, 0) ) and for some group of observations, I observe that the estimated parameters don't respect the constraints, espacially the first. Where's the problem please ? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] optimization with inequalities
Florent Bresson [EMAIL PROTECTED] writes: I have to estimate the following model for several group of observations : y(1-y) = p[1]*(x^2-y) + p[2]*y*(x-1) + p[3]*(x-y) with constraints : p[1]+p[3] = 1 p[1]+p[2]+p[3]+1 = 0 p[3] = 0 I use the following code : func - sum((y(1-y) - p[1]*(x^2-y) + p[2]*y*(x-1) + p[3]*(x-y))^2) estim - optim( c(1,0,0),func, method=L-BFGS-B , lower=c(1-p[3], -p[1]-p[3]-1, 0) ) and for some group of observations, I observe that the estimated parameters don't respect the constraints, espacially the first. Where's the problem please ? If you think the boundaries in lower=c() are recomputed as the iteration progresses, you're wrong. L-BGFS-B does box constraints only. Instead parametrize using q1=p1+p3 q2=p1+p2+p3 q3=p3 which is easily inverted to get the p's from the q's. Then optimize as a function of q1..q3, substituting the inversion in the expression for func (which btw needs to be a _function_), using the relevant box constraints. -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] optimization with inequalities
On Mon, 28 Nov 2005, Florent Bresson wrote: I have to estimate the following model for several group of observations : y(1-y) = p[1]*(x^2-y) + p[2]*y*(x-1) + p[3]*(x-y) with constraints : p[1]+p[3] = 1 p[1]+p[2]+p[3]+1 = 0 p[3] = 0 I use the following code : func - sum((y(1-y) - p[1]*(x^2-y) + p[2]*y*(x-1) + p[3]*(x-y))^2) estim - optim( c(1,0,0),func, method=L-BFGS-B , lower=c(1-p[3], -p[1]-p[3]-1, 0) ) and for some group of observations, I observe that the estimated parameters don't respect the constraints, espacially the first. Where's the problem please ? User mis-reading the help page! L-BFGS-B handles `box constraints', not linear inequality constraints. You can reparametrize to make these box constraints: use p[3], p[1]+p[3] and p[1]+p[2]+p[3] are variables. -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] optimization with inequalities
If I understand this correctly the variables over which you are optimizing are p[1], p[2] and p[3] whereas x and y are fixed and known during the optimization. In that case its a linear programming problem and you could use the lpSolve library which would allow the explicit specification of the constraints. On 11/28/05, Florent Bresson [EMAIL PROTECTED] wrote: I have to estimate the following model for several group of observations : y(1-y) = p[1]*(x^2-y) + p[2]*y*(x-1) + p[3]*(x-y) with constraints : p[1]+p[3] = 1 p[1]+p[2]+p[3]+1 = 0 p[3] = 0 I use the following code : func - sum((y(1-y) - p[1]*(x^2-y) + p[2]*y*(x-1) + p[3]*(x-y))^2) estim - optim( c(1,0,0),func, method=L-BFGS-B , lower=c(1-p[3], -p[1]-p[3]-1, 0) ) and for some group of observations, I observe that the estimated parameters don't respect the constraints, espacially the first. Where's the problem please ? __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] optimization problem in R ... can this be done?
Gregory Gentlemen wrote: Spencer: Thank you for the helpful suggestions. I have another question following some code I wrote. The function below gives a crude approximation for the x of interest (that value of x such that g(x,n) is less than 0 for all n). # // btilda optimize g(n,x) for some fixed x, and then approximately finds that g(n,x) #such that abs(g(n*,x)=0 // btilda - function(range,len) { # range: over which to look for x bb - seq(range[1],range[2],length=len) OBJ - sapply(bb,function(x) {fixed - c(x,100,692,50,1600,v1,217227); return(optimize(g,c(1,1000),maximum=TRUE,tol=0.001,x=fixed)$objective)}) tt - data.frame(b=bb,obj=OBJ) tt$absobj - abs(tt$obj) d - tt[order(tt$absobj),][1:3,] return(as.vector(d)) } For instance a run of btilda(c(20.55806,20.55816),1) returns: b objabsobj 5834 20.55812 -0.00049428480.0004942848 5833 20.55812 0.0011715433 0.0011715433 5835 20.55812-0.00216011400.0021601140 My question is how to improve the precision of b (which is x) here. It seems that seq(20.55806,20.55816,length=1 ) is only precise to 5 or so digits, and thus, is equivalent for numerous succesive Why do you think so? It is much more accurate! See ?.Machine Uwe Ligges values. How can I get around this? Spencer Graves [EMAIL PROTECTED] wrote: Part of the R culture is a statement by Simon Blomberg immortalized in library(fortunes) as, This is R. There is no if. Only how. I can't see now how I would automate a complete solution to your problem in general. However, given a specific g(x, n), I would start by writing a function to use expand.grid and contour to make a contour plot of g(x, n) over specified ranges for x = seq(0, x.max, length=npts) and n = seq(0, n.max, npts) for a specified number of points npts. Then I'd play with x.max, n.max, and npts until I got what I wanted. With the right choices for x.max, n.max, and npts, the solution will be obvious from the plot. In some cases, nothing more will be required. If I wanted more than that, I would need to exploit further some specifics of the problem. For that, permit me to restate some of what I think I understood of your specific problem: (1) For fixed n, g(x, n) is monotonically decreasing in x0. (2) For fixed x, g(x, n) has only two local maxima, one at n=0 (or n=eps0, esp arbitrarily small) and the other at n2(x), say, with a local minimum in between at n1(x), say. With this, I would write functions to find n1(x) and n2(x) given x. I might not even need n1(x) if I could figure out how to obtain n2(x) without it. Then I'd make a plot with two lines (using plot and lines) of g(x, 0) and g(x, n2(x)) vs. x. By the time I'd done all that, if I still needed more, I'd probably have ideas about what else to do. hope this helps. spencer graves Gregory Gentlemen wrote: Im trying to ascertain whether or not the facilities of R are sufficient for solving an optimization problem I've come accross. Because of my limited experience with R, I would greatly appreciate some feedback from more frequent users. The problem can be delineated as such: A utility function, we shall call g is a function of x, n ... g(x,n). g has the properties: n 0, x lies on the real line. g may take values along the real line. g is such that g(x,n)=g(-x,n). g is a decreasing function of x for any n; for fixed x, g(x,n) is smooth and intially decreases upon reaching an inflection point, thereafter increasing until it reaches a maxima and then declinces (neither concave nor convex). My optimization problem is to find the largest positive x such that g(x,n) is less than zero for all n. In fact, because of the symmetry of g around x, we need only consider x 0. Such an x does exists in this problem, and of course g obtains a maximum value of 0 at some n for this value of x. my issue is writing some code to systematically obtain this value. Is R capable of handling such a problem? (i.e. through some sort of optimization fucntion, or some sort of grid search with the relevant constraints) Any suggestions would be appreciated. Gregory Gentlemen [EMAIL PROTECTED] The following is a sketch of an optimization problem I need to solve. __ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] optimization problem in R ... can this be done?
The precision is not a problem, only the display, as Uwe indicated. Consider the following: (seq(25.5,25.6,length=20)-25.5)[c(1, 2, 19, 20)] [1] 0.00e+00 5.25e-07 9.50e-02 1.00e-01 ?options options(digits=20) seq(25.5,25.6,length=20)[c(1, 2, 19, 20)] [1] 25.5 25.50525 25.59475 25.6 spencer graves Gregory Gentlemen wrote: Okay let me attempt to be clear: if i construct the following sequence in R: seq(25.5,25.6,length=20) For instance, the last 10 elements of the sequence are all 26, and the preceding 20 are all 25.5. Presumably some rounding up is being done. How do I adjust the precision here such that each element is distinct? Thanks in advance guys, Gregory [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] Uwe Ligges [EMAIL PROTECTED] wrote: Gregory Gentlemen wrote: Spencer: Thank you for the helpful suggestions. I have another question following some code I wrote. The function below gives a crude approximation for the x of interest (that value of x such that g(x,n) is less than 0 for all n). # // btilda optimize g(n,x) for some fixed x, and then approximately finds that g(n,x) # such that abs(g(n*,x)=0 // btilda - function(range,len) { # range: over which to look for x bb - seq(range[1],range[2],length=len) OBJ - sapply(bb,function(x) {fixed - c(x,100,692,50,1600,v1,217227); return(optimize(g,c(1,1000),maximum=TRUE,tol=0.001,x=fixed)$objective)}) tt - data.frame(b=bb,obj=OBJ) tt$absobj - abs(tt$obj) d - tt[order(tt$absobj),][1:3,] return(as.vector(d)) } ! For instance a run of btilda(c(20.55806,20.55816),1) returns: b obj absobj 5834 20.55812 -0.0004942848 0.0004942848 5833 20.55812 0.0011715433 0.0011715433 5835 20.55812 -0.0021601140 0.0021601140 My question is how to improve the precision of b (which is x) here. It seems that seq(20.55806,20.55816,length=1 ) is only precise to 5 or so digits, and thus, is equivalent for numerous succesive Why do you think so? It is much more accurate! See ?.Machine Uwe Ligges values. How can I get around this? Spencer Graves wrote: Part of the R culture is a statement by Simon Blomberg immortalized in library(fortunes) as, This is R. There is no if. Only how. I can't see now how I would automate a complete solution to your problem in general. However, given a specific ! g(x, n), I would start by writing a function to use expand.grid and contour to make a contour plot of g(x, n) over specified ranges for x = seq(0, x.max, length=npts) and n = seq(0, n.max, npts) for a specified number of points npts. Then I'd play with x.max, n.max, and npts until I got what I wanted. With the right choices for x.max, n.max, and npts, the solution will be obvious from the plot. In some cases, nothing more will be required. If I wanted more than that, I would need to exploit further some specifics of the problem. For that, permit me to restate some of what I think I understood of your specific problem: (1) For fixed n, g(x, n) is monotonically decreasing in x0. (2) For fixed x, g(x, n) has only two local maxima, one at n=0 (or n=eps0, esp arbitrarily small) and the other at n2(x), say, with a local minimum in betwee! n at n1(x), say. With this, I would write functions to find n1(x) and n2(x) given x. I might not even need n1(x) if I could figure out how to obtain n2(x) without it. Then I'd make a plot with two lines (using plot and lines) of g(x, 0) and g(x, n2(x)) vs. x. By the time I'd done all that, if I still needed more, I'd probably have ideas about what else to do. hope this helps. spencer graves Gregory Gentlemen wrote: Im trying to ascertain whether or not the facilities of R are sufficient for solving an optimization problem I've come accross. Because of my limited experience with R, I would greatly appreciate some feedback from more frequent users. The problem can be delineated as such: A utility function, we shall call g is a function of x, n ... g(x,n)! . g has the properties: n 0, x lies on the real line. g may take values along the real line. g is such that g(x,n)=g(-x,n). g is a decreasing function of x for any n; for fixed x, g(x,n) is smooth and intially decreases upon reaching an inflection point, thereafter increasing until it reaches a maxima and then declinces (neither concave nor convex). My
Re: [R] optimization problem in R ... can this be done?
Part of the R culture is a statement by Simon Blomberg immortalized in library(fortunes) as, This is R. There is no if. Only how. I can't see now how I would automate a complete solution to your problem in general. However, given a specific g(x, n), I would start by writing a function to use expand.grid and contour to make a contour plot of g(x, n) over specified ranges for x = seq(0, x.max, length=npts) and n = seq(0, n.max, npts) for a specified number of points npts. Then I'd play with x.max, n.max, and npts until I got what I wanted. With the right choices for x.max, n.max, and npts, the solution will be obvious from the plot. In some cases, nothing more will be required. If I wanted more than that, I would need to exploit further some specifics of the problem. For that, permit me to restate some of what I think I understood of your specific problem: (1) For fixed n, g(x, n) is monotonically decreasing in x0. (2) For fixed x, g(x, n) has only two local maxima, one at n=0 (or n=eps0, esp arbitrarily small) and the other at n2(x), say, with a local minimum in between at n1(x), say. With this, I would write functions to find n1(x) and n2(x) given x. I might not even need n1(x) if I could figure out how to obtain n2(x) without it. Then I'd make a plot with two lines (using plot and lines) of g(x, 0) and g(x, n2(x)) vs. x. By the time I'd done all that, if I still needed more, I'd probably have ideas about what else to do. hope this helps. spencer graves Gregory Gentlemen wrote: Im trying to ascertain whether or not the facilities of R are sufficient for solving an optimization problem I've come accross. Because of my limited experience with R, I would greatly appreciate some feedback from more frequent users. The problem can be delineated as such: A utility function, we shall call g is a function of x, n ... g(x,n). g has the properties: n 0, x lies on the real line. g may take values along the real line. g is such that g(x,n)=g(-x,n). g is a decreasing function of x for any n; for fixed x, g(x,n) is smooth and intially decreases upon reaching an inflection point, thereafter increasing until it reaches a maxima and then declinces (neither concave nor convex). My optimization problem is to find the largest positive x such that g(x,n) is less than zero for all n. In fact, because of the symmetry of g around x, we need only consider x 0. Such an x does exists in this problem, and of course g obtains a maximum value of 0 at some n for this value of x. my issue is writing some code to systematically obtain this value. Is R capable of handling such a problem? (i.e. through some sort of optimization fucntion, or some sort of grid search with the relevant constraints) Any suggestions would be appreciated. Gregory Gentlemen [EMAIL PROTECTED] The following is a sketch of an optimization problem I need to solve. __ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com http://www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] optimization problem in R ... can this be done?
Spencer: Thank you for the helpful suggestions. I have another question following some code I wrote. The function below gives a crude approximation for the x of interest (that value of x such that g(x,n) is less than 0 for all n). # // btilda optimize g(n,x) for some fixed x, and then approximately finds that g(n,x) #such that abs(g(n*,x)=0 // btilda - function(range,len) { # range: over which to look for x bb - seq(range[1],range[2],length=len) OBJ - sapply(bb,function(x) {fixed - c(x,100,692,50,1600,v1,217227); return(optimize(g,c(1,1000),maximum=TRUE,tol=0.001,x=fixed)$objective)}) tt - data.frame(b=bb,obj=OBJ) tt$absobj - abs(tt$obj) d - tt[order(tt$absobj),][1:3,] return(as.vector(d)) } For instance a run of btilda(c(20.55806,20.55816),1) returns: b objabsobj 5834 20.55812-0.00049428480.0004942848 5833 20.55812 0.00117154330.0011715433 5835 20.55812-0.00216011400.0021601140 My question is how to improve the precision of b (which is x) here. It seems that seq(20.55806,20.55816,length=1 ) is only precise to 5 or so digits, and thus, is equivalent for numerous succesive values. How can I get around this? Spencer Graves [EMAIL PROTECTED] wrote: Part of the R culture is a statement by Simon Blomberg immortalized in library(fortunes) as, This is R. There is no if. Only how. I can't see now how I would automate a complete solution to your problem in general. However, given a specific g(x, n), I would start by writing a function to use expand.grid and contour to make a contour plot of g(x, n) over specified ranges for x = seq(0, x.max, length=npts) and n = seq(0, n.max, npts) for a specified number of points npts. Then I'd play with x.max, n.max, and npts until I got what I wanted. With the right choices for x.max, n.max, and npts, the solution will be obvious from the plot. In some cases, nothing more will be required. If I wanted more than that, I would need to exploit further some specifics of the problem. For that, permit me to restate some of what I think I understood of your specific problem: (1) For fixed n, g(x, n) is monotonically decreasing in x0. (2) For fixed x, g(x, n) has only two local maxima, one at n=0 (or n=eps0, esp arbitrarily small) and the other at n2(x), say, with a local minimum in between at n1(x), say. With this, I would write functions to find n1(x) and n2(x) given x. I might not even need n1(x) if I could figure out how to obtain n2(x) without it. Then I'd make a plot with two lines (using plot and lines) of g(x, 0) and g(x, n2(x)) vs. x. By the time I'd done all that, if I still needed more, I'd probably have ideas about what else to do. hope this helps. spencer graves Gregory Gentlemen wrote: Im trying to ascertain whether or not the facilities of R are sufficient for solving an optimization problem I've come accross. Because of my limited experience with R, I would greatly appreciate some feedback from more frequent users. The problem can be delineated as such: A utility function, we shall call g is a function of x, n ... g(x,n). g has the properties: n 0, x lies on the real line. g may take values along the real line. g is such that g(x,n)=g(-x,n). g is a decreasing function of x for any n; for fixed x, g(x,n) is smooth and intially decreases upon reaching an inflection point, thereafter increasing until it reaches a maxima and then declinces (neither concave nor convex). My optimization problem is to find the largest positive x such that g(x,n) is less than zero for all n. In fact, because of the symmetry of g around x, we need only consider x 0. Such an x does exists in this problem, and of course g obtains a maximum value of 0 at some n for this value of x. my issue is writing some code to systematically obtain this value. Is R capable of handling such a problem? (i.e. through some sort of optimization fucntion, or some sort of grid search with the relevant constraints) Any suggestions would be appreciated. Gregory Gentlemen [EMAIL PROTECTED] The following is a sketch of an optimization problem I need to solve. __ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html -- Spencer Graves, PhD Senior Development Engineer PDF Solutions, Inc. 333 West San Carlos Street Suite 700 San Jose, CA 95110, USA [EMAIL PROTECTED] www.pdf.com Tel: 408-938-4420 Fax: 408-280-7915 __ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting
[R] optimization problem in R ... can this be done?
Im trying to ascertain whether or not the facilities of R are sufficient for solving an optimization problem I've come accross. Because of my limited experience with R, I would greatly appreciate some feedback from more frequent users. The problem can be delineated as such: A utility function, we shall call g is a function of x, n ... g(x,n). g has the properties: n 0, x lies on the real line. g may take values along the real line. g is such that g(x,n)=g(-x,n). g is a decreasing function of x for any n; for fixed x, g(x,n) is smooth and intially decreases upon reaching an inflection point, thereafter increasing until it reaches a maxima and then declinces (neither concave nor convex). My optimization problem is to find the largest positive x such that g(x,n) is less than zero for all n. In fact, because of the symmetry of g around x, we need only consider x 0. Such an x does exists in this problem, and of course g obtains a maximum value of 0 at some n for this value of x. my issue is writing some code to systematically obtain this value. Is R capable of handling such a problem? (i.e. through some sort of optimization fucntion, or some sort of grid search with the relevant constraints) Any suggestions would be appreciated. Gregory Gentlemen [EMAIL PROTECTED] The following is a sketch of an optimization problem I need to solve. __ [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Optimization of constrained linear least-squares problem
Thanx Dimitris, Patrick and Berwin! For other people interested in this problem, here are two valid solutions that work. 1) Re-parameterize e.g., EM - c(100,0,0,0,100,0,0,0,100) W - array(EM, c(3,3)) d - c(10, 20, 70) fn - function(x){ x - exp(x) / sum(exp(x)) r - W%*%x - d crossprod(r, r)[1,1] } opt - optim(rnorm(3), fn) res - exp(opt$par) / sum(exp(opt$par)) res The first line of the `fn()' function is just a re-pameterization of your problem, i.e., if `y' is a vector of real numbers, then it is straightforward to see that `x = exp(y) / sum(exp(y))' will be real numbers in (0, 1) for which `sum(y)=1'. So instead of finding xs that minimize your function under the constraint (which is more difficult) you just find the ys using the above transformation. (I owe you a drink Dimitris !!!) 2) Or minimize it as a quadratic function under a linear constraint: EM - c(100,0,0,0,100,0,0,0,100) W - array(EM, c(3,3)) d - c(10, 20, 70) library(quadprog) Dmat - crossprod(W,W) dvec - crossprod(d,W) A -matrix(c(1,1,1),ncol=1) bvec - 1 solve.QP(Dmat, dvec, A, bvec, meq=1) This is based on the objective function (i.e. the thing you want to minimise) : min x'C'Cx - 2 d'Cx + d'd where sum(x) = 1 (Thanx Berwin!!) Kind regards, Stef __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Optimization of constrained linear least-squares problem
Dear R-ians, I want to perform an linear unmixing of image pixels in fractions of pure endmembers. Therefore I need to perform a constrained linear least-squares problem that looks like : min || Cx - d || ² where sum(x) = 1. I have a 3x3 matrix C, containing the values for endmembers and I have a 3x1 column vector d (for every pixel in the image). In theory my x values should all be in the (0,1) interval but I don't want to force them so I can check the validity of my solution. I just want to calculate the x values. Can anyone help me with this problem? I've been checking the optim, optimize, constrOptim and nlm help files, burt I don't understand it very well. Wich function should I use for my problem? I did a first test using optim: # Make my C matrix EM- c(4.5000,6.,10.5000,5.,27.,20.7500,16.7500,23.,38.7500) C - array(EM, c(3,3)) # Take an arbitrary d d-c(10, 20, 20) # Define the function fr - function(x) { C[1,]*x=d C[2,]*x=d C[3,]*x=d sum(x)=1} # Perform the optimization optim(c(0.25,0.25,0.25),fr) But it did not work. I got the eror couldn't find function. Can anyone tell me what functyion I should use for my problem and how should I program it? Thanx in advance, Stef __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Optimization of constrained linear least-squares problem
you could re-parameterize, e.g., EM - c(4.5000,6.,10.5000,5.,27.,20.7500,16.7500,23.,38.7500) W - array(EM, c(3,3)) d - c(10, 20, 20) ##33 fn - function(x){ x - exp(x) / sum(exp(x)) r - W%*%x - d crossprod(r, r)[1,1] } opt - optim(rnorm(3), fn) res - exp(opt$par) / sum(exp(opt$par)) res I hope it helps. Best, Dimitris Dimitris Rizopoulos Ph.D. Student Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/16/336899 Fax: +32/16/337015 Web: http://www.med.kuleuven.ac.be/biostat/ http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm - Original Message - From: Stefaan Lhermitte [EMAIL PROTECTED] To: r-help@stat.math.ethz.ch Sent: Thursday, March 17, 2005 3:13 PM Subject: [R] Optimization of constrained linear least-squares problem Dear R-ians, I want to perform an linear unmixing of image pixels in fractions of pure endmembers. Therefore I need to perform a constrained linear least-squares problem that looks like : min || Cx - d || ² where sum(x) = 1. I have a 3x3 matrix C, containing the values for endmembers and I have a 3x1 column vector d (for every pixel in the image). In theory my x values should all be in the (0,1) interval but I don't want to force them so I can check the validity of my solution. I just want to calculate the x values. Can anyone help me with this problem? I've been checking the optim, optimize, constrOptim and nlm help files, burt I don't understand it very well. Wich function should I use for my problem? I did a first test using optim: # Make my C matrix EM- c(4.5000,6.,10.5000,5.,27.,20.7500,16.7500,23.,38.7500) C - array(EM, c(3,3)) # Take an arbitrary d d-c(10, 20, 20) # Define the function fr - function(x) { C[1,]*x=d C[2,]*x=d C[3,]*x=d sum(x)=1} # Perform the optimization optim(c(0.25,0.25,0.25),fr) But it did not work. I got the eror couldn't find function. Can anyone tell me what functyion I should use for my problem and how should I program it? Thanx in advance, Stef __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] R Optimization 101 for R Newcomers: An extended example
To new R-Programmers: This long-winded note gives an example of optimizing in R that should only be of interest to newcomers to the language. Others should ignore. My hope is that it might help illuminate some basic notions of code improvement, looping, and vectorization in R. However, I welcome comments from anyone that enlarge, correct, or improve it. I remain a novice R programmer. For the curious, the hardware/opsys details are a 1.5 ghz Windows 2000 machine, R1.9.0. The task is this: I have a data set, mydat, consisting of values =0, and a fixed constant, K, that is much larger than any of the values. I want to simulate the distribution of the number of bootstrap draws (random sampling with replacement) from mydat that is needed to make the sum of the values drawn just exceed K. That is, repeatedly bootstrap sample from mydat until the sum of the results exceeds K. Suppose this requires 7 draws. Then 7 is a random sample of size 1 from the distribution I'm interested in. I want to repeat this procedure a large number, nsim, of times to estimate the distribution. Before reading on, I invite R novices to consider how -- or better yet write code -- to do this. Approach 1 (Naive/ straightforward): Perform the simulation exactly as described by using a nested loop. The outer loop goes from 1 to nsim to record the number of draws needed. The inner loop repeatedly draws a sample of size 1 from mydat and accumulates and checks the sum until it exceeds K, returning the number of times it took to do this to the outer loop. This is easy to program up using a while loop for the inner loop and takes about a short look's time out my window (a few seconds) to run for modest nsim (5-10K). But surely one can do better... Approach 2: Was essentially the same, except all the sampling is done at once. That is, each execution of the inner loop in (1) required R to call the sample() function to get a bootstrap sample of size 1. This is very inefficient, because there is a high overhead of making that function call at each iteration. It can be avoided by calling sample() only once -- or at most a small number of times -- to get a large sample and then using indexing to work one's way through the large sample. This turns out to be a little more complicated than approach 1 because you first have to figure how large a large sample should be and then have to be a little careful about setting up your indexing. But it's not very difficult (especially for any C programmer who can manipulate pointers), and random sampling is so fast that it's no problem to generate a sample WAYYY bigger than you need (5 or 10 times as large, even) just to be safe. Or generate a not so big version and just check at each outer loop iteration whether you need to generate some more. Doing it this way -- now using indexing, not function calls in the inner loop -- made a considerable improvement (5 - 10 fold). It now took about one quick glance out the window time to generate the distribution (less than a second). This was more than adequate for my needs under any conceivable situation. But still, self-respecting R programmers are supposed to eschew looping, and here was this ugly(?) loop within a loop! So let us persevere. Approach 3: So the question is -- how to remove that inner loop? Again, I invite novice R programmers to think about that before continuing on... The key here is that the inner loop is doing a very simple calculation: just accumulating a sum. Forsooth! -- the answer fairly leaps out: R already has a cumsum() function that does this (looping in the underlying C code), so one only has to figure out how to make use of it. The essential idea to do this is to get a small sample from the large mydat sample that is guaranteed to be bigger than one needs to total up to more than K. Again, one can be profligate: if I need about 10 mydat values on average to sum up to K, if I sample 30 or 40 or some such thing, I'm sure to have enough (and can always add a check or use try() to make sure if I am faint of heart). If smallsamp is this sample of the next 30 or 40 values from my large bootstrap sample, then the following code vectorizes the inner loops: count - sum(cumsum(smallsamp)K)+1 Note a trick here: The K produces a vector of logical TRUE's for all cumsum values K. This is silently treated as numeric 1's when added by sum(). A moment's thought will make clear why 1 must be added at the end.. Sure enough, this vectorization reduced the time about another twofold -- to an eyeblink -- for my modest nsim sizes. Approach 4: The last stage, of course, is to get rid of the outer loop, which fills the vector for the distribution one element at a time, again another R no-no. One way to do this -- the ONLY way I could think of (so a challenge to smarter R programmers) -- is to make use of apply(), which basically always can remove loops. Unfortunately, this is an illusion, because what the apply() functions actually do is hide the
[R] Optimization failed in fitting mixture 3-parameter Weibulldistri bution using fitdistr()
Dear All; I tried to use fitdistr() in the MASS library to fit a mixture distribution of the 3-parameter Weibull, but the optimization failed. Looking at the source code, it seems to indicate the error occurs at if (res$convergence 0) stop(optimization failed). The procedures I tested are as following: w3den - function(x, a,b,c) {c/b*((x -a)/b)^(c-1)*exp(-((x-a)/b)^c)}# define 3-parameter Weibull density w3den - function(x, a,b,c) {c/b*((x -a)/b)^(c-1)*exp(-((x-a)/b)^c)} set.seed(123) x3 - rweibull(100, shape = 4, scale = 100) # Distribution 1 fitdistr(x3, w3den, start= list(a = 8.8, b = 90.77, c = 3.678)) # Fitting 3-parameter abc 8.9487415 90.66677123.6722124 (15.1462445) (16.0657103) ( 0.7582376) Warning messages: 1: NaNs produced in: log(x) 2: NaNs produced in: log(x) x4 - rweibull(50, shape = 3, scale = 90) # Distribution 2 fitdistr(x4, w3den, start=list(a = 5.0, b = 90, c = 3)) abc -67.572244 159.020786 5.835588 (105.975870) (107.842846) ( 4.288019) Warning messages: 1: NaNs produced in: log(x) 2: NaNs produced in: log(x) mweib - function(x, p, a,b,c,a1, b1, c1) {p*(c/b*((x - a)/b)^(c-1)*exp(-((x-a)/b)^c))+ + +(1-p)*(c1/b1*((x -a1)/b1)^(c1-1)*exp(-((x-a1)/b1)^c1))} # define mixture distribution x5 - c(x3, x4) fitdistr(x5, mweib, start=list(p = 0.7, a = 8.9, b=90.77, c = 3.68, a1 = -67.57, b1 = 159.02, c1 = 5.83)) Error in fitdistr(x5, mweib, start = list(p = 0.7, a = 8.9, b = 90.77, : optimization failed In addition: There were 14 warnings (use warnings() to see them). I tested the same procedures with a mixture 2-parameter Weibull distribution without problems. With a limited experience with the fitdistr(), it seems to me, similar to all nonlinear or optimization procedures, the starting values are critical for convergence. Any suggestions for solving the optimization problem? TIA Richard Yang Northern Forestry Centre /Centre de foresterie du Nord Canadian Forest Service/Service canadien des forêts Natural Resources Canada /Ressources naturelles Canada 5320-122 Street/5320, rue 122 Edmonton (Alberta) Canada T6H 3S5 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] Optimization failed in fitting mixture 3-parameter Weibulldistri bution using fitdistr()
I don't think that is the right density: haven't you forgotten I(x a)? So you need a constraint on a in the optimization, or at least to return density 0 if a = min(x_i) (but I suspect the mle may well occur at the boundary). Without that constraint you don't have a valid optimization problem. On Mon, 28 Jul 2003, Yang, Richard wrote: Dear All; I tried to use fitdistr() in the MASS library to fit a mixture distribution of the 3-parameter Weibull, but the optimization failed. Looking at the source code, it seems to indicate the error occurs at if (res$convergence 0) stop(optimization failed). The procedures I tested are as following: w3den - function(x, a,b,c) {c/b*((x -a)/b)^(c-1)*exp(-((x-a)/b)^c)}# define 3-parameter Weibull density w3den - function(x, a,b,c) {c/b*((x -a)/b)^(c-1)*exp(-((x-a)/b)^c)} set.seed(123) x3 - rweibull(100, shape = 4, scale = 100) # Distribution 1 fitdistr(x3, w3den, start= list(a = 8.8, b = 90.77, c = 3.678)) # Fitting 3-parameter abc 8.9487415 90.66677123.6722124 (15.1462445) (16.0657103) ( 0.7582376) Warning messages: 1: NaNs produced in: log(x) 2: NaNs produced in: log(x) x4 - rweibull(50, shape = 3, scale = 90) # Distribution 2 fitdistr(x4, w3den, start=list(a = 5.0, b = 90, c = 3)) abc -67.572244 159.020786 5.835588 (105.975870) (107.842846) ( 4.288019) Warning messages: 1: NaNs produced in: log(x) 2: NaNs produced in: log(x) mweib - function(x, p, a,b,c,a1, b1, c1) {p*(c/b*((x - a)/b)^(c-1)*exp(-((x-a)/b)^c))+ + +(1-p)*(c1/b1*((x -a1)/b1)^(c1-1)*exp(-((x-a1)/b1)^c1))} # define mixture distribution x5 - c(x3, x4) fitdistr(x5, mweib, start=list(p = 0.7, a = 8.9, b=90.77, c = 3.68, a1 = -67.57, b1 = 159.02, c1 = 5.83)) Error in fitdistr(x5, mweib, start = list(p = 0.7, a = 8.9, b = 90.77, : optimization failed In addition: There were 14 warnings (use warnings() to see them). I tested the same procedures with a mixture 2-parameter Weibull distribution without problems. With a limited experience with the fitdistr(), it seems to me, similar to all nonlinear or optimization procedures, the starting values are critical for convergence. Any suggestions for solving the optimization problem? TIA Richard Yang Northern Forestry Centre / Centre de foresterie du Nord Canadian Forest Service /Service canadien des forêts Natural Resources Canada / Ressources naturelles Canada 5320-122 Street /5320, rue 122 Edmonton (Alberta) Canada T6H 3S5 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help