Re: problems with sqp constraints?

[fork]

fork <forkandwait <at> gmail.com> writes:

> 
> Hi,
> 
> I got sqp to run with an entropy maximization set up ...

Here is my self-teaching sqp() script.  It lacks inequality constraints, 
and only works with centered data, and only uses linear constraints, 
but it might still be useful as an example to someone.  Sorry to post such a big
thing, but I don't have a blog ;)

Basically, it finds a probability distribution ("x") for a random
variable ("val") that meets some moment criteria for that distribution. 

Corrections to the math are welcome, but be gentle... ;)

(The bigger point is that I am trying to figure out a way to integrate
a lot of different demographic information, at widely different scales and 
levels of accuracy, and come up with a most likely social configuration 
"underneath" these disparate data. I will model the data as 
constraints (persons_per_household = 2.8, occupancy = 0.92, etc), 
and let the objective function float.  I am not there yet, though...)


## the random variable with associated probabilities
val = -10:0.1:10;

## ##############################
## negative entropy function
## ##############################
function H = entr(p)
    p = p(:);    
    _h = (p .* (log(p)));
    _pm = (p == 0);
    h = ifelse(_pm, 0, _h);
    H = real(sum(h));
endfunction

## ##############################
## gradient of negative entropy function
## ##############################
function dh = entr_gr(p)
    dh = log(p) + 1;
endfunction

## ##############################
## hessian of negative entropy function 
## ##############################
function ddh = entr_h(p)
    ddh = diag(1 ./ p);
endfunction

## ##############################
## equality constr...
## ##############################
function f = mk_g (A, b)
    f = @(x) (A*x - b);
endfunction;
clear A b;

## probability definition
A(1,:)  = ones(1, length(val));
b(1,1)  = 1;

## set the mean to be something
A(2,:)  = val;
b(2,1)  = -0.10;

## set the variance
A(3,:)  = val .^ 2;
b(3,1)  = 25.5;

## set the skewness (ish -- 3rd moment)
A(4,:)  = val .^ 3;
b(4,1)  = 5;

## put it all in a matrix and curry it up into a function
g = mk_g(A, b);


## ##############################
## inequality constr -- nothing yet ...
## ##############################
function f = mk_h (A, b)
    f = @(x) (A*x - b);
endfunction;
clear A b;

## ##############################
## do everything
## ##############################
x0 = repmat(.01, length(val), 1);
lb = repmat(  0, length(val), 1);
ub = repmat(  1, length(val), 1);

## [X, OBJ, INFO, ITER, NF, LAMBDA] = sqp (X0, PHI, G,  H,  LB, UB, MAXITER, 
TOL)
##  g(x)==0, h(x)>=0

[x, obj, info, iter, nf, lambda] =    sqp (x0, address@hidden; @entr_gr; 
@entr_h}, g,
[], lb, ub);