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
|
include <mach.h>
include <ctype.h>
include <error.h>
include <gset.h>
include <pkg/rg.h>
include "pdm.h"
define NUMTRIES 100
# PDM_SIGNIF -- Calculate the significance of the theta statistic for
# a certain period. Use the "Method of Randomization".
double procedure pdm_signif (pdmp, period)
pointer pdmp # pointer to PDM data structure
double period # period at which to find significance
int lesscount, i, npt
double otheta, theta, pdm_theta(), pdm_thetaran()
long seed
pointer rg, pt, inuse, oinuse, rg_xrangesd(), sp
errchk pdm_statistics, pdm_theta(), pdm_ranperm, pdm_thetaran()
begin
# Do NUMTRIES random permutations on the data, calculate the theta
# statistic on the scrambled data for this period. Return the
# fraction of these permutaions which yield thetas less than
# the unmixed data.
call smark (sp)
npt = PDM_NPT(pdmp)
# Make sure the statistics are up to date.
call pdm_statistics (pdmp)
# Allocate a temporary array for the scrambled data, allocate a
# temporary copy of the inuse array and copy the real inuse array
# into it.
call salloc (pt, npt, TY_DOUBLE)
call salloc (inuse, npt, TY_INT)
call salloc (oinuse, npt, TY_INT)
call amovi (PDM_INUSE(pdmp,1), Memi[inuse], npt)
lesscount = 0
# Calculate the ranges information from the sample string.
rg = rg_xrangesd (PDM_SAMPLE(pdmp), PDM_X(pdmp,1), npt)
otheta = pdm_theta (pdmp, rg, period)
seed = 1.0
do i = 1, NUMTRIES {
call pdm_ranperm (PDM_DY(pdmp,1), Memi[inuse], Memd[pt],
Memi[oinuse], npt, seed)
theta = pdm_thetaran (pdmp, pt, oinuse, rg, period)
if (theta < otheta)
lesscount = lesscount + 1
}
call sfree (sp)
return (1.0d+0 - (double(lesscount)/double(NUMTRIES)))
end
|
