R help - Sep 2005 - regression with restrictions - optimization problem

Dear WizaRds!

I am sorry to ask for some help, but I have come to a complete stop in 
my efforts. I hope, though, that some of you might find the problem 
quite interesting to look at.

I have been trying to estimate parameters for lotteries, the so called 
utility of chance, i.e. the "felt" probability compared to a rational 
given probability. A real brief example: Given is a lottery payoff 
matrix of the type

x1    x2 ...     x10     median
1000    5000 ... 5000    3750
0    1000 ... 5000    2250 etc.

The actual data frame consists of 11 columns and 28 rows.

Each entry x1 ... x10 gives the amount of money resp. the utility of 
that amount you receive playing the lottery. The probability for each 
column is 10%. The median represents the empirical answers of players 
where the person is indifferent if they prefer to receive the lottery or 
the sum of money as a sure payoff.

I try to determine the probability people feel instead of the known 10% 
probability of each column payoff entry. But here's the catch:

People also give different utilities to each amount of money, which 
basically gives us some sort of function like this:
u(x1...x10) = u(x1)*pi(p1) + u(x2)*pi(p2) +...+u(x10)*pi(p10)=y
u() - unknown utility function
pi() - unknown probability function
y - empirical answer
p1..p10 - probabilities, here always 0.1

To keep it simple, I set u(0)=0 and u(5000)=5000 and vary u(1000) 
between a start and end point. On each cycle R computes the regression 
coefficients that serve as the pi(p) estimators for every 10% step.
Then I minimize the residual sum of squares which should give the best 
estimators for every 10% step.

How can I possibly calculate a "smooth" pi(p) curve, a curve that
should
look like an "S", plotted against the cumulative 10% probabilities? I 
only have my ten estimators. How can I possibly tell R the necessary 
restrictions of nonnegative estimators and their sum to equal one? Here 
is my quite naive approach:

a70 <- matrix(c(1000,5000,5000,5000,2150, 0,1000,5000,5000,1750, 
0,0,1000,5000,1150, 0,0,0,1000,200, 1000,1000,5000,5000,2050, 
0,1000,1000,5000,1972), ncol=5, byrow=T)
colnames(a70)=c(paste("x", 1:4, sep=""), "med")
a70 <- as.data.frame(a70)

start=800; end=2000
step=10; u1000=start-step

u1000 <- u1000+step # varying the 1000 entry
a70[a70==1000] <- u1000
reg70 <- lm(a70$med ~ -1+x1+x2+x3+x4, data=a70)
res <- sum( (reg70$residuals^2) )

for (i in 1:( (end-start)/step) ){
         a70[a70==u1000]    <- u1000+step
     u1000 <- u1000+step
     reg70 <- lm(a70$med ~ -1+x1+x2+x3+x4, data=a70)
     if (res >= sum( (reg70$residuals^2) )) {
         res <- sum( (reg70$residuals^2) )
         print(paste("cycle", i, "u1000=", u1000,
"RSS=", res))
         final70 <- a70
         finalreg <- reg70
         }
}
print(final70)
summary(finalreg)

     Maybe a better approach works with optim(stats) or dfp(Bhat), but I 
have no idea how to correctly approach such a restricted optimization 
problem.
     Thank you su much for your help and support.

Mark Hempelmann

I have not seen any replies, so I will offer a comment:

	  1.  You speak of x1, x2, ..., x10, but your example includes only 
x1+x2+x3+x4.  I'm confused.  If you could still use help with this, 
could you please simplify your example further so there was only x1+x2, 
say?  Can you solve the problem with only x1?  If no, state your problem 
in terms only of x1.  The answer to that may include a fairly obvious 
generalization to x10.  If not, that can become a follow-on question.

	  2.  What is your objective function?  What do you want to maximize or 
minimize?  Are your y's ("med"?), e.g., some model plus normal
error?
If yes, then some kind of (nonlinear) regression might be appropriate. 
If no, then it might be best to start by writing an objective function. 
  If you have only one parameter to estimate to minimize the objective 
