Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

Commit [afe5a1] default Maximize Restore History

Altered the directory structure to match the package system

hauberg hauberg 2006-08-20

<< < 1 2 3 > >> (Page 2 of 3)
removed poisson_moments.m
removed poisson.m
removed average_moments.m
removed nls_obj.m
removed nls_obj_nodes.m
removed nls_estimate.m
removed prettyprint_c.m
removed gmm_example.m
removed mle_example.m
removed README
removed nls_example.m
removed sum_moments_nodes.m
removed mle_variance.m
removed delta_method.m
removed mle_obj.m
removed parameterize.m
removed prettyprint.m
copied gmm_variance.m -> inst/poisson.m
copied gmm_obj.m -> inst/mle_obj_nodes.m
copied gmm_variance_inefficient.m -> inst/gmm_estimate.m
copied mle_estimate.m -> inst/mle_results.m
poisson_moments.m
File was removed.
poisson.m
File was removed.
average_moments.m
File was removed.
nls_obj.m
File was removed.
nls_obj_nodes.m
File was removed.
nls_estimate.m
File was removed.
prettyprint_c.m
File was removed.
gmm_example.m
File was removed.
mle_example.m
File was removed.
README
File was removed.
nls_example.m
File was removed.
sum_moments_nodes.m
File was removed.
mle_variance.m
File was removed.
delta_method.m
File was removed.
mle_obj.m
File was removed.
parameterize.m
File was removed.
prettyprint.m
File was removed.
gmm_variance.m to inst/poisson.m
--- a/gmm_variance.m
+++ b/inst/poisson.m
@@ -1,4 +1,4 @@
-# Copyright (C) 2003,2004  Michael Creel <michael.creel@uab.es>
+# Copyright (C) 2005  Michael Creel <michael.creel@uab.es>
 # 
 # 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
@@ -13,13 +13,14 @@
 # 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 likelihood function (Poisson for count data) with score
 
