Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: pointing to a class member function (likelihood) with optimize() in mata


From   Pierre Carl Michaud <pcmichaud@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: pointing to a class member function (likelihood) with optimize() in mata
Date   Thu, 28 Oct 2010 11:40:45 -0400

Thanks. I think I was working along similar lines although I had not yet succeeded. I was thinking of decoupling the model class from the estimation class. If I understood correctly, the following could work following what you wrote

class model {
	specifies model
	provides eval function
}

class estimate {
	starts an instance of model, say m
	sets parameters in m
	has the wrapper eval calling m.eval
	uses optimize using &eval for formats output, covar, etc
}

Is that correct?

PC



Le 2010-10-28 à 11:19, Jeff Pitblado, StataCorp LP a écrit :

> Pierre-Carl Michaud <pcmichaud@gmail.com> is trying to use -optimize()- with
> an evaluator function implemented as a member function of a Mata class:
> 
>> In mata, I am trying to write a general class to estimate models. I  
>> have had no problem with optimize outside of a class. I want now to  
>> call optimize in the class and I have a problem with the line
>>   optimize_init_evaluator(s, &func())
>> The likelihood func() is a member function (method) in the class and  
>> as such referring to &func() does not find it. I've tried  
>> &classname::func() and no luck either. &this.func() does refer to  
>> something but then I get
>>   modelbase::func():  3001  expected 6 arguments but received 1
>> same error message if I define func() in the class definition as a  
>> static member function.
>> Anyone has an idea of how to point to a class function so that  
>> optimize() finds it and its arguments?
> 
> There is no mechanism in object oriented programming (OOP) for calling a
> class's public member function using it's address.  The user-programmer would
> be required to know too many things about how the class system is implemented,
> and the system would have to provide hooks to even make it possible.
> 
> Pierre-Carl can still use his general class, he will just need to provide
> -optimize()- with a simple wrap-around function that basically is a
> call-through to the class's evaluator function.  Here is a simple example,
> using the logistic regression model.
> 
> First, let's define a class called 'logit_model' that has two private data
> members, and three public function members.
> 
> 	class logit_model {
> 	private:
> 		real colvector	y
> 		real matrix	X
> 
> 	public:
> 		void set_y()
> 		void set_X()
> 		void eval()
> 	}
> 
> The 'set_y()' function simply provides a way for the programmer to supply the
> class instance with a dependent variable; similarly, 'set_X()' provides a way
> to supply the independent variables.  The 'eval()' function will compute the
> log-likelihood, gradient, and Hessian for the logistic regression model.
> 
> The 'set_*()' functions simply make a copy of the supplied argument.
> 
> 	void logit_model::set_y(real colvector y)
> 	{
> 		this.y = y
> 	}
> 
> 	void logit_model::set_X(real matrix X)
> 	{
> 		this.X = X
> 	}
> 
> The 'eval()' function computes a log-likelihood value, gradient vector , and
> Hessian matrix according to the requirements of -optimize()-'s "d2" evaluator
> functions:
> 
> 	void logit_model::eval(	real	scalar		todo,
> 				real	rowvector	beta,
> 				real	scalar		lnf,
> 				real	rowvector	g,
> 				real	matrix		H)
> 	{
> 		real colvector	pm
> 		real colvector	xb
> 		real colvector	lj
> 		real colvector	dllj
> 		real colvector	d2llj
> 
> 		pm	= 2*(this.y :!= 0) :- 1
> 		xb	= this.X*beta'
> 
> 		lj	= invlogit(pm:*xb)
> 		if (any(lj :== 0)) {
> 			lnf = .
> 			return
> 		}
> 		lnf = quadcolsum(ln(lj))
> 		if (todo == 0) return
> 
> 		dllj	= pm :* invlogit(-pm:*xb)
> 		if (missing(dllj)) {
> 			lnf = .
> 			return
> 		}
> 		g	= quadcross(dllj, X)
> 		if (todo == 1) return
> 
> 		d2llj	= abs(dllj) :* lj
> 		if (missing(d2llj)) {
> 			lnf = .
> 			return
> 		}
> 		H	= - quadcross(X, d2llj, X)
> 	}
> 
> The only other function we need is the wrap-around function.  This wrap-around
> function will take the arguments of -optimize()-'s "d2" evaluators, with the
> addition of a single argument for the class instance.  We named the
> wrap-around function 'logit_eval()', and here it is:
> 
> 	void logit_eval(	real	scalar		todo,
> 				real	rowvector	beta,
> 				class	logit_model	M,
> 				real	scalar		lnf,
> 				real	rowvector	g,
> 				real	matrix		H)
> 	{
> 		M.eval(todo,beta,lnf,g,H)
> 	}
> 
> Here is an example using the auto data, showing how this class and function
> are used.
> 
> 	sysuse auto
> 	gen one = 1
> 
> First we need to declare an instance of our class, and supply it with some
> data for 'y' and 'X'.
> 
> 	mata
> 	u
> 	M = logit_model()
> 	M.set_y(st_data(., "foreign"))
> 	M.set_X(st_data(., "mpg turn trunk one"))
> 
> Now we can work with -optimize()-.
> 
> First we initialize an -optimize()- handle,
> 
> 	S = optimize_init()
> 
> then declare the wrap-around evaluator function and its type,
> 
> 	optimize_init_evaluator(S, &logit_eval())
> 	optimize_init_evaluatortype(S, "d2")
> 
> then declare our class instance as an extra argument for the evaluator,
> 
> 	optimize_init_argument(S, 1, M)
> 
> provide starting values, and finally perform the optimization.
> 
> 	optimize_init_params(S, J(1,4,0))
> 	optimize(S)
> 	end
> 
> Here is a log of the above example.  We added a call to Stata's -logit-
> command to show that the results from -optimize()- match.
> 
> ***** BEGIN:
> . sysuse auto
> (1978 Automobile Data)
> 
> . gen one = 1
> 
> . 
> . mata:
> ------------------------------------------------- mata (type end to exit) -----
> : 
> : class logit_model {
>> private:
>>        real colvector  y
>>        real matrix     X
>> 
>> public:
>>        void set_y()
>>        void set_X()
>>        void eval()
>> }
> 
> : 
> : void logit_model::set_y(real colvector y)
>> {
>>        this.y = y
>> }
> 
> : 
> : void logit_model::set_X(real matrix X)
>> {
>>        this.X = X
>> }
> 
> : 
> : void logit_model::eval( real    scalar          todo,
>>                        real    rowvector       beta,
>>                        real    scalar          lnf,
>>                        real    rowvector       g,
>>                        real    matrix          H)
>> {
>>        real colvector  pm
>>        real colvector  xb
>>        real colvector  lj
>>        real colvector  dllj
>>        real colvector  d2llj
>> 
>>        pm      = 2*(this.y :!= 0) :- 1
>>        xb      = this.X*beta'
>> 
>>        lj      = invlogit(pm:*xb)
>>        if (any(lj :== 0)) {
>>                lnf = .
>>                return
>>        }
>>        lnf = quadcolsum(ln(lj))
>>        if (todo == 0) return
>> 
>>        dllj    = pm :* invlogit(-pm:*xb)
>>        if (missing(dllj)) {
>>                lnf = .
>>                return
>>        }
>>        g       = quadcross(dllj, X)
>>        if (todo == 1) return
>> 
>>        d2llj   = abs(dllj) :* lj
>>        if (missing(d2llj)) {
>>                lnf = .
>>                return
>>        }
>>        H       = - quadcross(X, d2llj, X)
>> }
> 
> : 
> : void logit_eval(        real    scalar          todo,
>>                        real    rowvector       beta,
>>                        class   logit_model     M,
>>                        real    scalar          lnf,
>>                        real    rowvector       g,
>>                        real    matrix          H)
>> {
>>        M.eval(todo,beta,lnf,g,H)
>> }
> 
> : 
> : M = logit_model()
> 
> : M.set_y(st_data(., "foreign"))
> 
> : M.set_X(st_data(., "mpg turn trunk one"))
> 
> : 
> : S = optimize_init()
> 
> : optimize_init_evaluator(S, &logit_eval())
> 
> : optimize_init_evaluatortype(S, "d2")
> 
> : optimize_init_argument(S, 1, M)
> 
> : optimize_init_params(S, J(1,4,0))
> 
> : optimize(S)
> Iteration 0:   f(p) = -51.292891  
> Iteration 1:   f(p) = -25.856762  
> Iteration 2:   f(p) = -25.685523  
> Iteration 3:   f(p) = -25.685148  
> Iteration 4:   f(p) = -25.685148  
>                  1              2              3              4
>    +-------------------------------------------------------------+
>  1 |  -.0776544663   -.6146979766   -.0214041449    24.53154516  |
>    +-------------------------------------------------------------+
> 
> : 
> : end
> -------------------------------------------------------------------------------
> 
> . 
> . logit for mpg turn trunk
> 
> Iteration 0:   log likelihood =  -45.03321  
> Iteration 1:   log likelihood = -28.007564  
> Iteration 2:   log likelihood = -25.802293  
> Iteration 3:   log likelihood = -25.685271  
> Iteration 4:   log likelihood = -25.685148  
> Iteration 5:   log likelihood = -25.685148  
> 
> Logistic regression                               Number of obs   =         74
>                                                  LR chi2(3)      =      38.70
>                                                  Prob > chi2     =     0.0000
> Log likelihood = -25.685148                       Pseudo R2       =     0.4296
> 
> ------------------------------------------------------------------------------
>     foreign |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
> -------------+----------------------------------------------------------------
>         mpg |  -.0776545   .0708997    -1.10   0.273    -.2166154    .0613065
>        turn |   -.614698   .1599212    -3.84   0.000    -.9281377   -.3012583
>       trunk |  -.0214041   .1120528    -0.19   0.849    -.2410236    .1982153
>       _cons |   24.53155   6.785011     3.62   0.000     11.23317    37.82992
> ------------------------------------------------------------------------------
> ***** END:
> 
> --Jeff					-- Hua
> jpitblado@stata.com			hpeng@stata.com
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index