[R] Optimization under an absolute value constraint

2007-09-07 Thread Phil Xiang
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

2007-09-07 Thread roger koenker
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

2007-09-07 Thread Paul Smith
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

2007-08-22 Thread Gabor Grothendieck
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

2007-08-21 Thread Alan Harrison
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

2007-08-21 Thread Gabor Grothendieck
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

2007-08-21 Thread Ben Bolker

  (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

2007-07-18 Thread Tobias Schlottmann
  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

2007-07-17 Thread Paul Smith
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

2007-07-17 Thread massimiliano\.talarico
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

2007-07-17 Thread Berwin A Turlach
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

2007-07-17 Thread Stephen Tucker
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

2007-07-17 Thread Stephen Tucker
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

2007-07-17 Thread Talarico Massimiliano \(UniCredit Xelion Banca\)
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

2007-07-17 Thread massimiliano\.talarico
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

2007-07-17 Thread Stephen Tucker
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

2007-07-17 Thread Ravi Varadhan
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

2007-07-17 Thread Moshe Olshansky
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

2007-07-17 Thread Berwin A Turlach
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

2007-07-17 Thread Moshe Olshansky
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

2007-07-17 Thread ecatchpole
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

2007-07-16 Thread massimiliano\.talarico
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

2007-06-19 Thread livia

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

2007-06-18 Thread livia

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

2007-06-18 Thread Robert A LaBudde
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

2007-06-18 Thread Uwe Ligges


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

2007-06-18 Thread Ted Harding
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

2007-06-18 Thread Prof Brian Ripley
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

2007-06-18 Thread Petr Klasterecky
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

2007-04-03 Thread Jin Huang
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

2006-11-01 Thread stat stat
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

2006-11-01 Thread Patrick Burns
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

2006-09-05 Thread Simone Vincenzi
 

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

2006-09-05 Thread Spencer Graves
 = 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

2006-09-04 Thread Spencer Graves
  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

2006-08-30 Thread Simone Vincenzi
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

2006-03-06 Thread Jasjeet S. Sekhon

 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

2006-03-04 Thread Spencer Graves
  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

2006-03-02 Thread nojhan
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

2006-03-01 Thread Mark
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

2006-03-01 Thread Berton Gunter
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

2006-03-01 Thread Gabor Grothendieck
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

2005-11-28 Thread Florent Bresson
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

2005-11-28 Thread Peter Dalgaard
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

2005-11-28 Thread Prof Brian Ripley
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

2005-11-28 Thread Gabor Grothendieck
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?

2005-06-26 Thread Uwe Ligges
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?

2005-06-26 Thread Spencer Graves
  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?

2005-06-25 Thread Spencer Graves
  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?

2005-06-25 Thread Gregory Gentlemen
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?

2005-06-24 Thread Gregory Gentlemen
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

2005-03-18 Thread Stefaan Lhermitte
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

2005-03-17 Thread Stefaan Lhermitte
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

2005-03-17 Thread Dimitris Rizopoulos
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

2004-05-19 Thread Berton Gunter
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()

2003-07-28 Thread Yang, Richard
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()

2003-07-28 Thread Prof Brian Ripley
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