-
-# GMM variance, which assumes weights are optimal
-function V = gmm_variance(theta, data, weight, moments, momentargs)
-	D = numgradient("average_moments", {theta, data, moments, momentargs});
-	D = D';
-	m = feval(moments, theta, data, momentargs); # find out how many obsns. we have
-	n = rows(m);
-	V = (1/n)*inv(D*weight*D');
-endfunction
+function [log_density, score] = poisson(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);
+	if (otherargs{1} == 1) score = "na"; endif # provide analytic score or not?
+endfunction	
gmm_obj.m to inst/mle_obj_nodes.m
--- a/gmm_obj.m
+++ b/inst/mle_obj_nodes.m
@@ -1,32 +1,34 @@
-# Copyright (C) 2003,2004,2005  Michael Creel <michael.creel@uab.es>
-# 
-# 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,2004,2005  Michael Creel <michael.creel@uab.es>
+##
+## 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
 
-
-# The GMM objective function, for internal use by gmm_estimate
-# This is scaled so that it converges to a finite number.
-# To get the chi-square specification
-# test you need to multiply by n (the sample size)
-function obj_value = gmm_obj(theta, data, weight, moments, momentargs)
-
-	m = average_moments(theta, data, moments, momentargs);
+function contrib = mle_obj_nodes(theta, data, model, modelargs, nn)
+	global NEWORLD NSLAVES
 	
-	obj_value = m' * weight *m;
-
-	if (((abs(obj_value) == Inf)) || (isnan(obj_value)))
-		obj_value = realmax;
+	# Who am I?
+	[info, rank] = MPI_Comm_rank(NEWORLD); 
+	if rank == 0 # Do this if I'm master
+		startblock = NSLAVES*nn + 1;
+		endblock = rows(data);
+	else	# this is for the slaves
+		startblock = rank*nn-nn+1;
+		endblock = rank*nn;
 	endif	
 
-endfunction	
+	data = data(startblock:endblock,:);
+	contrib = feval(model, theta, data, modelargs);
+	contrib = sum(contrib);
+
+endfunction
gmm_variance_inefficient.m to inst/gmm_estimate.m
--- a/gmm_variance_inefficient.m
+++ b/inst/gmm_estimate.m
@@ -1,28 +1,76 @@
-# Copyright (C) 2003,2004  Michael Creel <michael.creel@uab.es>
-#
+# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es>
+# 
 # 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
 
-# GMM variance, which assumes weights are not optimal
-function V = gmm_variance_inefficient(theta, data, weight, omega, moments, momentargs)
-	D = numgradient("average_moments", {theta, data, moments, momentargs});
-	D = D';
+# usage: [theta, obj_value, convergence, iters] = 
+#           gmm_estimate(theta, data, weight, moments, momentargs, control, nslaves)
+#
+# inputs:
+#      theta: column vector initial parameters
+#       data: data matrix
+#     weight: the GMM weight matrix
+#    moments: name of function computes the moments
+#	      (should return nXg matrix of contributions)
+# momentargs: (cell) additional inputs needed to compute moments.
+# 	      May be empty ("")
+#    control: (optional) BFGS or SA controls (see bfgsmin and samin).
+#             May be empty (""). 
+#    nslaves: (optional) number of slaves if executed in parallel
+#             (requires MPITB)
+#
+# outputs:
+# theta: GMM estimate of parameters
+# obj_value: the value of the gmm obj. function
+# convergence: return code from bfgsmin
+#              (1 means success, see bfgsmin for details)
+# iters: number of BFGS iteration used
+#
+# please type "gmm_example" while in octave to see an example
 
-	n = rows(data);
+# call the minimizing routine
+function [theta, obj_value, convergence, iters] = gmm_estimate(theta, data, weight, moments, momentargs, control, nslaves)
 
-	J = D*weight*D';
-	J = inv(J);
-	I = D*weight*omega*weight*D';
-	V = (1/n)*J*I*J;
-endfunction
+	if nargin < 5 error("gmm_estimate: 5 arguments required"); endif
+
+	if nargin < 6 control = {Inf,0,1,1}; endif # default controls
+	if !iscell(control) control = {Inf,0,1,1}; endif # default controls if receive placeholder
+	if nargin < 7 nslaves = 0; endif
+
+ 	if nslaves > 0
+		global NSLAVES PARALLEL NEWORLD NSLAVES TAG;
+		LAM_Init(nslaves);
+		# Send the data to all nodes
+		NumCmds_Send({"data", "weight", "moments", "momentargs"}, {data, weight, moments, momentargs});
+	endif
+
+	# bfgs or sa?
+	if (size(control,1)*size(control,2) == 0) # use default bfgs if no control
+		control = {Inf,0,1,1};	
+	  method = "bfgs";
+	elseif (size(control,1)*size(control,2) < 11)
+		method = "bfgs";
+	else method = "sa";	
+	endif
+	
+
+	if strcmp(method, "bfgs")
+		[theta, obj_value, convergence, iters] = bfgsmin("gmm_obj", {theta, data, weight, moments, momentargs}, control);
+	elseif strcmp(method, "sa") 	
+	  	[theta, obj_value, convergence] = samin("gmm_obj", {theta, data, weight, moments, momentargs}, control);
+	endif	
+
+	if nslaves > 0 LAM_Finalize; endif # clean up		
+
+endfunction	
mle_estimate.m to inst/mle_results.m
--- a/mle_estimate.m
+++ b/inst/mle_results.m
@@ -1,27 +1,30 @@
-## Copyright (C) 2003,2004,2005  Michael Creel <michael.creel@uab.es>
-##
-## 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,2004,2005  Michael Creel <michael.creel@uab.es>
+#  
+# 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
 
-# usage: 
-# [theta, obj_value, conv, iters] = mle_estimate(theta, data, model, modelargs, control, nslaves) 
+# usage: [theta, V, obj_value, infocrit] = 
+#    mle_results(theta, data, model, modelargs, names, title, unscale, control)
 #
 # inputs:
 # theta: column vector of model parameters
 # data: data matrix
 # model: name of function that computes log-likelihood
 # modelargs: (cell) additional inputs needed by model. May be empty ("")
+# names: vector of parameter names, e.g., use names = str2mat("param1", "param2");
+# title: string, describes model estimated
+# unscale: (optional) cell that holds means and std. dev. of data (see scale_data) 
 # control: (optional) BFGS or SA controls (see bfgsmin and samin). May be empty (""). 
 # nslaves: (optional) number of slaves if executed in parallel (requires MPITB)
 #
@@ -30,45 +33,61 @@
 # obj_value: the value of the log likelihood function at ML estimate
 # conv: return code from bfgsmin (1 means success, see bfgsmin for details)
 # iters: number of BFGS iteration used
-#
-# please see mle_example.m for examples of how to use this
-function [theta, obj_value, convergence, iters] = mle_estimate(theta, data, model, modelargs, control, nslaves)
 
 
-	if nargin < 3
-		error("mle_estimate: 3 arguments required");
+##
+## Please see mle_example for information on how to use this
+ 
+# report results
+function [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control, nslaves)
+	if nargin < 9 nslaves = 0; endif # serial by default
+	if nargin < 8
+		[theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, "", nslaves);
+	else
+		[theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control, nslaves);
+	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 > 6)
+    		if iscell(unscale) # don't try it if unscale is simply a placeholder
+			[theta, V] = unscale_parameters(theta, V, unscale);
+    		endif
 	endif
 
-	if nargin < 4 modelargs = {}; endif # create placeholder if not used
-	if !iscell(modelargs) modelargs = {}; endif # default controls if receive placeholder
-	if nargin < 5 control = {Inf,0,1,1}; endif # default controls and method
-	if !iscell(control) control = {Inf,0,1,1}; endif # default controls if receive placeholder
-	if nargin < 6 nslaves = 0; endif
-	if nslaves > 0
-		global NSLAVES PARALLEL NEWORLD TAG;
-		LAM_Init(nslaves);
-		# Send the data to all nodes
-		NumCmds_Send({"data", "model", "modelargs"}, {data, model, modelargs});
+	[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 Log-L: %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", "t-stat", "p-value");
+
+	printf("\n");
+	if names !=0 prettyprint(a, names, clabels);
+	else prettyprint_c(a, clabels);
 	endif
 
-	# bfgs or sa?
-	if (size(control,1)*size(control,2) == 0) # use default bfgs if no control
-		control = {Inf,0,1,1};	
-	  method = "bfgs";
-	elseif (size(control,1)*size(control,2) < 11)
-		method = "bfgs";
-	else method = "sa";	
-	endif
-	
-
-	if strcmp(method, "bfgs")
-	  [theta, obj_value, convergence, iters] = bfgsmin("mle_obj", {theta, data, model, modelargs, nslaves}, control);
-	elseif strcmp(method, "sa") 	
-	  [theta, obj_value, convergence] = samin("mle_obj", {theta, data, model, modelargs, nslaves}, control);
-	endif	
-		
-	if nslaves > 0
-		LAM_Finalize;
-	endif # cleanup			
-	obj_value = - obj_value; # recover from minimization rather than maximization
-endfunction	
+	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
<< < 1 2 3 > >> (Page 2 of 3)