rotatemat.ado

来自「是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到」· ADO 代码 · 共 1,657 行 · 第 1/3 页

*! version 1.0.0  06apr2005
program rotatemat, rclass
	version 9

	#del ;
	syntax  anything(id=matrix name=A)
	[,
	    // ----------------------- display options
	     
		BLanks(real 0)         // display small loadings as blanks
		noDISPlay              // suppresses all output
		noLOADing              // suprresses display rotated loadings
		noROTAtion             // supresses display rotation matrix
		FORmat(str)            // display format for matrices
		MATName(string)        // name for matrix
		COLNames(string)       // name for columns of matrix (plural)
		
	    //  ---------------------- minimization options
	    
		INIT(string)           // initial rotation matrix (dflt = I)
		RANdom                 // use random initial matrix (dflt = I)
		LOG                    // display iteration log
		TRAce                  // display trace 
		ITERate(passthru)      // max number of iterations
		PROTect(str)           // #trials (protect agnst local conv)
		TOLerance(passthru)    // tolerance in rotation matrix T
		LTOLerance(passthru)   // tolerance for optim criterion 
		GTOLerance(passthru)   // tolerance for projected gradient
		MAXSTEP(passthru)      // max #step halvings within iter
		MAXABSCORR(passthru)   // max abs(correlation), oblique only
		                       //
		PERturb(passthru)      // undocumented
		
	    //  ---------------------- rotation class
	    
		NORMalize              // apply horst normalization 
		Horst                  // synonym for normalize
	  	OBLIQue                // normal rotations
		ORTHOgonal             // orthonormal rotations
		
	    //  ---------------------- rotation criterion
	    
		BENtler
	  	CF(str)		       // no default argument with cf()
	  	ENTRopy
	  	EQUAmax
	  	OBLIMAx
	  	OBLIMIn                // short for oblimin(0)
	  	OBLIMInn(string)
	  	PARSImax
	  	PARTIAL(string)
	  	Promax                 // short for promax(3)
	  	Promaxn(string)        // abbrev for backward compatibility
	  	QUARTIMAx
	  	QUARTIMIn
	  	TANDEM1
	  	TANDEM2
	  	TARget(string)
	  	Varimax                // abbrev for backward compatibility
	  	VGPF
	] ;
	#del cr

