In der Arbeit [HH] wurden keine Hypothesen vor der Durchführung des Versuchs festgelegt. Die Autoren leiteten vielmehr aus der Analyse der Daten erst ihre Aussagen ab. Dagegen ist im Prinzip nichts einzuwenden - nur handelt es sich dann nicht um den Nachweis eines Effekts, sondern um die Vorermittlung eines Zusammenhangs, dessen Nachweis durch eine Wiederholungsstudie möglicherweise interessant sein könnte.
Als ein Maß, wie lohnend derartige Replikationsstudien sein könnten, kann der Vergleich mit dem Zufall dienen. Ist das Ergebnis "signifikant", d.h.: Heben sich die Aussagen über den Zusammenhang mit Mondphasen hinreichend deutlich von zufällig gezogenen "Beratungsterminen" ab, die sich über den Beobachtungszeitraum ungefähr gleich verteilen wie der Originaldatensatz?
Diesen Fragen ging Volker Guiard in seiner Arbeit [VG] nach. Der vorliegende Artikel versteht sich als eine Ergänzung zu dieser Arbeit unter Verwendung einer anderen Methode.
Hirschmann und Hirschmann schreiben: "Der Kurvenverlauf ist eindeutig. Die Abhängigkeit von den Mondphasen ist gegeben. Das Maximum liegt hier etwa 8 Tage vor Vollmond, das Minimum etwa sieben Tage nach Vollmond." Hieraus könnte man zwei zu prüfende Hypothesen entnehmen:
Wie liesse sich messen, ob eine Zufallsprobe einen ebenso deutlichen oder gar noch deutlicheren Zusammenhang aufweist? Antwort:
| 1967 | 1968 | 1969 | 1970 | 1971 | 1972 | 1973 | 1974 | 1975 | 1976 |
| 7 | 9 | 24 | 49 | 13 | 24 | 62 | 93 | 131 | 81 |
| 1977 | 1978 | 1979 | 1980 | 1981 | 1982 | 1983 | 1984 | 1985 | 1986 |
| 138 | 62 | 39 | 111 | 219 | 140 | 67 | 89 | 63 | 24 |
| 1987 | 1988 | 1989 | 1990 | 1991 | 1992 | 1993 | 1994 | 1995 | 1996 |
| 8 | 33 | 25 | 26 | 18 | 14 | 28 | 79 | 15 | 35 |
| 1997 | 1998 | ||||||||
| 10 | 30 |
| So | Mo | Di | Mi | Do | Fr | So |
| 382 | 281 | 169 | 223 | 198 | 199 | 314 |
| Juli | August | September | Oktober |
| 230 | 417 | 749 | 370 |
Man kann (Pseudo-)Zufallsdaten mit diesen Häufigkeiten wie folgt erzeugen:
Man erzeugt nicht N=1766, sondern N=7*382 zufällig verteilte Daten. Die Liste wird in sieben Unterlisten nach den Wochentagen zerlegt. Nun werden aus der Liste der Sonntage Null Tage herausgenommen, aus der Liste der Montage 382-281 = 101 Elemente, etc. Welche Elemente herauszunehmen sind, wird wiederum durch Zufallsgeneratoren ermittelt. Im Ergebnis erhält man 1766 zufällige Daten, die ungefähr so wie die Vorgabe auf die Wochentage verteilt sind.
Um auch die Verteilung auf die Monate zu erreichen, startet man mit einer noch größeren Datenmenge, nämlich N = 7*382*( 4*749 / 1766 ) Daten. Die Daten werden in vier Listen nach den Monaten zerlegt. Aus den auf den Juli entfallenden Daten entfernt man nach dem Zufallsprinzip 7*382*(749-230)/1766 Daten, usf. für August und Oktober. Der verbleibende Satz von 7*382 Daten wie wird oben beschrieben hinsichtlich der Wochentage bearbeitet. Man verbleibt mit 1766 zufälligen Daten, die wie vorgegeben auf Wochentage und Monate verteilt sind.
Schließlich kann man auch der Ungleichverteilung der Protokolle auf die Jahre gerecht werden, indem man als Ausgangsmenge für die beschriebenen Ausdünnungen im Jahr x eine Anzahl von N * r(x) Daten erzeugt, wobei r(x) die relative Häufigkeit von Protokollen im Jahre x bedeutet. Da man nur ganzzahlige Datenmengen erzeugen kann, entstehen hierbei kleine Fehler in der Gesamtzahl, die am Ende wieder durch Hinzufügen weiterer Zufallsdaten zu korrigieren sind.
Die Winkelabstände zwischen Sonne und Mond werden, beginnend mit dem Neumond, in 30 "Phasentage" a 12 Grad geteilt. Für jedes Sample wird die Verteilung der Daten auf die 30 Phasentage ermittelt. An die so erhaltenen Anzahlen werden Korrekturfaktoren angebracht, die aus der Ungleichverteilung der Phasentage über den Untersuchungszeitraum resultieren. Die gleichen Faktoren werden auch an die Originalprobe von Hirschmann und Hirschmann angebracht. Statt der Messgröße "Anzahl von Beratungen pro Phasentag" werden wir nach Anbringung dieses Normierungsfaktors von der "normierten Anzahl von Beratungen je Phasentag" sprechen.
Diese Daten werden nun in ein Mathematica Notebook eingelesen. Es wird jeweils ein Least Squares Fit mit einem sinusförmig von den Mondphasen abhängenden Lauf geschätzt und die resultierende mittlere quadratische Abweichung berechnet. Außerdem wird abgezählt, wieviele Daten auf zunehmenden Mond, also auf die Phasentage 0 bis 14 fallen. Schließlich wird abgezählt, in wievielen Samples eine bessere quadratische Abweichung als bei Hirschmann und Hirschmann erzielt wird, und in wie vielen Samples mindestens ebenso viele Daten auf den zunehmendem Mond fielen wie bei Hirschmann und Hirschmann Beratungen.
|
Die links stehende Tabelle enthält die relevanten
Daten, die sich aus der Erhebung von Hirschmann und Hirschmann
ergeben: Die normierte Anzahl von Beratungen je Phasentag.
Diese gilt es zu analysieren. Die Zahlen haben sich - im
Vergleich zu den von [HH] verwendeten - durch die
Normierung an manchen Stellen deutlich verändert. Der Höhepunkt
um den 8. Phasentag verliert durch die Normierung deutlich an Gewicht,
während das "Zwischenhoch" zum Vollmond nun den Höchstwert
annimmt. Diese Verschiebungen ändern jedoch nichts an der
Gestalt der Glättungskurve, die einen Berg in der ersten
und ein Tal in der zweiten Hälfte des Mondlaufs ergibt. Dieser
Trend ist bereits deutlich bei der Zwei-Tage-Glättung ersichtlich,
und nicht bloß bei der von [HH] herangezogenen Sieben-Tage-Glättung.
In den folgenden Abbildungen habe ich die unnormierten Daten nach [HH] zusammen mit ihrer Zwei-Tage-Glättung dargestellt (in der also jeder Funktionswert der Mittelwert aus fünf Säulenhöhen des Balkendiagramms ist). Die normierten Daten wurden zusammen mit ihrer Regressionskurve nach dem Modell f(x) = a sin( Pi x/15 + b) + c
|
| Number of Samples | 1000 |
| Hirschmann/Hirschmann Residual Sum of Squares (RSSHH) | 2431.22 |
| Number of Samples with RSS < RSSHH | 918 |
| Consultations at Moon Wax for Hirschmann/Hirschmann | 933.848 |
| Number of samples with WaxSum < WaxSumHH | 983 |
Hieraus ist zu entnehmen:
Anhand dieses Ergebnisses bin ich geneigt, den sinusförmigen Zusammenhang für ein Zufallsprodukt zu halten, während die Aussage, daß mehr Beratungen auf zunehmenden Mond fallen, erfolgversprechender aussieht. Es könnte möglicherweise lohnend sein, eine Überprüfung dieser zweiten Aussage mit unabhängigem Datenmaterial vozunehmen.
Referenz
[HH] Ursula und Fritz Hirschmann: Unser Mond und sein
Einfluß auf das Pilzwachstum - Lassen sich Zusammenhänge nachweisen?,
Natur und Mensch, Jahresmitteilungen der Naturhistorischen Gesellschaft
Nürnberg e.V., 1999, S. 79-82.
[VG] Volker Guiard: Auswertung einer Datenerhebung
bezüglich des Mondeinflusses auf das Pilzwachstum,
zu erscheinen in Nr. 3/2002 der
Zeitschrift für Anomalistik
(* Load random samples *)
pd = ReadList["Sample.txt", Table[Number,{i,1,30}]];
(* Load Hirschmann/Hirschmann dataset *)
Print["Number of Samples : ", Length[pd]];
hh = ReadList["B_PT.txt", Number];
(* Compute Residual Sum of Squares for HH *)
<< Statistics`NonlinearFit`
f = NonlinearFit[Table[{i,hh[[i]]},{i,1,30}],
a Sin[ Pi x / 15 + b ] + c,
{x},
{a,b,c}];
RSSHH = Sum[(hh[[i]]-f/.x->i)^2, {i,1,30}];
Print["Hirschmann/Hirschmann RSS : ", RSSHH];
(* Now compute the RSS for the samples *)
RSS = Table[0,{i,1,Length[pd]}];
For[j=1, j <= Length[pd], j++,
f = NonlinearFit[Table[{i,pd[[j]][[i]]},{i,1,30}],
a Sin[ Pi x / 15 + b ] + c,
{x},
{a,b,c}];
RSS[[j]] = Sum[ (pd[[j]][[i]]-f/.x->i)^2, {i,1,30}]
];
RSS = Sort[RSS];
For[ i=1, ( RSS[[i]] < RSSHH ) && ( i < Length[pd] ), i++];
Print["Number of Samples with RSS < RSSHH : ", i];
(* The waxing/waning moon hypothesis *)
WaxSumHH = Sum[hh[[j]],{j,1,15}];
Print["Consultations at Moon Wax for Hirschmann/Hirschmann",
"(WaxSumHH) : ", WaxSumHH];
WaxSumSample = Sort[Table[Sum[pd[[i]][[j]],{j,1,15}],{i,1,Length[pd]}]];
For[ i=1, ( WaxSumSample[[i]] < WaxSumHH ) && ( i < Length[pd] ), i++];
Print["Number of samples with WaxSum < WaxSumHH : ",i];
|
//---------------------------------------------------------------------------
// Sample generator for the Moon-Mushroom series
//---------------------------------------------------------------------------
// Generate a series of N=1766 dates in the years 1967-1998
// within the months July-October,
// and with prescribed approximated frequencies of
// - number of dates per day of the week
// - number of dates per month
// - number of dates per year
//---------------------------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#pragma hdrstop
#define s1(dummy) sin((dummy)*1.745329251994329576923691e-2)
#define c1(dummy) cos((dummy)*1.745329251994329576923691e-2)
#define fract(dummy) (dummy - floor(dummy))
//---------------------------------------------------------------------------
#define Pi2 6.28318530717958647692529e0
#define Conv 1.745329251994329576923691e-2
//---------------------------------------------------------------------------
// Ecliptical longitude of the sun, trig series expansion due to S. Newcomb
//---------------------------------------------------------------------------
void newcomb (ti,x)
double ti,*x;
{
double t,t1,dlp,l0,l01,g,g2,g4,g5,g6,d,a_,dl,dl2,dl4,dl5,dl6;
double dlm,l,lx;
t1 = ti - (double)2415020.;
t = t1/(double)36525.;
// --- long periodic disturbation
dlp = ((double)1.882e0-(double).016*t)*
s1((double)57.24 + (double)150.27*t)
+ (double) 6.4 *s1(231.19e0 + 20.20*t)
+ (double)0.266 *s1( 31.8e0 + 119.00*t);
// --- mean longitude sun
l0 = 279.6966778e0 + 360.*fract(t1/(double)365.25);
l01 = 2768.13*t + 1.089*t*t + 0.202*s1(315.6e0+893.3*t)+dlp;
l0 = l0 + l01/3600.e0;
// --- mean anomaly sun
g = 358.475833e0
+ 360.e0*fract(t1*(double).0027104722792607802874743)
+ t1*(double).009828884325804243668
+ ( 179.1*t - .54*t*t + dlp )/3600.e0;
// --- mean anomaly venus,mars,jupiter,saturn
g2 = 212.45e0 + 360.*fract(t1*(double).004435318275154004107)
+ t1*(double).0054070636650308;
g4 = 319.58e0 + t1*(double).524024010951403;
g5 = 225.28e0 + t1*(double).0830823545516769336
+ (double).3611111111111111111* s1(133.775e0 + 39.804*t);
g6 = 175.6e0 + t1*(double).03345089664613278576;
// --- moon data
d = 350.737486e0 + t1*(double).00840832900752909
+ 360.*fract(t1*(double).033839835728952772);
a_ = 296.104608e0 + t1*(double).005444191868583162
+ 360.*fract(t1*(double).0362765229295003422);
// --- terms of two body motion
dl = (6910.057e0 - 17.240e0*t)*s1( g )
+ ( 72.338e0 - .361e0*t)*s1(2.e0*g)
+ 1.054e0*t *s1(3.e0*g );
// --- perturbations in longitude
dl2 = 4.838*c1(299.102e0 + g2 - g)
+ 0.116*c1(148.900e0 + 2.*g2 - g)
+ 5.526*c1(148.313e0 + 2.*g2 - 2.*g)
+ 2.497*c1(315.943e0 + 2.*g2 - 3.*g)
+ 0.666*c1(177.710e0 + 3.*g2 - 3.*g)
+ 1.559*c1(345.243e0 + 3.*g2 - 4.*g)
+ 1.024*c1(318.150e0 + 3.*g2 - 5.*g)
+ 0.210*c1(206.200e0 + 4.*g2 - 4.*g)
+ 0.144*c1(195.400e0 + 4.*g2 - 5.*g)
+ 0.152*c1(343.800e0 + 4.*g2 - 6.*g)
+ 0.123*c1(195.300e0 + 5.*g2 - 7.*g)
+ 0.154*c1(359.600e0 + 5.*g2 - 8.*g);
// ---
dl4 = 0.273*c1(217.700e0 - g4 + g)
+ 2.043*c1(343.888e0 - 2.*g4 + 2.*g)
+ 1.770*c1(200.402e0 - 2.*g4 + g)
+ 0.129*c1(294.200e0 - 3.*g4 + 3.*g)
+ 0.425*c1(338.800e0 - 3.*g4 + 2.*g)
+ 0.500*c1(105.180e0 - 4.*g4 + 3.*g)
+ 0.585*c1(334.060e0 - 4.*g4 + 2.*g)
+ 0.204*c1(100.800e0 - 5.*g4 + 3.*g)
+ 0.154*c1(227.400e0 - 6.*g4 + 4.*g)
+ 0.101*c1( 96.300e0 - 6.*g4 + 3.*g)
+ 0.106*c1(222.700e0 - 7.*g4 + 4.*g);
// ---
dl5 = 0.163*c1(198.600e0 - g5 + 2.*g)
+ 7.208*c1(179.532e0 - g5 + g)
+ 2.600*c1(263.217e0 - g5 )
+ 2.731*c1( 87.145e0 - 2.*g5 + 2.*g)
+ 1.610*c1(109.493e0 - 2.*g5 + g)
+ 0.164*c1(170.500e0 - 3.*g5 + 3.*g)
+ 0.556*c1( 82.650e0 - 3.*g5 + 2.*g)
+ 0.210*c1( 98.500e0 - 3.*g5 + g);
// ---
dl6 = 0.419*c1(100.580e0 - g6 + g)
+ 0.320*c1(269.460e0 - g6 )
+ 0.108*c1(290.600e0 - 2.*g6 + 2.*g)
+ 0.112*c1(293.600e0 - 2.*g6 + g);
// ---
dlm = 6.454*s1(d)
+ 0.177*s1(d+a_)
- 0.424*s1(d-a_)
+ 0.172*s1(d-g);
// --- now sum up to true longitude
l = dl + dl2 + dl4 + dl5 + dl6 + dlm ;
lx = fmod(l/.36e4 + l0,.36e3);
*x = lx;
}
//---------------------------------------------------------------------------
// lunar ephemerides computed to within 1' (brown/i.l.e.)
//---------------------------------------------------------------------------
void brown ( xjd,l)
double xjd,*l;
{
double w,xnut_,day,xl,yl,xp,yp,xn,al,bl,ap,bp,an,
u,v,f,d,dl;
double floor(),sin(),cos(),fmod();
w=(xjd-2415020.e0)/36525.e0;
xnut_ = -17.2327/1296000.*sin((259.18-1934.142*w)*Conv);
day=xjd-2415020.e0;
xl=.751206e0+day*.0366011014634e0;
yl=.776935e0+day*.0027379092649e0;
xp=.928693e0+day*.0003094557786e0;
yp=.781169e0+day*.0000001307457e0;
xn=.719954e0-day*.0001470942283e0;
al=(xl-floor(xl))*Pi2;
bl=(yl-floor(yl))*Pi2;
ap=(xp-floor(xp))*Pi2;
bp=(yp-floor(yp))*Pi2;
an=(xn-floor(xn))*Pi2;
u=al-ap;
v=bl-bp;
f=al-an;
d=al-bl;
dl=22639*sin(u)-4586*sin(u-d-d)+2370*sin(d+d)+769*sin(u+u)
-668*sin(v)-412*sin(f+f)-212*sin(u+u-d-d)-206*sin(u+v-d-d)
+192*sin(u+d+d)-165*sin(v-d-d)+148*sin(u-v)-125*sin(d)
-110*sin(u+v)-55*sin(f+f-d-d)-45*sin(u+f+f)+40*sin(u-f-f)
-38*sin(u-4*d)+36*sin(3*u)-31*sin(u+u-4*d)+28*sin(u-v-d-d)
-24*sin(v+d+d)+19*sin(u-d)+18*sin(v+d)+15*sin(u-v+d+d)
+14*sin(4*d)+14*sin(u+u+d+d)-13*sin(3*u-d-d);
xl += dl/1296000.+xnut_;
*l =(xl-floor(xl))*360.e0;
}
//---------------------------------------------------------------------------
// Calendar conversion civil -> julian
//---------------------------------------------------------------------------
void date1 (int year,
int month,
int day,
double et,
double *xjd)
{
double y,m,b,floor();
// --- calendar date --> julian date
y = (double) year;
m = (double) month;
if (m <= 2.) {
m = m + 12.;
y = y - 1.;
}
b = floor(y/400.)- floor(y/100.);
*xjd = floor(365.25*y) + floor(30.6001*(m+1)) + b +
1720996.5e0 + day + et/24.;
return;
}
//---------------------------------------------------------------------------
// Calendar conversion julian -> civil
//---------------------------------------------------------------------------
void date2 (year,month,day,et,xjd)
int *year,*month,*day;
double *et,xjd;
{
double a,b,c,d,e,f,floor();
// --- julian date --> calendar date
a = floor(xjd + .5);
if (a < 2299161.)
c = a + 1524.;
else {
b = floor((a-1867216.25)/36524.25);
c = a + b - floor(b/4.) + 1525.;
}
d = floor((c-122.1)/365.25);
e = floor(365.25*d);
f = floor((c-e)/30.6001);
*day = c - e - floor(30.6001*f);
*et = (xjd + 0.5 - floor(xjd+0.5))*24.;
*month = f - 1 - (int)(f/14)*12;
*year = d - 4715 - (int)((7+*month)/10);
return;
}
//---------------------------------------------------------------------------
// Compute the moon's phase value
//---------------------------------------------------------------------------
int phase(double juldat) {
double ls, lm, v;
newcomb( juldat, &ls);
brown( juldat, &lm );
v = (lm - ls)/360.;
return (int) floor((v - floor(v))*30.);
}
//---------------------------------------------------------------------------
// Get random Julian date
//---------------------------------------------------------------------------
double randomdate() {
int y,m,d;
static double et = 12.e0, xjd;
// Convert random number into calendar date
d = rand() % 3936;
y = div(d,123).quot + 1967;
d = div(d,123).rem + 1;
if (d < 32) m = 7;
else if (d < 62) { m = 8; d -= 31; }
else if (d < 92) { m = 9; d -= 62; }
else { m = 10; d -= 92; }
// Convert to julian day number
date1(y,m,d,et,&xjd);
return xjd;
}
//---------------------------------------------------------------------------
// Get random Julian date for a single year
//---------------------------------------------------------------------------
double randomdate_year(int year) {
int m,d;
static double et = 12.e0, xjd;
// Convert random number into calendar date
d = rand() % 123 + 1;
if (d < 32) m = 7;
else if (d < 62) { m = 8; d -= 31; }
else if (d < 92) { m = 9; d -= 62; }
else { m = 10; d -= 92; }
// Convert to julian day number
date1(year,m,d,et,&xjd);
return xjd;
}
//---------------------------------------------------------------------------
// Leave the program due to an error
//---------------------------------------------------------------------------
void drop_out() {
fgetc(stdin);
exit(EXIT_FAILURE);
}
#pragma hdrstop
//---------------------------------------------------------------------------
// Main program
//---------------------------------------------------------------------------
#pragma argsused
int main(int argc, char* argv[])
{
int fatal = 0;
FILE *out, *factors;
double xjd0, et, cdpd[30], factor;
long int n, n1;
int i,j,k,i0,
n0 = 0,
iter,
d,m,y,
weekday_fr[7] = { 281, 169, 223, 198, 199, 314, 382 },
month_fr[4] = { 230, 417, 749, 370 },
year_fr[32] = { 7,9,24,49,13,24,62,93,131,81,138,62,39,
111,219,140,67,89,63,24,8,33,25,26,18,
14,28,79,15,35,10,30 },
max_week=0,
max_month=0,
data[5000],
subset[7][2000],
index[7],
phaseday[30]
;
// Initializations
randomize();
// Load normalization factors from disk
factors = fopen("C:\\MondPilz\\CD_PD.txt","r");
i = 0;
while (fscanf(factors, "%lf", &(cdpd[i++])) != EOF) {}
fclose(factors);
for (i=0;i<7;i++) {
if (weekday_fr[i] > max_week) max_week=weekday_fr[i];
n0 += weekday_fr[i];
}
for (i=0;i<4;i++) if (month_fr[i] > max_month) max_month=month_fr[i];
// Start date for all mushroom series
date1( 1967, 7, 0, 12.e0, &xjd0);
out = fopen("C:\\MondPilz\\Sample.txt","w");
factor = 7*max_week*4*max_month/1766./1766.;
// Main Loop: How many samples are to be created
for (iter = 0; iter < 1000; iter++ ) {
// Random date generation, subtract from 0.7.1967
n = 0;
for (i = 0; i<32; i++)
{
for (j=0; j <= factor*year_fr[i]; j++)
data[n++] = (int) floor(randomdate_year(i+1967) - xjd0);
}
// First step - filter out the months
for (i=0; i<4; i++) index[i] = 0;
// Form 4 subsets for the months
for (i=0; i<n; i++) {
date2( &y, &m, &d, &et, xjd0 + data[i] );
if ((m >= 7) && (m<11)) {
subset[m-7][index[m-7]++] = data[i];
if (index[m-7] == 1999) {
printf("Datata overflow, month %d\n",m);
fatal = 1;
break;
}
}
else {
fatal = 1;
printf("Fatal program error! Month %d\n",m);
break;
}
}
if (fatal) drop_out();
// Now we have the 4 subsets - thin them out
for ( i = 0; i < 4; i++) {
n1 = (max_month - month_fr[i])*max_week*7/1766;
i0 = index[i];
if (i0 == 0) {
fatal = 1;
printf("Program error: index[%d] = 0\n",i);
drop_out();
}
for (j = 0; j < n1; j++ ) {
// Use negative number to mark them, skip numbers already cancelled
while ( subset[i][k = (rand() % index[i])] < 0 ) {}
subset[i][k] = -1;
}
}
// Now rejoin the subsets
n = 0;
for ( i = 0; i < 4; i++) {
for (j = 0; j < index[i]; j++)
if (subset[i][j] >= 0) data[n++] = subset[i][j];
}
// Second step - filter out the days of the week
for (i=0; i<7; i++) index[i] = 0;
for (i=0; i<n; i++) {
d = div((int)floor(data[i]+xjd0),7).rem;
subset[d][index[d]++] = data[i];
if (index[d] == 1999) {
printf("Data overflow, month %d\n",m);
fatal = 1;
break;
}
}
if (fatal) drop_out();
// Now we have the 7 subsets - thin them out
for ( i = 0; i < 7; i++) {
n1 = max_week - weekday_fr[i];
i0 = index[i];
if (i0 == 0) {
fatal = 1;
printf("Program error: index[%d] = 0\n",i);
drop_out();
}
for (j = 0; j < n1; j++ ) {
// Use negative number to mark them, skip numbers already cancelled
while ( subset[i][k = (rand() % index[i])] < 0 ) {}
subset[i][k] = -1;
}
}
// Now rejoin the subsets
n = 0;
for ( i = 0; i < 7; i++) {
for (j = 0; j < index[i]; j++)
if (subset[i][j] >= 0) data[n++] = subset[i][j];
}
// Add some more random dates to achieve prescribed size
while (n < 1766) {
data[n++] = (int) floor(randomdate(i+1967) - xjd0);
}
// ... or thin the data out if necessary
if (n > 1766) {
for (i=0; i < n - 1766; i++) {
while(data[j=rand()%(n-1766)] == -1){}
data[j] = -1;
}
j = 0;
for (i = 0; i < n; i++) {
if (data[i] < 0) continue;
if (j < i) data[j++] = data[i];
}
if (j != 1766) {
printf("Program error - j = %d\n",j);
drop_out();
}
else n = 1766;
}
// Now: compute phaseday[i], i=1,..,30
for (i=0; i<30; i++) phaseday[i] = 0;
for (i=0; i<n; i++)
phaseday[phase(xjd0+data[i])]++;
for (i = 0; i < 30; i++) fprintf(out,"%9.4lf",phaseday[i]*cdpd[i]);
fprintf(out,"\n");
} // End of main loop
fclose(out);
return 0;
}
//---------------------------------------------------------------------------
|
Nussloch, im Februar 2002
Rüdiger Plantiko
