File restructure #1
This commit is contained in:
@ -1,198 +0,0 @@
|
||||
*! version 1.0.1 TJS 9jun2000
|
||||
program define nctprob, rclass
|
||||
version 6.0
|
||||
args t delta df extra
|
||||
|
||||
if "`df'" == "" | "`extra'" != "" {
|
||||
di in gr "Syntax for " in wh "nctprob" _c
|
||||
di in gr ", the cumulative non-central t"
|
||||
di in gr "distribution from negative infinity to t', is: " _n
|
||||
di in wh " nctprob " in gr "t' delta df" _n
|
||||
di in gr " where " in wh "t' " in gr "is the observed t"
|
||||
di in wh " delta " in gr "is the noncentrality parameter"
|
||||
di in wh " df " in gr "is the degrees of freedom" _n
|
||||
di in wh " nctprob " in gr "computes " in wh "p" _c
|
||||
di in gr " such that P(t<=" in wh "t'" in gr "| " in wh "delta" _c
|
||||
di in gr ", " in wh "df" in gr ") = " in wh "p"
|
||||
di in gr " and returns the value in result " _c
|
||||
di in wh "r(p) " in gr "and global " in wh "S_1" in gr "."
|
||||
global S_1 = .
|
||||
return scalar p = .
|
||||
exit 9
|
||||
}
|
||||
|
||||
if `df' != int(`df') | `df' < 1 {
|
||||
di in re "degrees of freedom must be a positive integer"
|
||||
global S_1 = .
|
||||
return scalar p = .
|
||||
exit 498
|
||||
}
|
||||
|
||||
local even = mod(`df',2) == 0
|
||||
|
||||
/* numerical calculation of C(h,a) requires -preserve-
|
||||
but C(h,a) is only needed for odd df's */
|
||||
if !`even' { preserve }
|
||||
|
||||
tempname A B h a p C M0 M1 M2 ak Mo Me k Mk
|
||||
tempvar x y
|
||||
|
||||
scalar `A' = `t' / sqrt(`df')
|
||||
scalar `B' = `df' / (`df' + (`t')^2)
|
||||
scalar `h' = `delta' * sqrt(`B')
|
||||
|
||||
if `even' { scalar `p' = normprob(-(`delta')) }
|
||||
else {
|
||||
scalar `p' = normprob(-(`delta' * sqrt(`B')))
|
||||
scalar `a' = abs(`A')
|
||||
qui range `x' 0 `a' 1001
|
||||
gen `y' = exp(-((`h')^2 / 2) * (1 + (`x')^2)) / (1 + (`x')^2)
|
||||
qui integ `y' `x'
|
||||
scalar `C' = r(integral)/ (2 * _pi)
|
||||
scalar `p' = `p' + 2 * `C' * sign(`A')
|
||||
}
|
||||
if `df' == 1 {
|
||||
di _n in gr " p =" in ye %10.6f `p'
|
||||
di _n in gr " P(t <= " `t' " | delta = " _c
|
||||
di in gr `delta' ", df = " `df' _c
|
||||
di in gr ") = " in ye `p'
|
||||
global S_1 = `p'
|
||||
return scalar p = `p'
|
||||
exit
|
||||
}
|
||||
|
||||
scalar `M0' = `A' * sqrt(`B') * normd(`delta' * sqrt(`B'))
|
||||
scalar `M0' = `M0' * normprob(`delta' * `A' * sqrt(`B'))
|
||||
if `df' == 2 {
|
||||
scalar `p' = `p' + sqrt(2 * _pi) * `M0'
|
||||
di _n in gr " p =" in ye %10.6f `p'
|
||||
di _n in gr " P(t <= " `t' " | delta = " _c
|
||||
di in gr `delta' ", df = " `df' _c
|
||||
di in gr ") = " in ye `p'
|
||||
global S_1 = `p'
|
||||
return scalar p = `p'
|
||||
exit
|
||||
}
|
||||
|
||||
scalar `M1' = `A' * normd(`delta') / sqrt(2 * _pi)
|
||||
scalar `M1' = `B' * (`delta' * `A' * `M0' + `M1')
|
||||
if `df' == 3 {
|
||||
scalar `p' = `p' + 2 * `M1'
|
||||
di _n in gr " p =" in ye %10.6f `p'
|
||||
di _n in gr " P(t <= " `t' " | delta = " _c
|
||||
di in gr `delta' ", df = " `df' _c
|
||||
di in gr ") = " in ye `p'
|
||||
global S_1 = `p'
|
||||
return scalar p = `p'
|
||||
exit
|
||||
}
|
||||
|
||||
scalar `M2' = `B' * ( `delta' * `A' * `M1' + `M0') / 2
|
||||
if `df' == 4 {
|
||||
scalar `p' = `p' + sqrt(2 * _pi) * (`M0' + `M2')
|
||||
di _n in gr " p =" in ye %10.6f `p'
|
||||
di _n in gr " P(t <= " `t' " | delta = " _c
|
||||
di in gr `delta' ", df = " `df' _c
|
||||
di in gr ") = " in ye `p'
|
||||
global S_1 = `p'
|
||||
return scalar p = `p'
|
||||
exit
|
||||
}
|
||||
|
||||
* calculate Mk's for k = 3 to `df'-2 and sum odds and evens
|
||||
scalar `ak' = 1
|
||||
scalar `Mo' = `M1'
|
||||
scalar `Me' = `M0' + `M2'
|
||||
scalar `k' = 3
|
||||
while `k' <= `df' - 2 {
|
||||
scalar `ak' = 1 / ((`k' - 2) * `ak')
|
||||
scalar `Mk' = `ak' * `delta' * `A' * `M2' + `M1'
|
||||
scalar `Mk' = (`k'-1) * `B' * `Mk' / `k'
|
||||
if mod(`k',2) == 0 { scalar `Me' = `Me' + `Mk' }
|
||||
else { scalar `Mo' = `Mo' + `Mk' }
|
||||
scalar `M1' = `M2'
|
||||
scalar `M2' = `Mk'
|
||||
scalar `k' = `k' + 1
|
||||
}
|
||||
if `even' { scalar `p' = `p' + sqrt(2 * _pi) * `Me' }
|
||||
else { scalar `p' = `p' + 2 * `Mo' }
|
||||
|
||||
di _n in gr " p =" in ye %10.6f `p'
|
||||
di _n in gr " P(t <= " `t' " | delta = " _c
|
||||
di in gr `delta' ", df = " `df' _c
|
||||
di in gr ") = " in ye `p'
|
||||
|
||||
global S_1 = `p'
|
||||
return scalar p = `p'
|
||||
|
||||
exit
|
||||
end
|
||||
|
||||
/* ------------------------------------------------------------
|
||||
|
||||
Note: formula implemented above is from
|
||||
D. B. Owen (Technometrics 10(3):445-478, 1968)
|
||||
|
||||
|
||||
Let
|
||||
t be the observed t-value
|
||||
d be the non-centrality parameter
|
||||
v be the degrees of freedom
|
||||
|
||||
G(z) be the cumulative standard Normal distribution
|
||||
G'(z) be the standard Normal density function
|
||||
|
||||
Define
|
||||
|
||||
A = t / sqrt(v)
|
||||
|
||||
B = v / (v + t^2)
|
||||
|
||||
then
|
||||
|
||||
M = 0
|
||||
-1
|
||||
|
||||
M = A * sqrt(B) * G'(d * sqrt(B)) * G(d * A * sqrt(B))
|
||||
0
|
||||
|
||||
M = B * [ d * A * M + A * G'(d) / sqrt(2 * pi)]
|
||||
1 0
|
||||
|
||||
M = B * [ d * A * M + M ]
|
||||
2 0 1
|
||||
|
||||
and, for k>= 3,
|
||||
|
||||
M = (k - 1) * B * [ a * d * A * M + M ] / k
|
||||
k k k-1 k-2
|
||||
|
||||
where a = 1 / [(k - 2) * a ] and a = 1
|
||||
k k-1 2
|
||||
|
||||
|
||||
Finally, for even df's,
|
||||
|
||||
P{T <= t | d, v} =
|
||||
G(-d) + sqrt(2 * pi) * [M + M + ... + M ]
|
||||
0 2 v-2
|
||||
|
||||
and for odd df's,
|
||||
|
||||
P{T <= t | d, v} =
|
||||
G(-d * sqrt(B)) + 2 * C(d * sqrt(B), A) + 2 * [M + M + ... + M ]
|
||||
1 3 v-2
|
||||
|
||||
|
||||
where
|
||||
abs(a)
|
||||
1 / exp[-h^2 / 2 * (1 + x^2)]
|
||||
C(h,a) = sign(a) -------- | --------------------------- dx
|
||||
2 * pi / 1 + x^2
|
||||
x=0
|
||||
|
||||
|
||||
Note: integral C(h,a) is computed numerically in this program.
|
||||
|
||||
------------------------------------------------------------ */
|
||||
|
||||
Reference in New Issue
Block a user