From: Michael Creel <mcreel@us...>  20041020 09:12:50

Update of /cvsroot/octave/octaveforge/main/econometrics In directory sc8prcvs1.sourceforge.net:/tmp/cvsserv30326/main/econometrics Added Files: delta_method.m mle_estimate.m mle_example.m mle_obj.m mle_results.m mle_variance.m parameterize.m prettyprint_c.m prettyprint.m scale_data.m unscale_parameters.m Log Message: initial commit of econometrics functions  NEW FILE: scale_data.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Standardizes and normalizes data matrix, # primarily for use by BFGS function [zz, scalecoefs] = scale_data(z); n = rows(z); k = columns(z); # Scale data s = std(z)'; test = s != 0; s = s + (1  test); # don't scale if column is a constant (avoid div by zero) A = diag(1 ./ s); # Demean all variables except constant, if a constant is present test = std(z(:,1)) != 0; if test warning("ScaleData: no constant present  only scale  no demean"); bb = zeros(n,k); b = zeros(k,1); else b = mean(z)'; test = std(z)' != 0; # don't take out mean if the column is a constant, to preserve identification b = b .* test; b = A*b; bb = dmult(b, ones(k,n))'; endif zz = z*A + bb; scalecoefs = {A,b}; endfunction  NEW FILE: prettyprint.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # this prints matrices with row and column labels function prettyprint(mat, rlabels, clabels) # left pad the column labels a = columns(rlabels); for i = 1:a printf(" "); endfor printf(" "); # print the column labels clabels = [" ";clabels]; # pad to 8 characters wide clabels = strjust(clabels,"right"); k = columns(mat); for i = 1:k printf("%s ",clabels(i+1,:)); endfor # now print the row labels and rows printf("\n"); k = rows(mat); for i = 1:k if isstr(rlabels(i,:)) printf(rlabels(i,:)); else printf("%i", rlabels(i,:)); endif printf(" %10.3f", mat(i,:)); printf("\n"); endfor endfunction  NEW FILE: mle_example.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # 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);  NEW FILE: mle_results.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # report results function [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control) if nargin < 8 [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs); else [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control); endif V = mle_variance(theta, data, model, modelargs); # unscale results if argument has been passed # this puts coefficients into scale corresponding to the original modelargs if (nargin > 5) if iscell(unscale) # don't try it if unscale is simply a placeholder [theta, V] = unscale_parameters(theta, V, unscale); endif endif [theta, V] = delta_method("parameterize", theta, {data, model, modelargs}, V); n = rows(data); k = rows(V); se = sqrt(diag(V)); if convergence == 1 convergence="Normal convergence"; elseif convergence == 2 convergence="No convergence"; elseif convergence == 1 convergence = "Max. iters. exceeded"; endif printf("\n\n******************************************************\n"); disp(title); printf("\nMLE Estimation Results\n"); printf("BFGS convergence: %s\n\n", convergence); printf("Average LogL: %f\n", obj_value); printf("Observations: %d\n", n); a =[theta, se, theta./se, 2  2*normal_cdf(abs(theta ./ se))]; clabels = str2mat("estimate", "st. err", "tstat", "pvalue"); printf("\n"); if names !=0 prettyprint(a, names, clabels); else prettyprint_c(a, clabels); endif printf("\nInformation Criteria \n"); caic = 2*n*obj_value + rows(theta)*(log(n)+1); bic = 2*n*obj_value + rows(theta)*log(n); aic = 2*n*obj_value + 2*rows(theta); infocrit = [caic, bic, aic]; printf("CAIC : %8.4f Avg. CAIC: %8.4f\n", caic, caic/n); printf(" BIC : %8.4f Avg. BIC: %8.4f\n", bic, bic/n); printf(" AIC : %8.4f Avg. AIC: %8.4f\n", aic, aic/n); printf("******************************************************\n"); endfunction  NEW FILE: delta_method.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Computes Delta method mean and covariance of a nonlinear # transformation defined by "func" function [theta_transf, vartheta] = delta_method(func, theta, otherargs, vartheta) theta_transf = feval(func, theta, otherargs); D = numgradient(func, {theta, otherargs}); var_transf = D * vartheta * D'; endfunction  NEW FILE: unscale_parameters.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Unscales parameters that were estimated using scaled data # primarily for use by BFGS function [theta_us, vartheta_us] = unscale_parameters(theta, vartheta, scalecoefs); k = rows(theta); A = nth(scalecoefs, 1); b = nth(scalecoefs, 2); kk = rows(b); B = zeros(kk1,kk); B = [b'; B]; C = A + B; # allow for parameters that aren't associated with x if (k > kk) D = zeros(kk, k  kk); C = [C, D; D', eye(k  kk)]; endif; theta_us = C*theta; vartheta_us = C * vartheta * C'; endfunction  NEW FILE: mle_obj.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # this takes a general model and calculates average log likelihood function [obj_value, score] = mle_obj(theta, data, model, modelargs) [obj_value, score] = feval(model, theta, data, modelargs); obj_value =  mean(obj_value); # let's bulletproof this in case the model goes nuts if (((abs(obj_value) == Inf))  (isnan(obj_value))) obj_value = realmax; endif if isnumeric(score) score =  mean(score)'; endif endfunction  NEW FILE: parameterize.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # # This is an empty function, provided so that # delta_method will work as is. Replace it with # the parameter transformations your models use. # Note: you can let "otherargs" contain the model # name so that this function can do parameterizations # for a variety of models function theta = parameterize(theta, otherargs) endfunction  NEW FILE: mle_estimate.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # call the minimizing routine function [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control) if nargin == 3 [theta, obj_value, convergence] = bfgsmin("mle_obj", {theta, data, model, modelargs}); else [theta, obj_value, convergence] = bfgsmin("mle_obj", {theta, data, model, modelargs}, control); endif obj_value =  obj_value; # recover from minimization rather than maximization endfunction  NEW FILE: prettyprint_c.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # this prints matrices with column labels but no row labels function prettyprint_c(mat, clabels) printf(" "); # print the column labels clabels = [" ";clabels]; # pad to 8 characters wide clabels = strjust(clabels,"right"); k = columns(mat); for i = 1:k printf("%s ",clabels(i+1,:)); endfor # now print the row labels and rows printf("\n"); k = rows(mat); for i = 1:k printf(" %8.3f", mat(i,:)); printf("\n"); endfor endfunction  NEW FILE: mle_variance.m  # Copyright (C) 2003,2004 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # Copyright (C) 2003 Michael Creel michael.creel@... # 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 021111307 USA # sandwich form of varcov matrix function [V,scorecontribs,J_inv] = mle_variance(theta, data, model, modelargs) scorecontribs = numgradient(model, {theta, data, modelargs}); n = rows(scorecontribs); I = scorecontribs'*scorecontribs / n; J = numhessian("mle_obj", {theta, data, model, modelargs}); J_inv = inverse(J); V = J_inv*I*J_inv/n; endfunction 