## [ab826e]: mle_example.m  Maximize  Restore  History

### 88 lines (72 with data), 3.4 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 ## Copyright (C) 2003,2004 Michael Creel ## ## 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 ## Example to show how to use MLE functions # 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"; modelargs = {0}; # if this is zero the function gives analytic score, otherwise not # 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 = {50,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); modelargs = {1}; # set the switch for no score [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control); # Example doing estimation in parallel on a cluster (requires MPITB) if exist("MPI_Init") printf("to do estimation in parallel, you need to have MPITB installed\n... and your computer must be lambooted. If this is not the case press... CRTL-C to abort. Pausing 10 seconds\n"); pause(10); theta = zeros(3,1); nslaves = 1; title = "MLE estimation done in parallel"; [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control, nslaves); else printf("sorry, MPITB is not installed, can't do estimation in parallel\n"); endif