Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: st: query regarding thin plate spline method
From
behailu ayele <[email protected]>
To
[email protected]
Subject
Re: st: query regarding thin plate spline method
Date
Tue, 6 Dec 2011 09:27:01 +0000
Dear statalisters,
I am trying to do interpoliations to my climate data using the thin
plate spline method. I searched in statalist and I found an ado file
for that (which I paste below). I have now faced two problems: when I
install THINPLATE and try to recall it in stata, I get error message
saying THINPLATE is not found. The second problem is when I just copy
the ado file and run it straigt in stata, the lines run fine but I get
no output. My programming knowledge in stata is poor so any help on
this is much appreciated.
Also if there are easier ways of doing interpoliations in stata, you
may sugget them to me as well.
Thank you in advance.
Behailu
*! Thin-plate base spline: v.1.1
program define thinplate
version 7
gettoken action 0 : 0
if "`action'"=="generate" { TPGen `0' }
if "`action'"=="cluster" { TPClust `0' }
if "`action'"=="graph" { TPGraph `0' }
end
program define TPGraph
* syntax: depvar x y, cluster(how many) level
*
#delimit ;
syntax varlist (num min=3 max=3) [if] [in], LEVel(numlist)
[CLuster(int 0) Grid(int 20) Repeat(int 1)
Symbol(str) Connect(str) T1title(str) T2title(str) L1title(str) *]
;
#delimit cr
qui {
if "`symbol'"=="" { local symbol T }
tokenize `symbol', parse("[ ]")
local grsym `2'
marksample touse
tokenize `varlist'
local depvar `1'
local x `2'
local y `3'
if `cluster'==0 {
count if `touse'
local `cluster'=int(log(r(N)))
}
local nlevels: word count `level'
if `nlevels'>5 {
di as err "cannot draw more than 4 levels"
exit 198
}
preserve
keep `varlist' `touse' `grsym'
keep if `touse'
* 1. create the grid x grid for the graph
sum `x'
local xmin=r(min)
local xmax=r(max)
local xh=(`xmax'-`xmin')/(`grid'-1)
sum `y'
local ymin=r(min)
local ymax=r(max)
local yh=(`ymax'-`ymin')/(`grid'-1)
count
local nobswas = r(N)
local nobswas1 = `nobswas'+1
local nobswas2 = `nobswas'+2
local nobs = `nobswas' + `grid'
set obs `nobs'
replace `x'=`xmin' in `nobswas1'
replace `x'=`x'[_n-1]+`xh' in `nobswas2'/l
expand `grid' in `nobswas1'/l
sort `x' `touse'
by `x' `touse': replace `y'=`ymin' if _n==1 & `touse'==.
by `x' `touse': replace `y'=`y'[_n-1]+`yh' if _n>1 & `touse'==.
forvalues i=1/`repeat' {
* 2. create the thin plate splines
cap drop _TPS*
* if !_rc { di as error _n "Sorry, your _TPS* variables are
dropped..." }
TPClust _TPS = `x' `y' if `touse', k(`cluster') start(prandom)
* 3. predict the values
reg `depvar' _TPS* if `touse'
tempname fit`i'
predict `fit`i'' if `touse'==.
local fitlist `fitlist' `fit`i''
}
tempvar avfit labfit oldsitey
egen `avfit'=rmean(`fitlist')
tokenize `level'
g str1 `labfit' = " " if `touse'==.
replace `labfit' = "." if `touse'==. & `avfit'<`1'
replace `labfit' = "-" if `touse'==. & `avfit'>=`1'
local legend . < `1' < -
if "`2'"~="" {
replace `labfit' = "+" if `touse'==. & `avfit'>=`2'
local legend `legend' < `2' < +
}
if "`3'"~="" {
replace `labfit' = "o" if `touse'==. & `avfit'>=`3'
local legend `legend' < `3' < o
}
if "`4'"~="" {
replace `labfit' = "X" if `touse'==. & `avfit'>=`4'
local legend `legend' < `4' < X
}
local legend `legend' in terms of `depvar'
local ldepvar : var label `depvar'
if "`ldepvar'"=="" { local ldepvar `depvar' }
g `oldsitey' = `y' if `touse'==1
local laby : var label `y'
graph `y' `oldsitey' `x', title(Fit of `ldepvar' by thin plate
splines) t1(`legend') t2(Triangles: observed sites) l1(`laby')
s([`labfit']`symbol') `options'
* 4. cleanup
drop _TPS*
}
* end of quietly
end
program define TPClust
* syntax:
* TPClust newvariables = x y if in, cluster kmeans options
gettoken prefix 0: 0
gettoken eqs 0: 0
if "`eqs'"~="=" { error 198 }
syntax varlist(num min=2 max=2) [if] [in], K(int) [START(passthru)]
marksample touse
tempvar clust
qui {
if "`start'"=="" { local start start(prandom) }
cluster kmeans `varlist' if `touse', k(`k') gen(`clust') l2 `start'
tokenize `varlist'
forvalues i=1/`k' {
sum `1' if `clust'==`i', meanonly
local x=r(mean)
sum `2' if `clust'==`i', meanonly
local y=r(mean)
local clist `clist' `x' `y'
}
} /* end of quietly */
TPGen `prefix' = `varlist', knots(`clist')
end
program define TPGen
* syntax:
* TPGen new variables = x y, knots(list of centers)
gettoken prefix 0: 0
gettoken eqs 0: 0
if "`eqs'"~="=" { error 198 }
syntax varlist(num min=2 max=2) [if] [in] , Knots(numlist
min=2) [DOUBLE]
marksample touse
tokenize `varlist'
local xvar `1'
local yvar `2'
tokenize `knots'
local nc : word count `knots'
if `nc'-int(`nc'/2)*2 {
di as err "Coordinates of the knots are not valid"
exit 198
}
local nc=`nc'/2
forvalues i=1/`nc' {
local kx = 2*`i'-1
local ky = 2*`i'
local x ``kx''
local y ``ky''
tempname spline`i'
qui g `double'
`spline`i''=(`x'-`xvar')*(`x'-`xvar')+(`y'-`yvar')*(`y'-`yvar')
qui replace
`spline`i''=`spline`i''*log(`spline`i'')/(16*_pi) if `spline`i''>0 &
`spline`i''~=.
}
nobreak {
forvalues i=1/`nc' {
ren `spline`i'' `prefix'`i'
}
qui count if `prefix'1==.
if r(N) { di as text r(N) " missing values generated" }
if `nc'==1 { rename `prefix'1 `prefix' }
}
end
/*
History:
04.03.2001 thinplate generate | cluster | graph, help file
*/
On 12/6/11, behailu ayele <[email protected]> wrote:
> Dear statalisters,
>
> I am trying to do interpoliations to my climate data using the thin
> plate spline method. I searched in statalist and I found an ado file
> for that (which I paste below). I have now faced two problems: when I
> install THINPLATE and try to recall it in stata, I get error message
> saying THINPLATE is not found. The second problem is when I just copy
> the ado file and run it straigt in stata, the lines run fine but I get
> no output. My programming knowledge in stata is poor so any help on
> this is much appreciated.
>
> Also if there are easier ways of doing interpoliations in stata, you
> may sugget them to me as well.
>
> Thank you in advance.
>
> Behailu
> *
> * 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/