Section 7.1.1: Counting problems with Poisson distribution

% Boyd & Vandenberghe "Convex Optimization"
% Joëlle Skaf - 04/24/08
%
% The random variable y is nonnegative and integer valued with a Poisson
% distribution with mean mu > 0. In a simple statistical model, the mean mu
% is modeled as an affine function of a vector u: mu = a'*u + b.
% We are given a number of observations which consist of pairs (u_i,y_i),
% i = 1,..., m, where y_i is the observed value of y for which the value of
% the explanatory variable is u_i. We find a maximum likelihood estimate of
% the model parameters a and b from these data by solving the problem
%           maximize    sum_{i=1}^m (y_i*log(a'*u_i + b) - (a'*u_i + b))
% where the variables are a and b.

% Input data
rand('state',0);
n = 10;
m = 100;
atrue = rand(n,1);
btrue = rand;

u = rand(n,m);
mu = atrue'*u + btrue;

% Generate random variables y from a Poisson distribution
% (The distribution is actually truncated at 10*max(mu) for simplicity)
L  = exp(-mu);
ns = ceil(max(10*mu));
y  = sum(cumprod(rand(ns,m))>=L(ones(ns,1),:));

% Maximum likelihood estimate of model parameters
cvx_begin
    variables a(n) b(1)
    maximize sum(y.*log(a'*u+b) - (a'*u+b))
cvx_end
 
Successive approximation method to be employed.
   For improved efficiency, sedumi is solving the dual problem.
   sedumi will be called several times to refine the solution.
   Original size: 276 variables, 103 equality constraints
   92 exponentials add 736 variables, 460 equality constraints
-----------------------------------------------------------------
          Errors   
Act Centering    Conic    Status
-----------------------------------
92  6.171e-01  3.305e-02  Solved
92  8.713e-02  6.259e-04  Solved
92  1.348e-02  1.565e-05  Solved
92  1.754e-03  5.349e-07  Solved
92  2.200e-04  2.265e-07  Solved
92  2.791e-05  2.160e-07S Solved
92  3.407e-06  2.158e-07S Solved
92  3.030e-04S 2.188e-08  Solved
92  8.477e-05  3.812e-09  Solved
92  1.063e-05  3.362e-09  Solved
92  1.324e-06  3.389e-09S Solved
92  3.283e-06S 1.596e-09  Inaccurate/Solved
92  2.130e-05S 0.000e+00  Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +102.57