// check/parse input

	if `"`matname'"' == "" {
		local matname `"`A'"'
	}
	
	if `"`colnames'"' == "" {
		local colnames "factors"
	}
	
	if `"`format'"' != "" { 
		local junk : display `format' 0.3
	}
	else {
		local format %9.5f     // default mentioned in help file
	}	

	confirm matrix `A'
	local p = rowsof(`A')   // number of rows/variables
	local k = colsof(`A')   // number of columns/factors
	if (`p' < 2) | (`k' < 1) {
		dis as err "too few variables or factors"
		exit 503 
	}
	if `p' < `k' {
		dis as err /// 
		    "number of variables should exceed number of factors" 
		exit 503
	}	
	
	NoMiss   `A' `"`matname'"'
	FullRank `A' `"`matname'"'

	// row-wise normalization ?
	if "`normalize'" != "" {
		if "`horst'" != "" { 
			dis as err "options horst and normalize are synonyms"
			exit 198
		}
		local horst horst
	}

	// gpopt := options for GPFcmd
	local gp_opt  `maxstep'   `maxabscorr' `iterate'  `tolerance' ///
	              `gtolerance' `ltolerance'

	// options for protected-optimization
	local nprotect  `protect'
	local protect   protect(`protect')
	local pargs     `protect' `perturb' `oblique' `log' `trace'
	if (`k'==1) & ("`nprotect'"!="") {
		dis as txt "(protect() ignored)"
		local nprotect
		local protect
		local pargs
	}	

	Check_args , `pargs' `gp_opt'

// specification of rotation class 

	if "`orthogonal'" != "" & "`oblique'" != "" {
		opts_exclusive "orthogonal oblique"
	}
	
// specification of rotation criterion

	if "`promax'" != "" & "`promaxn'" != "" {
		dis as err "promax and promax() are exclusive"
		exit 198
	}
	else if "`promax'" != "" { 
		local promaxn = 3
	}
	
	if "`oblimin'" != "" & "`obliminn'" != "" { 
		dis as err "oblimin and oblimin() are exclusive"
		exit 198
	}
	else if "`oblimin'" != "" { 
		local obliminn = 0
	}	

	local critn `bentler' `entropy' `equamax' `oblimax' `parsimax' ///
	   `quartimax' `quartimin' `tandem1' `tandem2' `varimax' `vgpf' 
	   
	if ("`cf'"       != "")   local critn `critn' cf
	if ("`obliminn'" != "")   local critn `critn' oblimin
	if ("`partial'"  != "")   local critn `critn' partial
	if ("`target'"   != "")   local critn `critn' target

	if "`promax'`promaxn'" != "" {
		local critn `critn' promax
		if "`orthogonal'" != "" {
			dis as err "promax is always oblique"
			exit 198
		}
		
		local oblique   // NOTE: later switched on again
		
		if "`promaxn'" != "" {
			confirm number `promaxn'
			if `promaxn' < 1 {
				dis as err "promax() should be at least 1"
				exit 198
			}
			else if `promaxn' > 4 {
				dis as txt "(promax() > 4 is ill advised)"
			}
		}
		else {
			local promaxn = 3
		}
	}

	if `"`critn'"' == "" {
		local critn   "varimax"
		local varimax "varimax"
		// dis as txt "(orthogonal varimax rotation assumed)"
	}
	
	if `:list sizeof critn' > 1 {
		dis as err "multiple rotation criteria specified"
		exit 198
	}

// INTERNAL CODE //////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

if inlist("`critn'","varimax","promax") {
	tempname AT f Gq iter rc T

	scalar `f' = .
	scalar `iter' = .
	
	local mxopts `init' `random' `nprotect' `log' `trace' `iterate' /// 
	        `tolerance' `ltolerance' `gtolerance' `maxstep' `maxabscorr'
	local nmax : word count `mxopts'
	if `nmax' > 0 {
		dis as error `"{p}maximize options for rotatemat do not "' ///
		    `"apply for the varimax or promax criterion{p_end}"'
		exit 198
	}
	
	if "`critn'" == "varimax" {
		NoOblique varimax `oblique'
		
		capture noisily {
			_rotate_mat `A' `AT' `T' , /// 
			   varimax `horst' 
		}
		scalar `rc' = _rc
		if colsof(`AT') < colsof(`A') { 
			dis as err "degenerate solution for varimax rotation"
			exit 198
		}
		
		local ctitle   "varimax"
		local ctitle12 "varimax"
		local class    "orthogonal"
	}
	else {
		// promax
		capture noisily { 
			_rotate_mat `A' `AT' `T' , /// 
			   promax(`promaxn') `horst' 
		}	
		scalar `rc' = _rc
		if colsof(`AT') < colsof(`A') { 
			dis as err "degenerate solution for promax rotation"
			exit 198
		}
		
		// _rotate_mat returns the varimax rotation
		// we solve for T in AT = A*inv(T')
		// note that T is not -normal- 
		
		matrix `T' = (`A''*`A')* inv(`AT''*`A')  
		
		local ctitle   "promax(`promaxn')"
		local ctitle12 "promax(`promaxn')"
		if length("`ctitle12'") > 12 {
			local ctitle12 "promax(..)"
		}
		local class   "oblique"
		local oblique oblique 
		local carg    `promaxn'
	}
}

// GRADIENT PROJECTION ALGORITHM //////////////////////////////////////////////
// code is indented to the left; all comments are prededed with [GPF]
///////////////////////////////////////////////////////////////////////////////

else {

//[GPF] initialize rotation T
	if "`init'"!="" & "`random'"!="" {
		di as err "only one of init() or random are allowed"
		exit 198
	}
	if "`nprotect'"!="" & "`init'"!="" {
		di as error "init() cannot be used with protect(); " ///
		   "protect() implies random"
		exit 198
	}

	tempname T
	Init `T' `A' "`init'" "`random'" "`oblique'"
	
