-
Notifications
You must be signed in to change notification settings - Fork 44
/
Copy pathlreg.gs
151 lines (129 loc) · 3.87 KB
/
lreg.gs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
function lreg( args )
_version = '0.01b1'
rc = gsfallow("on")
if( args = '' )
help()
return
endif
***** Arguement *****
'vx = 'subwrd( args, 1 )
'vy = 'subwrd( args, 2 )
* optional
slope_in = subwrd( args, 3 )
icept_in = subwrd( args, 4 )
cor_in = subwrd( args, 5 )
xmin = qdims( xmin ) ; xmax = qdims( xmax )
ymin = qdims( ymin ) ; ymax = qdims( ymax )
zmin = qdims( zmin ) ; zmax = qdims( zmax )
tmin = qdims( tmin ) ; tmax = qdims( tmax )
emin = qdims( emin ) ; emax = qdims( emax )
str_head = ''
str_tail = ''
if( xmin != xmax )
str_head = 'ave(' % str_head
str_tail = str_tail % ',x='xmin',x='xmax')'
endif
if( ymin != ymax )
str_head = 'ave(' % str_head
str_tail = str_tail % ',y='ymin',y='ymax')'
endif
if( zmin != zmax )
str_head = 'ave(' % str_head
str_tail = str_tail % ',z='zmin',z='zmax')'
endif
if( tmin != tmax )
str_head = 'ave(' % str_head
str_tail = str_tail % ',t='tmin',t='tmax')'
endif
if( emin != emax )
str_head = 'ave(' % str_head
str_tail = str_tail % ',e='emin',e='emax')'
endif
'set x 1'
'set y 1'
'set z 1'
'set t 1'
'set e 1'
* avx, avy: ave(vx), ave(vy)
'tmp = ' % str_head % 'vx' % str_tail
'd tmp'
avx = subwrd( result, 4 )
say 'avx = ' % str_head % 'vx' % str_tail % ' = ' % avx
'tmp = ' % str_head % 'vy' % str_tail
'd tmp'
avy = subwrd( result, 4 )
say 'avy = ' % str_head % 'vy' % str_tail % ' = ' % avy
* avxy: ave(vx*vy)
'tmp = ' % str_head % 'vx*vy' % str_tail
'd tmp'
avxy = subwrd( result, 4 )
say 'avxy = ' % str_head % 'vx*vy' % str_tail % ' = ' % avxy
* avxx: ave(vx*vx)
'tmp = ' % str_head % 'vx*vx' % str_tail
'd tmp'
avxx = subwrd( result, 4 )
say 'avxx = ' % str_head % 'vx*vx' % str_tail % ' = ' % avxx
* sgx, sgy: sigma(vx), sigma(vy)
'tmp = sqrt(' % str_head % '(vx-'avx')*(vx-'avx')' % str_tail % ')'
'd tmp'
sgx = subwrd( result, 4 )
say 'sgx = sqrt(' % str_head % '(vx-avx)*(vx-avx)' % str_tail % ' = ' % sgx
'tmp = sqrt(' % str_head % '(vy-'avy')*(vy-'avy')' % str_tail % ')'
'd tmp'
sgy = subwrd( result, 4 )
say 'sgy = sqrt(' % str_head % '(vy-avy)*(vy-avy)' % str_tail % ' = ' % sgy
* linear regression coefficient
slope = ( avxy - avx * avy ) / ( avxx - avx * avx )
icept = avy - slope * avx
say 'y = ' % slope % ' x + ' % icept
* correlation coefficient
'tmp = ' % str_head % '(vx-'avx')*(vy-'avy')' % str_tail % ' / ('sgx'*'sgy')'
'd tmp'
cor = subwrd( result, 4 )
say 'correlation: ' % cor
if( slope_in != '' )
slope_in' = 'slope
endif
if( icept_in != '' )
icept_in' = 'icept
endif
if( cor_in != '' )
cor_in' = 'cor
endif
ret = slope % ' ' % icept % ' ' % cor
'set x 'xmin' 'xmax
'set y 'ymin' 'ymax
'set z 'zmin' 'zmax
'set t 'tmin' 'tmax
'set e 'emin' 'emax
return ret
*
* help
*
function help()
say ' Name:'
say ' lreg '_version' - linear regression and correlation for any varying dimensions'
say ' '
say ' Usage:'
say ' lreg vx vy [slope [icept [cor]]]'
say ''
say ' vx,vy : input field variables'
say ' slope : output field variable name for slope parameter'
say ' icept : output field variable name for intercept parameter'
say ' cor : output field variable name for correlation coefficient'
say ' Output as "result"'
say '1. slope parameter in linear regression'
say '2. intercept parameter in linear regression'
say '3. correlation coefficient'
say ''
say ' Note:'
say ' This script is under construction and its syntax will perhaps be changed in near future.'
say ' lon/lat is weighted...'
say ''
say ' [arg-name] : specify if needed'
say ' (arg1 | arg2) : arg1 or arg2 must be specified'
say ''
say ' Copyright (C) 2012-2012 Chihiro Kodama'
say ' Distributed under GNU GPL (http://www.gnu.org/licenses/gpl.html)'
say ''
return