function, then you can just compute the objective function over a range 
for the parameter, create a plot, and be done.  If you want some 
refinement of that, please look at "optimize".  If you only two 
parameters, you can create contour plots, and use "optim" to refine
the
result.  For more unknowns, "optim" can be fairly useful.

	  Good luck.
	  spencer graves

Mark Hempelmann wrote:
> Dear WizaRds!
> 
> I am sorry to ask for some help, but I have come to a complete stop in 
> my efforts. I hope, though, that some of you might find the problem 
> quite interesting to look at.
> 
> I have been trying to estimate parameters for lotteries, the so called 
> utility of chance, i.e. the "felt" probability compared to a
rational
> given probability. A real brief example: Given is a lottery payoff 
> matrix of the type
> 
> x1    x2 ...     x10     median
> 1000    5000 ... 5000    3750
> 0    1000 ... 5000    2250 etc.
> 
> The actual data frame consists of 11 columns and 28 rows.
> 
> Each entry x1 ... x10 gives the amount of money resp. the utility of 
> that amount you receive playing the lottery. The probability for each 
> column is 10%. The median represents the empirical answers of players 
> where the person is indifferent if they prefer to receive the lottery or 
> the sum of money as a sure payoff.
> 
> I try to determine the probability people feel instead of the known 10% 
> probability of each column payoff entry. But here's the catch:
> 
> People also give different utilities to each amount of money, which 
> basically gives us some sort of function like this:
> u(x1...x10) = u(x1)*pi(p1) + u(x2)*pi(p2) +...+u(x10)*pi(p10)=y
> u() - unknown utility function
> pi() - unknown probability function
> y - empirical answer
> p1..p10 - probabilities, here always 0.1
> 
> To keep it simple, I set u(0)=0 and u(5000)=5000 and vary u(1000) 
> between a start and end point. On each cycle R computes the regression 
> coefficients that serve as the pi(p) estimators for every 10% step.
> Then I minimize the residual sum of squares which should give the best 
> estimators for every 10% step.
> 
> How can I possibly calculate a "smooth" pi(p) curve, a curve that
should
> look like an "S", plotted against the cumulative 10%
probabilities? I
> only have my ten estimators. How can I possibly tell R the necessary 
> restrictions of nonnegative estimators and their sum to equal one? Here 
> is my quite naive approach:
> 
> a70 <- matrix(c(1000,5000,5000,5000,2150, 0,1000,5000,5000,1750, 
> 0,0,1000,5000,1150, 0,0,0,1000,200, 1000,1000,5000,5000,2050, 
> 0,1000,1000,5000,1972), ncol=5, byrow=T)
> colnames(a70)=c(paste("x", 1:4, sep=""),
"med")
> a70 <- as.data.frame(a70)
> 
> start=800; end=2000
> step=10; u1000=start-step
> 
> u1000 <- u1000+step # varying the 1000 entry
> a70[a70==1000] <- u1000
> reg70 <- lm(a70$med ~ -1+x1+x2+x3+x4, data=a70)
> res <- sum( (reg70$residuals^2) )
> 
> for (i in 1:( (end-start)/step) ){
>          a70[a70==u1000]    <- u1000+step
>      u1000 <- u1000+step
>      reg70 <- lm(a70$med ~ -1+x1+x2+x3+x4, data=a70)
>      if (res >= sum( (reg70$residuals^2) )) {
>          res <- sum( (reg70$residuals^2) )
>          print(paste("cycle", i, "u1000=", u1000,
"RSS=", res))
>          final70 <- a70
>          finalreg <- reg70
>          }
> }
> print(final70)
> summary(finalreg)
> 
>      Maybe a better approach works with optim(stats) or dfp(Bhat), but I 
> have no idea how to correctly approach such a restricted optimization 
> problem.
>      Thank you su much for your help and support.
> 
> Mark Hempelmann
> 
> ______________________________________________
> R-help at 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

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915