//[GPF] normalization A --> L D

	if "`horst'" != "" {
		tempname L D
		Normalize `A' `L' `D'
	}
	else {
		local L `A'
	}

//[GPF] minimze orthogonal/oblique rotation L over suitable class of rotations

	// default criterion title, may be overwritten below
	local ctitle `critn'
	local ctitle12 `critn'
	if length("`ctitle12'") > 12 {
		local ctitle12 = substr("`ctitle12'",1,12)
	}
	
	if "`bentler'" != "" {
		local ctitle    "Bentler's invariant pattern simplicity"
		local ctitle12  "Bentler"
		local copt      "crit(Bentler)"
	}
	else if `"`cf'"' != "" {
		tempname carg
		capture scalar `carg' = `cf'
		if _rc {
			dis as err "argument of cf() invalid"
			exit 198
		}
		if `carg' < 0 | `carg' > 1 {
			dis as err "cf() must be between 0 and 1"
			exit 198
		}
		local ctitle    "Crawford-Ferguson(`cf')"
		local ctitle12  "cf(`cf')" 
		if length("`ctitle12'") > 12 {
			local ctitle12 "cf(..)"
		}
		local copt      "crit(Cf) critopt(`carg')"
		local carg = `carg'
	}
	else if "`entropy'" != "" {
		NoOblique entropy `oblique'
		local ctitle    "minimum entropy"
		local ctitle12  "min. entropy"
		local copt      "crit(Entropy)"  
	}
	else if "`equamax'" != "" {
		NoOblique equamax `oblique'
		local gam = `p'/2
		local copt "crit(Oblimin) critopt(`gam')"
	}
	else if "`oblimax'" != "" {
		local copt "crit(Oblimax)" 
	}
	else if `"`obliminn'"' != "" {
		tempname carg
		capture scalar `carg' = `obliminn'
		if _rc {
			dis as err "argument of oblimin() invalid"
			exit 198
		}
		local ctitle  "oblimin(`obliminn')"
		local ctitle12 `ctitle'
		if length("`ctitle12'") > 12 {
			local ctitle12 "cf(..)"
		}
		local copt    "crit(Oblimin) critopt(`carg')" 
		local carg = `carg'
	}
	else if "`parsimax'" != "" {
		NoOblique parsimax `oblique'
	  	local kap = (`k'-1) / (`p'+`k'-2)
	  	local copt  "crit(Cf) critopt(`kap')" 
	}
	else if `"`partial'"' != "" {
		if `:word count `partial'' != 2 {
			dis as err "partial() expects two matrices"
			exit 198
		}
		gettoken H W : partial
		CheckPartial `L' `H' `W'
		local ctitle   "partial target (target=`H' weight=`W')"
		local ctitle12 "ptarget(..)"
		local carg     "`H' `W'"
		local copt     "crit(PartialTarget) critopt(`carg')"
	}
	else if "`quartimax'" != "" {
		NoOblique quartimax `oblique'
		local copt  "crit(Quartimax) "
	}
	else if "`quartimin'" != "" {
		local copt  "crit(Quartimin)"
	}
	else if "`tandem1'" != "" {
		NoOblique tandem1 `oblique'
		local ctitle   "Comrey's principle 1"
		local ctitle12 "tandem-1"
		local copt     "crit(Tandem1)"
	}
	else if "`tandem2'" != "" {
		NoOblique tandem2 `oblique'
		local ctitle   "Comrey's principle 2"
		local ctitle12 "tandem-2"
		local copt     "crit(Tandem2)"
	}
	else if `"`target'"' != "" {
		CheckTarget `A' `target'
		local ctitle   "target(`target')"
		local ctitle12 "target(..)"
		local carg     `target'
		local copt     "crit(Target) critopt(`carg')"
	}
	else if "`vgpf'" != "" {
		NoOblique varimax `oblique'
		local ctitle   varimax
		local ctitle12 varimax
		local critn    varimax
		local gam = 1
		local copt  "crit(Oblimin) critopt(`gam')"
	}
	else {
		_stata_internalerror rotatemat
	}

	
	local GPFcmd = cond("`oblique'"!="", "GPFoblique", "GPFortho")
		
	Minimize `"`GPFcmd' `L', `copt' `gp_opt' colnames(`colnames')"' ///
	         `T'   "`pargs'"

//[GPF] retrieve results

	tempname AT f fmin iter nnconv rc T

	local  class  "`r(class)'"
	scalar `iter' = r(iter)
	scalar `f'    = r(f)
	matrix `T'    = r(Th)
	matrix `AT'   = r(Lh)

	if "`nprotect'" != "" {
		scalar `nnconv' = r(nnconv)
		matrix `fmin'   = r(fmin)
	}
	scalar `rc' = r(rc)

//[GPF] invert normalization r(AT)

	if "`horst'" != "" {
		InvertNormalize `AT' `D'
	}
}
///// END-OF-GPF //////////////////////////////////////////////////////////////

// unique representation wrt direction and L1 norm	

	Order `AT' `T' 
	
// set column and row names

	matrix rownames `AT' = `:rownames `A''
	matrix colnames `AT' = `:colnames `A''
	matrix rownames  `T' = `:colnames `A''
	matrix colnames  `T' = `:colnames `A''

// display results ////////////////////////////////////////////////////////////

if "`display'" == "" {
	dis _n as txt `"Rotation of `matname'[`p',`k']"' _n

	dis as txt "    Criterion               " as res `"`ctitle'"'
	dis as txt "    Rotation class          " as res `"`class'"'
	dis as txt "    Horst normalization     " as res ///
	                                        cond("`horst'"!="","on","off")
	    
	if `f' != . {
		local lf : display %9.0g `f' 
		dis as txt "    Criterion value         " as res `lf' 
	}
	
	if `iter' != . {
		dis as txt "    Number of iterations    " as res `iter'
	}
	
	if `rc' != 0 {
		dis as txt "    Warning: convergence not achieved"
	}

	if "`loading'" == "" { 
		if `blanks' == 0 {
			dis _n as txt `"Rotated `colnames'"'
			matlist `AT', format(`format') border(row) left(4)
		}
		else {
			dis _n as txt `"Rotated `colnames' "' /// 
			              `"(blanks represent abs()<`blanks')"'
			     
			tempname AT0
			matrix `AT0' = `AT'
			_small2dotz `AT0' `blanks'
			matlist `AT0', /// 
			   format(`format') border(row) left(4) nodotz
			matrix drop `AT0'
		}
	}
	if "`rotation'" == "" {
		dis _n as txt proper("`class'") " rotation"
		matlist `T', format(`format') border(row) left(4)
	}	

	if "`nprotect'" != "" {
		dis
		if rowsof(`fmin') > 1 {
			dis as txt /// 
			    "    Warning: convergence to multiple minima"
			dis _n as txt "  criterion       freq" 
			mat list `fmin', nonames noheader
		}
		else {
			dis as txt "{pstd}criterion minimization from " ///
			   "`nprotect' random initial rotations converged " /// 
			   "to the same minimum{p_end}"
		}
	}
}

// return results ////////////////////////////////////////////////////////////

	return clear

	return scalar rc   = `rc'                // return code
	return scalar iter = `iter'              // number of GPF iterations
	return scalar f    = `f'                 // crit. f(L.T), L = norm(A)
	return matrix T    = `T'                 // optimal transformation T
	return matrix AT   = `AT'                // optimal AT = A.T

	if "`nprotect'" != "" {
		return scalar nnconv = `nnconv'  // #non-convergent trails
		return matrix fmin   = `fmin'    // minima found
	}

	return local ctitle    `"`ctitle'"'      // label of rotation method
	return local ctitle12  `"`ctitle12'"'    // title of at most 12 chars
	return local criterion `"`critn'"'       // criterion name (eg oblimin)
	return local class     `class'           // orthoginal or oblique
	return local carg      `carg'            // argument for criterion

	if "`horst'" != "" { 
		return local normalization horst
	}
	else {
		return local normalization none 
	}	

	return local cmd rotatemat
end

// ----------------------------------------------------------------------------
// criterion minimization commands
//
//    GPFortho   -- with respect to orthogonal matrices
//    GPFoblique -- with respect to normal matrices
//
// convergence is declared if all threse conditions are met:
//
//    npg   = norm(projected gradient)           < gtolerance
//    dcoef = mreldif(A(iter),A(iter-1))         < tolerance
//    dcrit = reldif(crit(iter), crit(iter-1))   < ltolerance
//
// return codes
//    rc = 0     convergence achieved
//    rc = 1001  stephalving failed
//    rc = 1002  maxiter reached
// ----------------------------------------------------------------------------

program GPFortho, rclass