[eb38ba]: mle_example.m  Maximize  Restore  History

Download this file

134 lines (113 with data), 5.1 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
# Copyright (C) 2003,2004 Michael Creel michael.creel@uab.es
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@uab.es
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@uab.es
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1;
# Example to show how to use MLE functions
# Example likelihood function, no score
function [log_density, score] = poisson_no_score(theta, data, otherargs)
y = data(:,1);
x = data(:,2:columns(data));
lambda = exp(x*theta);
log_density = -lambda + y .* (x*theta) - lgamma(y+1);
score = "na";
endfunction
# Example likelihood function, with score
function [log_density, score] = poisson_with_score(theta, data, otherargs)
y = data(:,1);
x = data(:,2:columns(data));
lambda = exp(x*theta);
log_density = -lambda + y .* (x*theta) - lgamma(y+1);
score = dmult(y - lambda,x);
endfunction
# Generate data
n = 1000; # how many observations?
# the explanatory variables: note that they have unequal scales
x = [ones(n,1) rand(n,1) randn(n,1)];
theta = 1:3; # true coefficients are 1,2,3
theta = theta';
lambda = exp(x*theta);
y = randp(lambda); # generate the dependent variable
####################################
# define arguments for mle_results #
####################################
# starting values
theta = zeros(3,1);
# data
data = [y, x];
# name of model to estimate
model = "poisson_with_score";
# placeholder, poisson model has no additional args
modelargs = {};
# parameter names
names = str2mat("beta1", "beta2", "beta3");
title = "Poisson MLE trial"; # title for the run
# controls for bfgsmin: 30 iterations is not always enough for convergence
control = {30,0,1,1};
# This displays the results
printf("\n\nanalytic score, unscaled data\n");
[theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, 0, control);
# This just calculates the results, no screen display
printf("\n\nanalytic score, unscaled data, no screen display\n");
theta = zeros(3,1);
[theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control);
printf("obj_value = %f, to confirm it worked\n", obj_value);
# This show how to scale data during estimation, but get results
# for data in original units (recommended to avoid conditioning problems)
# This usually converges faster, depending upon the data
printf("\n\nanalytic score, scaled data\n");
[scaled_x, unscale] = scale_data(x);
data = [y, scaled_x];
theta = zeros(3,1);
[theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control);
# Example using numeric score
printf("\n\nnumeric score, scaled data\n");
theta = zeros(3,1);
model = "poisson_no_score";
[theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control);

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks