_biplotmat.ado

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

*! version 1.0.0  25mar2005
program _biplotmat, rclass
	version 9.0

	syntax anything(id=matrix name=m)	///
	[,					///
		ALPha(numlist max=1 >=0 <=1)	///
		MULTANDDIV(numlist max=1)	/// undocumented
		STretch(numlist max=1 >0)	///
		TABle 				///
		SEParate			///
		XNEGate YNEGate 		///
		AUTOaspect			///
		ASPECTratio(str)		///
		noGraph				///
		touse(varname)			///
		*				///
	]
	
	if "`separate'" != "" {
		local getcombine getcombine
	}

	_parse combop options : options, option(ROWopts) opsin
	_parse combop options : options, option(COLopts) opsin
	_get_gropts , graphopts(`options') 			///
		getallowed(ROWopts COLopts) `getcombine'
		
	local goptions   `"`s(graphopts)'"'
	local coloptions `"`s(colopts)'"'
	local rowoptions `"`s(rowopts)'"'
	local gcopts   `"`s(combineopts)'"'
	
	if "`aspectratio'" != "" & "`autoaspect'" != "" {
		display as error				///
		    "options aspectratio() and autoaspect may not be combined"
		exit 198
	}

	_parse combop rowoptions : rowoptions, option(NAME) opsin rightmost
	_parse combop coloptions : coloptions, option(NAME) opsin rightmost

	local 0 , `coloptions'
	syntax , [ noLABel NAME(str) MLABPosition(passthru) 		///
		   MLABVPosition(passthru) MLabel(str) * ]
	local cname `name'
	local nolabel_cols `label'
	local coloptions `options' `mlabposition' `mlabvposition'
	
	if "`mlabel'" != "" {
		display as error "option mlabel() is not allowed with colopts()"
		exit 198
	}
	if "`mlabposition'" != "" | "`mlabvposition'" != "" {
		local colposition_flag "yes"
	}

	local 0 , `rowoptions'
	syntax , [ noLABel NAME(str) MLABPosition(passthru) 		///
		   MLABVPosition(passthru) MLabel(str) * ]
	local rname `name'
	local nolabel_rows `label'
	local rowoptions `options' `mlabposition' `mlabvposition'

	if "`mlabel'" != "" {
		local rowoptions `rowoptions' mlabel(`mlabel')
		local rlabvar `mlabel'
	}
	if "`mlabposition'" != "" | "`mlabvposition'" != "" {
		local rowposition_flag "yes"
	}
	
	.atk = .aspect_axis_toolkit.new
	.atk.setPreferredLabelCount 7
	
	_parse comma aspect_ratio placement : aspectratio
	if "`aspect_ratio'" != "" {
		confirm number `aspect_ratio'
		.atk.setPreferredAspect `aspect_ratio'
		.atk.setShowAspectTrue
	}
	
	if "`autoaspect'" != "" {
		.atk.setAutoAspectTrue
	}	

// check input is valid

	confirm matrix `m'
	local nr = rowsof(`m')
	local nc = colsof(`m')
	if `nr'<2 | `nc'<2 {
		dis as err "biplot requires at least 2 rows and columns"
		exit 503
	}

	if "`multanddiv'" != "" {
		if `multanddiv' == 0 {
			di as error "multanddiv() may not be zero"
			exit 198
		}
	}

	if "`alpha'" == "" {
		local alpha = 0.5
	}

	if (`"`rname'"'  == "")  local rname   Rows
	if (`"`cname'"'  == "")  local cname   Columns

// create biplot matrices

	tempname U V Vhold W rho rho1 rho2 w2
	if `nr' > `nc' {
		matrix svd `U' `W' `V' = `m'
	}
	else {
		tempname mt
		matrix `mt' = `m''
		matrix svd `V' `W' `U' = `mt'
		matrix drop `mt'
	}

	// singular values are not sorted.
	// select largest sv with index i1 and second largest with index i2

	local nw = colsof(`W')
	local i1  = 1
	local ni1 = 1
	forvalues i = 2 / `nw' {
		if `W'[1,`i'] > `W'[1,`i1'] {
			local i1  = `i'
			local ni1 = 1
		}
		else if `W'[1,`i'] == `W'[1,`i1'] {
			local ++ni1
		}
	}

	local i2  = 1 + (`i1'==1)
	local ni2 = 1
	forvalues i = 1 / `nw' {
		if inlist(`i',`i1',`i2') continue

		if (`W'[1,`i'] > `W'[1,`i2']) {
			local i2  = `i'
			local ni2 = 1
		}
		else if `W'[1,`i'] == `W'[1,`i2'] {
			local ++ni2
		}
	}

	if (`ni1' > 1) | (`ni2' > 1) {
		display as txt "(largest and second largest singular values not uniquely defined)"
		matlist `W', nonames
	}

	// For each column of U (and associated column of V) it is arbitrary
	// if it is +1 or -1 times it (i.e. sign flip of the column).  To
	// minimize differences across different platforms / runs, we take the
	// first nonzero element (starting from the top of a column of U) and
	// make it positive.  We arbitrarily use a threshold of 1e-10 for
	// determining if something is zero or not for this exercise (the
	// threshold makes little difference).  It should be very rare for
	// sign flip differences to appear between platforms using this rule
	// (e.g. first element is very close to our arbitrary threshold and
	// lands either side of it on the 2 platforms).

	// check it for column i1
	forvalues i = 1 / `=rowsof(`U')' { // rare to loop over more than a few
		if abs(`U'[`i',`i1']) > 1e-10 {
			if `U'[`i',`i1'] < 0 {
				local flip1 "-1*"
			}
			continue, break
		}
	}

	// check it for column i2
	forvalues i = 1 / `=rowsof(`U')' { // rare to loop over more than a few
		if abs(`U'[`i',`i2']) > 1e-10 {
			if `U'[`i',`i2'] < 0 {
				local flip2 "-1*"
			}
			continue, break
		}
	}

// select U2 and V2 with columns (i1,i2)

	matrix `U' = `flip1'`W'[1,`i1']^( `alpha' )*`U'[1...,`i1'] , ///
		     `flip2'`W'[1,`i2']^( `alpha' )*`U'[1...,`i2']

	matrix `V' = `flip1'`W'[1,`i1']^(1-`alpha')*`V'[1...,`i1'] , ///
		     `flip2'`W'[1,`i2']^(1-`alpha')*`V'[1...,`i2']

	matrix colnames `U' = dim1 dim2
	matrix colnames `V' = dim1 dim2

	local Unames : rownames `U'
	local Vnames : rownames `V'

// explained variance by (i1,i2)

	scalar `w2' = 0
	forvalues i = 1 / `=colsof(`W')' {
		scalar `w2' = `w2' + `W'[1,`i']^2
	}
	scalar `rho1' = `W'[1,`i1']^2 / `w2'
	scalar `rho2' = `W'[1,`i2']^2 / `w2'
	scalar `rho'  = `rho1' + `rho2'

// Deal with xnegate, ynegate, stretch, and multanddiv

	if "`xnegate'" != "" {
		mat `U' = (-`U'[1...,1]) , (`U'[1...,2])
		mat `V' = (-`V'[1...,1]) , (`V'[1...,2])
	}

	if "`ynegate'" != "" {
		mat `U' = (`U'[1...,1]) , (-`U'[1...,2])
		mat `V' = (`V'[1...,1]) , (-`V'[1...,2])
	}

	if "`multanddiv'" != "" {
		// multiply U by value   and    divide V by value
		mat `U' = `U' * `multanddiv'
		mat `V' = `V' / `multanddiv'	// previously verified not 0
	}

	if "`stretch'" != "" {
		// multiply V by stretch
		mat `Vhold' = `V'	// ... but hold copy of V for later
		mat `V' = `V' * `stretch'
	}
	
// override the rowlabel if specified
	if "`rlabvar'" != "" {
		capture confirm string var `rlabvar'
		if _rc == 0 {
			forvalues i = 1 / `c(N)' {
				if `touse'[`i'] {
					local rr = `rlabvar'[`i']
					local rr = subinstr("`rr'",".","",.)
					local rr = subinstr("`rr'"," ","_",.)
					local rownames `rownames' `rr'
				}
			}
			matrix rownames `U' = `rownames'
		}
	}


// turn biplot info into variables for plotting
//
// (Note: dim1/dim2 can be filled with faster code (generate) 
//  but I deem this code clearer)

if "`graph'" != "nograph" {

	local nrc = `nr'+`nc'
	if `nrc' > _N {
		preserve
		quietly set obs `nrc'
	}

	tempvar dim1 dim2 labvar

	quietly gen `dim1' = .
	quietly gen `dim2' = .
	quietly gen `labvar' = ""

	forvalues i = 1 / `nr' {
		gettoken lab Unames : Unames

		quietly replace `dim1'   = `U'[`i',1] in `i'
		quietly replace `dim2'   = `U'[`i',2] in `i'
		quietly replace `labvar' = `"`lab'"'  in `i'
	}

	forvalues j = 1 / `nc' {
		local jnr = `nr' + `j'
		gettoken lab Vnames : Vnames

		quietly replace `dim1'   = `V'[`j',1] in `jnr'
		quietly replace `dim2'   = `V'[`j',2] in `jnr'
		quietly replace `labvar' = `"`lab'"'  in `jnr'
	}

// output and plot
			
	summarize `dim1', meanonly
	local dim1_max = `r(max)'
	local dim1_min = `r(min)'

	summarize `dim2', meanonly
	local dim2_max = `r(max)'
	local dim2_min = `r(min)'

	.atk.getAspectAdjustedScales , xmin(`dim1_min') 	///
		xmax(`dim1_max') ymin(`dim2_min') ymax(`dim2_max')
		
	local aspect aspectratio(`s(aspectratio)'`placement')
		
	
	tempvar pos
	_plotpos `dim1' `dim2' , gen(`pos') arrows

	if "`nolabel_rows'" != "nolabel" {
		local rowoptions mlabel(`labvar') `rowoptions'
	}

	if "`nolabel_cols'" != "nolabel" {
		local coloptions mlabel(`labvar') `coloptions'
	}

	tempvar origin
	generate `origin' = 0
	
	if "`colposition_flag'" == "" {
		local col_mlabvpos "mlabvpos(`pos')"
	}
	if "`rowposition_flag'" == "" {
		local row_mlabvpos "mlabvpos(`pos')"
	}	

	if "`separate'" == "" {
		// row and column categories are overlaid

		twoway	pcarrow `origin' `origin' `dim2' `dim1' 	///
			in `=`nr'+1' / `nrc' ,				///
			headlabel `col_mlabvpos' `coloptions'		///
		    || 							///
			scatter `dim2' `dim1' in 1 / `nr' ,		///
			`row_mlabvpos' pstyle(p1) `rowoptions'		///
		    || ,						///
		    	`s(scales)' `aspect'				///
			xtitle(Dimension 1) ytitle(Dimension 2)		///
			legend(label(1 `cname') label(2 `rname') span) 	///
			title("Biplot")  `goptions'
	}
	else {

		// row and column categories in separate plots	
		tempname rplot cplot
		local title

		scatter `dim2' `dim1' in 1/`nr' , 			///
			`rowoptions'					///
		    || ,						///
		    	`s(scales)' `aspect'				///
			name(`rplot') nodraw				///
			xtitle(Dimension 1) ytitle(Dimension 2)		///
			title(`"`rname'"')  `goptions'

		twoway pcarrow `origin' `origin' `dim2' `dim1' 		///
			in `=`nr'+1'/`nrc' ,				///
			mlabvpos(`pos') headlabel `coloptions'		///
		    || ,						///
			name(`cplot') nodraw 				///
			`s(scales)' `aspect'				///
			xtitle(Dimension 1) ytitle(Dimension 2) 	///
			title(`"`cname'"') `goptions'

		graph combine `rplot' `cplot' , 			///
			ycommon xcommon title("Two panel biplot")	///
			`gcopts'
	}

} // end of if nograph
	
// header 

	display _n as txt `"Biplot of {res:`nr'} `=lower("`rname'")' and {res:`nc'} `=lower("`cname'")'"'
	display _n as txt  "    Explained variance by component 1 " as res %7.4f `rho1'
	display    as txt  "    Explained variance by component 2 " as res %7.4f `rho2'
	display    as txt  "             Total explained variance " as res %7.4f `rho'
	
// table output

	if "`table'" != "" { 
		local len1 = length("`rname'")
		local len2 = length("`cname'")
		local length = max(`len1', `len2', 12)
		
		display _n as txt "Biplot coordinates"
		
		matlist `U' , rowtitle(`rname') underscore 	/// 
			format(%8.4f) left(4) border(bot) twidth(`length') 
		   
		matlist `V' , rowtitle(`cname')			/// 
			format(%8.4f) left(4) border(bot) twidth(`length') 
		display
	}	

// save values

	return matrix U = `U'
	if "`stretch'" != "" {
		return matrix V = `Vhold'
		return matrix Vstretch = `V'
	}
	else {
		return matrix V = `V'
	}

	return scalar alpha = `alpha'

	return scalar rho   = `rho'
	return scalar rho1  = `rho1'
	return scalar rho2  = `rho2'
end

exit