1
$\begingroup$

The following is an assignment.

The matrix $\mathbf{x}$ containing the external factors has dimension $4\times1000$, and the vector $\mathbf{y}$ containing the categorical variable has dimension $1\times1000$.

I have to fit a GLM model of the type

$$ \pi(t) = \frac{1}{1 + exp{\left(\beta_0 + \sum_{l=1}^m \beta_l u_t^l\right)}}, $$

where $u_t^l$ is the $l$th external factor at time $t$.

The likelihood function I came up with is the following:

$$ L(\beta) = -\sum_{i=1}^N \left[ y_i \log \left( \frac{1}{1+\exp(-\beta^T\mathbf{x}_i)} \right) +\\ (1-y_i)\log \left( 1 - \frac{1}{1+\exp(-\beta^T\mathbf{x}_i)} \right) \right] .$$

In order to maximize it, I computed the gradient:

$$ \nabla_{\beta} L = \sum_{i=1}^N \left( y_i - \frac{1}{1+\exp(-\beta^T\mathbf{x}_i)} \right)\mathbf{x}_i ,$$

and the updating expression is:

$$ \beta = \beta + \nabla_{\beta} L .$$

I want to implement it in MATLAB and for doing so I would like to first transform this last expression in matrix form:

$$ \beta = \beta + \left( \mathbf{y} - \frac{1}{1+\exp(-\beta^T\mathbf{x})} \right)\mathbf{x} .$$

I've been trying for hours now but the result is not correct (too many $1$s or $0$s according to the initial values).

Here's the question: can somebody explain what the dimensions of the different elements for the matrix form should be?

$\endgroup$
1
$\begingroup$

It should be pretty straightforward to code:

function llik = fun(b, X, Y)  
num = X * b;  
prb = exp(num .* Y) ./ (1 + exp(num));  
llik = -sum(log(prb));  
end

Where: (Y) is a column vector (say 1000 x 1)
(X) is a matrix of predictors (say 1000 x 5)
Key is exp(num .* Y) that will be used to obtain both proba(Y==1) and proba(Y==0)

Given the relative simplicity of this model and the efficiency of Matlab optimisation routines (fmincon/fminunc), I don't think you need to code the gradient, but can easily be done if really needed.

This model has a closed-form solution and then the results should in principle not depend on the choice of starting values. Hope this helps.

=== SCRIPT TO ESTIMATE THE MODEL ===
Here I assume that you first have created a script including the log-likelihood function only - Let's say this script is called "LOGISTIC_LL".

% IMPORT DATA
A = importdata('path.csv')

% DATA SPECIFICATION
Y = A.data(:,1);
X = A.data(:,2:end);

% VARIABLE NAMES
Vnames = {'x1','x2',etc.};

% STARTING VALUES
b = zeros(length(Vnames),1);

% OPTIMISATION
options = optimoptions(@fminunc,'Display','iter','MaxIterations',1e3,'MaxFunctionEvaluations',1e5);
tic;
[paramhat,fval,~,~,grad,hessian] = fminunc(@(b)LOGISTIC_LL(b, X, Y), b, options);

Finally, you could create another function to reshape results and compute other statistics (SEs, Pval, etc).

$\endgroup$
  • $\begingroup$ Forgot to say that (b) corresponds to the vector of model parameters $\endgroup$ – Umka May 27 '17 at 17:04
  • $\begingroup$ I tried to use fminsearch for the likelihood function but I had the same problems. I actually do not understand what you are doing here. Is this the function to be passed to an optimization routines? $\endgroup$ – wrong_path May 27 '17 at 17:10
  • $\begingroup$ Yes. I usually create a first script in which I write my objective function (say log-likelihood for a logistic regression) and then I create another file in which I estimate the model - For instance in this last file I would import the data, declare the different objects (Y, X, etc.) and run the other script. $\endgroup$ – Umka May 27 '17 at 17:12
  • $\begingroup$ I've edited my answer accordingly (see "SCRIPT TO ESTIMATE THE MODEL") $\endgroup$ – Umka May 27 '17 at 17:18
  • 1
    $\begingroup$ This is because what you model in a logistic regression is proba(Y==1) - This is why you obtain values between 0 and 1. If you really want something in (0/1) you could transform your predicted proba into predicted outcomes (e.g., if P(Y==1) > 0.5 then predicted outcome = 1 and otherwise) - Finally you can compare the predicted outcomes with the actual (observed) ones and compute a measure of agreement (% of correctly predicted events) - This will tell you something about model performance $\endgroup$ – Umka May 27 '17 at 17:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.