dsource.org

Open Source Development for
the D Programming Language

Navigation

← Previous Changeset
Next Changeset →

Changeset 164

Timestamp:

06/19/10 18:23:09 (2 years ago)

Author:

dsimcha

Message:

Add AUROC, misc. 2.047-related updates.

Files:

trunk/alloc.d (modified) (2 diffs)
trunk/base.d (modified) (6 diffs)
trunk/cor.d (modified) (10 diffs)
trunk/distrib.d (modified) (49 diffs)
trunk/infotheory.d (modified) (3 diffs)
trunk/random.d (modified) (22 diffs)
trunk/regress.d (modified) (6 diffs)
trunk/tests.d (modified) (23 diffs)

Legend:

: Unmodified
: Added
: Removed
: Modified
: Copied
: Moved

trunk/alloc.d

r159	r164
1955	1955	* sorted on.*/
1956	1956	T find(U)(U whatToFind) {
1957		T* ptr = enforce( opIn_r!(U)(whatToFind),
	1957	T* ptr = dstatsEnforce( opIn_r!(U)(whatToFind),
1958	1958	"Item not found: " ~ to!string(whatToFind));
1959	1959	return *ptr;
…	…
2083	2083
2084	2084	// TreeRange!(T, mapFun) asRange() {
2085		// enforce(0, "Not implemented yet.");
	2085	// dstatsEnforce(0, "Not implemented yet.");
2086	2086	// }
2087	2087

trunk/base.d

r163	r164
65	65	void main (){}
66	66	}
	67
	68	class DstatsArgumentException : Exception {
	69	this(string msg) {
	70	super(msg);
	71	}
	72	}
	73
	74	T dstatsEnforce(T, string file = __FILE__, int line = __LINE__)
	75	(T value, lazy const(char)[] msg = null) {
	76	if(!value) {
	77	auto exceptMsg = (msg !is null) ? msg.idup : "Invalid argument.";
	78	throw new DstatsArgumentException(file ~ " (" ~ text(line) ~ ") : " ~
	79	exceptMsg);
	80	}
	81
	82	return value;
	83	}
	84
67	85
68	86	/** Tests whether T is an input range whose elements can be implicitly
…	…
110	128	* foreach(elem; T.init) {}.*/
111	129	template IterType(T) {
112		alias ReturnType!(
	130	alias ReturnType!(typeof(
113	131	{
114	132	foreach(elem; T.init) {
…	…
116	134	}
117	135	assert(0);
118		}) IterType;
	136	})) IterType;
119	137	}
120	138
…	…
151	169	}
152	170	return output;
	171	}
	172
	173	/**Given a tuple possibly containing forward ranges, returns a tuple where
	174	* save() has been called on all forward ranges.
	175	*/
	176	Tuple!T saveAll(T...)(T args) {
	177	Tuple!T ret;
	178
	179	foreach(ti, elem; args) {
	180	static if(isForwardRange!(typeof(elem))) {
	181	ret.field[ti] = elem.save;
	182	} else {
	183	ret.field[ti] = elem;
	184	}
	185	}
	186
	187	return ret;
153	188	}
154	189
…	…
670	705	}
671	706
	707	/**Finds the area under the ROC curve (a curve with sensitivity on the Y-axis
	708	* and 1 - specificity on the X-axis). This is a useful metric for
	709	* determining how well a test statistic discriminates between two classes.
	710	* The following assumptions are made in this implementation:
	711	*
	712	* 1. For some cutoff value c and test statistic T, your decision rule is of
	713	* the form "Class A if T > c, Class B if T < c".
	714	*
	715	* 2. In the case of ties, i.e. if class A and class B both have an identical
	716	* value, linear interpolation is used. This is because changing the
	717	* value of c infinitesimally will change both sensitivity and specificity
	718	* in these cases.
	719	*/
	720	double auroc(R1, R2)(R1 classATs, R2 classBTs)
	721	if(isNumeric!(ElementType!R1) && isNumeric!(ElementType!R2)) {
	722	mixin(newFrame);
	723	auto classA = tempdup(classATs);
	724	auto classB = tempdup(classBTs);
	725	qsort(classA);
	726	qsort(classB);
	727
	728	// Start cutoff at -infinity, such that we get everything in class A, i.e.
	729	// perfect specificity, zero sensitivity. We arbitrarily define class B
	730	// as our "positive" class.
	731	double tp = 0, tn = classA.length, fp = 0, fn = classB.length;
	732	double[2] lastPoint = 0;
	733
	734	CommonType!(ElementType!R1, ElementType!R2) currentVal;
	735
	736	ElementType!R1 popA() {
	737	tn--;
	738	fp++;
	739	auto ret = classA.front();
	740	classA.popFront();
	741	return ret;
	742	}
	743
	744	ElementType!R2 popB() {
	745	fn--;
	746	tp++;
	747	auto ret = classB.front();
	748	classB.popFront();
	749	return ret;
	750	}
	751
	752	double area = 0;
	753	while(!classA.empty && !classB.empty) {
	754	if(classA.front() < classB.front()) {
	755	currentVal = popA();
	756	} else {
	757	currentVal = popB();
	758	}
	759
	760	// Handle ties.
	761	while(!classA.empty && classA.front() == currentVal) {
	762	popA();
	763	}
	764
	765	while(!classB.empty && classB.front() == currentVal) {
	766	popB();
	767	}
	768
	769	double[2] curPoint;
	770	curPoint[0] = 1.0 - tn / (fp + tn);
	771	curPoint[1] = tp / (tp + fn);
	772
	773	immutable xDist = curPoint[0] - lastPoint[0];
	774	area += xDist * lastPoint[1]; // Rectangular part.
	775	area += xDist * 0.5 * (curPoint[1] - lastPoint[1]); // Triangular part.
	776	lastPoint[] = curPoint[];
	777	}
	778
	779	if(classA.length > 0 && classB.length == 0) {
	780	// Then we already have a sensitivity of 1, move straight to the right
	781	// to the point (1, 1).
	782
	783	immutable xDist = 1 - lastPoint[0];
	784	area += xDist * lastPoint[1]; // Rectangular part.
	785	area += xDist * 0.5 * (1 - lastPoint[1]); // Triangular part.
	786	}
	787
	788	return area;
	789	}
	790
	791	unittest {
	792	// Values worked out by hand on paper. If you don't believe me, work
	793	// them out yourself.
	794	assert(auroc([4,5,6], [1,2,3]) == 1);
	795	assert(approxEqual(auroc([8,6,7,5,3,0,9], [3,6,2,4,3,6]), 0.6904762));
	796	assert(approxEqual(auroc([2,7,1,8,2,8,1,8], [3,1,4,1,5,9,2,6]), 0.546875));
	797
	798	writeln("Passed auroc unittest.");
	799	}
	800
672	801	///
673	802	T sign(T)(T num) pure nothrow if(is(typeof(num < 0))) {
…	…
881	1010
882	1011	///
883		bool empty() ~~const pure nothrow~~ {
	1012	bool empty() @property {
884	1013	return nPerms == 0;
885	1014	}

trunk/cor.d

r161	r164
31	31	module dstats.cor;
32	32
33		import std.range, std.typecons, std.contracts, std.math, std.traits;
	33	import std.range, std.typecons, std.contracts, std.math, std.traits, std.typetuple;
34	34
35	35	import dstats.sort, dstats.base, dstats.alloc, dstats.regress : invert;
…	…
68	68	// properly aligned, this can result in about 2x speedups compared
69	69	// to simply submitting everything to a single PearsonCor struct.
70		enforce(input1.length == input2.length,
	70	dstatsEnforce(input1.length == input2.length,
71	71	"Ranges must be same length for Pearson correlation.");
72	72
…	…
120	120	}
121	121
122		enforce(input1.empty && input2.empty,
	122	dstatsEnforce(input1.empty && input2.empty,
123	123	"Ranges must be same length for Pearson correlation.");
124	124	}
…	…
336	336	auto ranks1 = spearmanCorRank(input1);
337	337	auto ranks2 = spearmanCorRank(input2);
338		enforce(ranks1.length == ranks2.length,
	338	dstatsEnforce(ranks1.length == ranks2.length,
339	339	"Ranges must be same length for Spearman correlation.");
340	340
…	…
429	429	if(isInputRange!(T) && isInputRange!(U)) {
430	430	static if(isArray!(T) && isArray!(U)) {
431		enforce(input1.length == input2.length,
	431	dstatsEnforce(input1.length == input2.length,
432	432	"Ranges must be same length for Kendall correlation.");
433	433	if(input1.length <= kendallSmallN) {
…	…
441	441	scope(exit) TempAlloc.free;
442	442
443		enforce(i1d.length == i2d.length,
	443	dstatsEnforce(i1d.length == i2d.length,
444	444	"Ranges must be same length for Kendall correlation.");
445	445
…	…
456	456	*/
457	457	double kendallCorDestructive(T, U)(T[] input1, U[] input2) {
458		enforce(input1.length == input2.length,
	458	dstatsEnforce(input1.length == input2.length,
459	459	"Ranges must be same length for Kendall correlation.");
460	460	return kendallCorDestructiveLowLevel(input1, input2).field[0];
…	…
569	569	m1++;
570	570	} else {
571		enforce(0, "Can't compute Kendall's Tau with NaNs.");
	571	dstatsEnforce(0, "Can't compute Kendall's Tau with NaNs.");
572	572	}
573	573	} else if(input2[i] < input2[j]) {
…	…
579	579	m1++;
580	580	} else {
581		enforce(0, "Can't compute Kendall's Tau with NaNs.");
	581	dstatsEnforce(0, "Can't compute Kendall's Tau with NaNs.");
582	582	}
583	583	} else if(input2[i] == input2[j]) {
…	…
589	589	m1++;
590	590	} else {
591		enforce(0, "Can't compute Kendall's Tau with NaNs.");
	591	dstatsEnforce(0, "Can't compute Kendall's Tau with NaNs.");
592	592	}
593	593
594	594	} else {
595		enforce(0, "Can't compute Kendall's Tau with NaNs.");
	595	dstatsEnforce(0, "Can't compute Kendall's Tau with NaNs.");
596	596	}
597	597	}

trunk/distrib.d

r163	r164
185	185	///
186	186	double uniformPDF(double X, double lower, double upper) {
187		enforce(X >= lower, "Can't have X < lower bound in uniform distribution.");
188		enforce(X <= upper, "Can't have X > upper bound in uniform distribution.");
	187	dstatsEnforce(X >= lower, "Can't have X < lower bound in uniform distribution.");
	188	dstatsEnforce(X <= upper, "Can't have X > upper bound in uniform distribution.");
189	189	return 1.0L / (upper - lower);
190	190	}
…	…
192	192	///
193	193	double uniformCDF(double X, double lower, double upper) {
194		enforce(X >= lower, "Can't have X < lower bound in uniform distribution.");
195		enforce(X <= upper, "Can't have X > upper bound in uniform distribution.");
	194	dstatsEnforce(X >= lower, "Can't have X < lower bound in uniform distribution.");
	195	dstatsEnforce(X <= upper, "Can't have X > upper bound in uniform distribution.");
196	196
197	197	return (X - lower) / (upper - lower);
…	…
200	200	///
201	201	double uniformCDFR(double X, double lower, double upper) {
202		enforce(X >= lower, "Can't have X < lower bound in uniform distribution.");
203		enforce(X <= upper, "Can't have X > upper bound in uniform distribution.");
	202	dstatsEnforce(X >= lower, "Can't have X < lower bound in uniform distribution.");
	203	dstatsEnforce(X <= upper, "Can't have X > upper bound in uniform distribution.");
204	204
205	205	return (upper - X) / (upper - lower);
…	…
208	208	///
209	209	double poissonPMF(ulong k, double lambda) {
210		enforce(lambda > 0, "Cannot have a Poisson with lambda <= 0 or nan.");
	210	dstatsEnforce(lambda > 0, "Cannot have a Poisson with lambda <= 0 or nan.");
211	211
212	212	return exp(cast(double) k * log(lambda) -
…	…
234	234	/*P(K <= k) where K is r.v./
235	235	double poissonCDF(ulong k, double lambda) {
236		enforce(lambda > 0, "Cannot have a poisson with lambda <= 0 or nan.");
	236	dstatsEnforce(lambda > 0, "Cannot have a poisson with lambda <= 0 or nan.");
237	237
238	238	return (max(k, lambda) >= POISSON_NORMAL) ?
…	…
286	286	/*P(K >= k) where K is r.v./
287	287	double poissonCDFR(ulong k, double lambda) {
288		enforce(lambda > 0, "Can't have a poisson with lambda <= 0 or nan.");
	288	dstatsEnforce(lambda > 0, "Can't have a poisson with lambda <= 0 or nan.");
289	289
290	290	return (max(k, lambda) >= POISSON_NORMAL) ?
…	…
329	329	* is closest to pVal is used.*/
330	330	uint invPoissonCDF(double pVal, double lambda) {
331		enforce(lambda > 0, "Cannot have a poisson with lambda <= 0 or nan.");
332		enforce(pVal >= 0 && pVal <= 1, "P-values must be between 0, 1.");
	331	dstatsEnforce(lambda > 0, "Cannot have a poisson with lambda <= 0 or nan.");
	332	dstatsEnforce(pVal >= 0 && pVal <= 1, "P-values must be between 0, 1.");
333	333
334	334	// Use normal approximation to get approx answer, then brute force search.
…	…
381	381	///
382	382	double binomialPMF(ulong k, ulong n, double p) {
383		enforce(k <= n, "k cannot be > n in binomial distribution.");
384		enforce(p >= 0 && p <= 1, "p must be between 0, 1 in binomial distribution.");
	383	dstatsEnforce(k <= n, "k cannot be > n in binomial distribution.");
	384	dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 in binomial distribution.");
385	385	return exp(logNcomb(n, k) + k * log(p) + (n - k) * log(1 - p));
386	386	}
…	…
416	416	///P(K <= k) where K is random variable.
417	417	double binomialCDF(ulong k, ulong n, double p) {
418		enforce(k <= n, "k cannot be > n in binomial distribution.");
419		enforce(p >= 0 && p <= 1, "p must be between 0, 1 in binomial distribution.");
	418	dstatsEnforce(k <= n, "k cannot be > n in binomial distribution.");
	419	dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 in binomial distribution.");
420	420
421	421	if(k == n) {
…	…
496	496	///P(K >= k) where K is random variable.
497	497	double binomialCDFR(ulong k, ulong n, double p) {
498		enforce(k <= n, "k cannot be > n in binomial distribution.");
499		enforce(p >= 0 && p <= 1, "p must be between 0, 1 in binomial distribution.");
	498	dstatsEnforce(k <= n, "k cannot be > n in binomial distribution.");
	499	dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 in binomial distribution.");
500	500
501	501	if(k == 0) {
…	…
570	570	* is closest to pVal is used.*/
571	571	uint invBinomialCDF(double pVal, uint n, double p) {
572		enforce(pVal >= 0 && pVal <= 1, "p-values must be between 0, 1.");
573		enforce(p >= 0 && p <= 1, "p must be between 0, 1 in binomial distribution.");
	572	dstatsEnforce(pVal >= 0 && pVal <= 1, "p-values must be between 0, 1.");
	573	dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 in binomial distribution.");
574	574
575	575	// Use normal approximation to get approx answer, then brute force search.
…	…
653	653	// both speed and accuracy are "good enough" for most practical purposes.
654	654	double hypergeometricCDF(long x, long n1, long n2, long n) {
655		enforce(x <= n, "x must be <= n in hypergeometric distribution.");
656		enforce(n <= n1 + n2, "n must be <= n1 + n2 in hypergeometric distribution.");
657		enforce(x >= 0, "x must be >= 0 in hypergeometric distribution.");
	655	dstatsEnforce(x <= n, "x must be <= n in hypergeometric distribution.");
	656	dstatsEnforce(n <= n1 + n2, "n must be <= n1 + n2 in hypergeometric distribution.");
	657	dstatsEnforce(x >= 0, "x must be >= 0 in hypergeometric distribution.");
658	658
659	659	ulong expec = (n1 * n) / (n1 + n2);
…	…
746	746	///P(X >= x), where X is random variable.
747	747	double hypergeometricCDFR(ulong x, ulong n1, ulong n2, ulong n) {
748		enforce(x <= n, "x must be <= n in hypergeometric distribution.");
749		enforce(n <= n1 + n2, "n must be <= n1 + n2 in hypergeometric distribution.");
750		enforce(x >= 0, "x must be >= 0 in hypergeometric distribution.");
	748	dstatsEnforce(x <= n, "x must be <= n in hypergeometric distribution.");
	749	dstatsEnforce(n <= n1 + n2, "n must be <= n1 + n2 in hypergeometric distribution.");
	750	dstatsEnforce(x >= 0, "x must be >= 0 in hypergeometric distribution.");
751	751
752	752	return hypergeometricCDF(n - x, n2, n1, n);
…	…
767	767
768	768	double hyperExact(ulong x, ulong n1, ulong n2, ulong n, ulong startAt = 0) {
769		enforce(x <= n, "x must be <= n in hypergeometric distribution.");
770		enforce(n <= n1 + n2, "n must be <= n1 + n2 in hypergeometric distribution.");
771		enforce(x >= 0, "x must be >= 0 in hypergeometric distribution.");
	769	dstatsEnforce(x <= n, "x must be <= n in hypergeometric distribution.");
	770	dstatsEnforce(n <= n1 + n2, "n must be <= n1 + n2 in hypergeometric distribution.");
	771	dstatsEnforce(x >= 0, "x must be >= 0 in hypergeometric distribution.");
772	772
773	773	immutable double constPart = logFactorial(n1) + logFactorial(n2) +
…	…
800	800	///
801	801	double chiSquarePDF(double x, double v) {
802		enforce(x >= 0, "x must be >= 0 in chi-square distribution.");
803		enforce(v >= 1.0, "Must have at least 1 degree of freedom for chi-square.");
	802	dstatsEnforce(x >= 0, "x must be >= 0 in chi-square distribution.");
	803	dstatsEnforce(v >= 1.0, "Must have at least 1 degree of freedom for chi-square.");
804	804
805	805	// Calculate in log space for stability.
…	…
839	839	*/
840	840	double chiSquareCDF(double x, double v) {
841		enforce(x >= 0, "x must be >= 0 in chi-square distribution.");
842		enforce(v >= 1.0, "Must have at least 1 degree of freedom for chi-square.");
	841	dstatsEnforce(x >= 0, "x must be >= 0 in chi-square distribution.");
	842	dstatsEnforce(v >= 1.0, "Must have at least 1 degree of freedom for chi-square.");
843	843	return gammaIncomplete( 0.5v, 0.5x);
844	844	}
…	…
846	846	///
847	847	double chiSquareCDFR(double x, double v) {
848		enforce(x >= 0, "x must be >= 0 in chi-square distribution.");
849		enforce(v >= 1.0, "Must have at least 1 degree of freedom for chi-square.");
	848	dstatsEnforce(x >= 0, "x must be >= 0 in chi-square distribution.");
	849	dstatsEnforce(v >= 1.0, "Must have at least 1 degree of freedom for chi-square.");
850	850	return gammaIncompleteCompl( 0.5Lv, 0.5Lx );
851	851	}
…	…
864	864	*/
865	865	double invChiSquareCDFR(double v, double p) {
866		enforce(v >= 1.0, "Must have at least 1 degree of freedom for chi-square.");
867		enforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
	866	dstatsEnforce(v >= 1.0, "Must have at least 1 degree of freedom for chi-square.");
	867	dstatsEnforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
868	868	return 2.0 * gammaIncompleteComplInv( 0.5*v, p);
869	869	}
…	…
878	878	///
879	879	double normalPDF(double x, double mean = 0, double sd = 1) {
880		enforce(sd > 0, "Standard deviation must be > 0 for normal distribution.");
	880	dstatsEnforce(sd > 0, "Standard deviation must be > 0 for normal distribution.");
881	881	double dev = x - mean;
882	882	return exp(-(dev * dev) / (2 * sd * sd)) / (sd * SQ2PI);
…	…
889	889	///P(X < x) for normal distribution where X is random var.
890	890	double normalCDF(double x, double mean = 0, double stdev = 1) {
891		enforce(stdev > 0, "Standard deviation must be > 0 for normal distribution.");
	891	dstatsEnforce(stdev > 0, "Standard deviation must be > 0 for normal distribution.");
892	892
893	893	// Using a slightly non-obvious implementation in terms of erfc because
…	…
907	907	///P(X > x) for normal distribution where X is random var.
908	908	double normalCDFR(double x, double mean = 0, double stdev = 1) {
909		enforce(stdev > 0, "Standard deviation must be > 0 for normal distribution.");
	909	dstatsEnforce(stdev > 0, "Standard deviation must be > 0 for normal distribution.");
910	910
911	911	double Z = (x - mean) / stdev;
…	…
1027	1027	*/
1028	1028	double invNormalCDF(double p, double mean = 0, double sd = 1) {
1029		enforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
1030		enforce(sd > 0, "Standard deviation must be > 0 for normal distribution.");
	1029	dstatsEnforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
	1030	dstatsEnforce(sd > 0, "Standard deviation must be > 0 for normal distribution.");
1031	1031
1032	1032	if (p == 0.0L) {
…	…
1105	1105	///
1106	1106	double logNormalPDF(double x, double mu = 0, double sigma = 1) {
1107		enforce(sigma > 0, "sigma must be > 0 for log-normal distribution.");
	1107	dstatsEnforce(sigma > 0, "sigma must be > 0 for log-normal distribution.");
1108	1108
1109	1109	immutable mulTerm = 1.0L / (x * sigma * SQ2PI);
…	…
1122	1122	///
1123	1123	double logNormalCDF(double x, double mu = 0, double sigma = 1) {
1124		enforce(sigma > 0, "sigma must be > 0 for log-normal distribution.");
	1124	dstatsEnforce(sigma > 0, "sigma must be > 0 for log-normal distribution.");
1125	1125
1126	1126	return 0.5L + 0.5L * erf((log(x) - mu) / (sigma * SQRT2));
…	…
1134	1134	///
1135	1135	double logNormalCDFR(double x, double mu = 0, double sigma = 1) {
1136		enforce(sigma > 0, "sigma must be > 0 for log-normal distribution.");
	1136	dstatsEnforce(sigma > 0, "sigma must be > 0 for log-normal distribution.");
1137	1137
1138	1138	return 0.5L - 0.5L * erf((log(x) - mu) / (sigma * SQRT2));
…	…
1147	1147	///
1148	1148	double weibullPDF(double x, double shape, double scale = 1) {
1149		enforce(shape > 0, "shape must be > 0 for weibull distribution.");
1150		enforce(scale > 0, "scale must be > 0 for weibull distribution.");
	1149	dstatsEnforce(shape > 0, "shape must be > 0 for weibull distribution.");
	1150	dstatsEnforce(scale > 0, "scale must be > 0 for weibull distribution.");
1151	1151
1152	1152	if(x < 0) {
…	…
1163	1163	///
1164	1164	double weibullCDF(double x, double shape, double scale = 1) {
1165		enforce(shape > 0, "shape must be > 0 for weibull distribution.");
1166		enforce(scale > 0, "scale must be > 0 for weibull distribution.");
	1165	dstatsEnforce(shape > 0, "shape must be > 0 for weibull distribution.");
	1166	dstatsEnforce(scale > 0, "scale must be > 0 for weibull distribution.");
1167	1167
1168	1168	double exponent = pow(x / scale, shape);
…	…
1176	1176	///
1177	1177	double weibullCDFR(double x, double shape, double scale = 1) {
1178		enforce(shape > 0, "shape must be > 0 for weibull distribution.");
1179		enforce(scale > 0, "scale must be > 0 for weibull distribution.");
	1178	dstatsEnforce(shape > 0, "shape must be > 0 for weibull distribution.");
	1179	dstatsEnforce(scale > 0, "scale must be > 0 for weibull distribution.");
1180	1180
1181	1181	double exponent = pow(x / scale, shape);
…	…
1204	1204	///
1205	1205	double studentsTCDF(double t, double df) {
1206		enforce(df > 0, "Student's T must have >0 degrees of freedom.");
	1206	dstatsEnforce(df > 0, "Student's T must have >0 degrees of freedom.");
1207	1207
1208	1208	double x = (t + sqrt(t * t + df)) / (2 * sqrt(t * t + df));
…	…
1237	1237	*/
1238	1238	double invStudentsTCDF(double p, double df) {
1239		enforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
1240		enforce(df > 0, "Student's T must have >0 degrees of freedom.");
	1239	dstatsEnforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
	1240	dstatsEnforce(df > 0, "Student's T must have >0 degrees of freedom.");
1241	1241
1242	1242	if (p==0) return -double.infinity;
…	…
1302	1302	*/
1303	1303	double fisherCDF(double x, double df1, double df2) {
1304		enforce(df1 > 0 && df2 > 0,
	1304	dstatsEnforce(df1 > 0 && df2 > 0,
1305	1305	"Fisher distribution must have >0 degrees of freedom.");
1306		enforce(x >= 0, "x must be >=0 for Fisher distribution.");
	1306	dstatsEnforce(x >= 0, "x must be >=0 for Fisher distribution.");
1307	1307
1308	1308	double a = cast(double)(df1);
…	…
1315	1315	/** ditto */
1316	1316	double fisherCDFR(double x, double df1, double df2) {
1317		enforce(df1 > 0 && df2 > 0,
	1317	dstatsEnforce(df1 > 0 && df2 > 0,
1318	1318	"Fisher distribution must have >0 degrees of freedom.");
1319		enforce(x >= 0, "x must be >=0 for Fisher distribution.");
	1319	dstatsEnforce(x >= 0, "x must be >=0 for Fisher distribution.");
1320	1320
1321	1321	double a = cast(double)(df1);
…	…
1345	1345	*/
1346	1346	double invFisherCDFR(double df1, double df2, double p ) {
1347		enforce(df1 > 0 && df2 > 0,
	1347	dstatsEnforce(df1 > 0 && df2 > 0,
1348	1348	"Fisher distribution must have >0 degrees of freedom.");
1349		enforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
	1349	dstatsEnforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
1350	1350
1351	1351	double a = df1;
…	…
1377	1377	///
1378	1378	double negBinomPMF(ulong k, ulong n, double p) {
1379		enforce(p >= 0 && p <= 1,
	1379	dstatsEnforce(p >= 0 && p <= 1,
1380	1380	"p must be between 0, 1 for negative binomial distribution.");
1381	1381
…	…
1414	1414	*/
1415	1415	double negBinomCDF(ulong k, ulong n, double p ) {
1416		enforce(p >= 0 && p <= 1,
	1416	dstatsEnforce(p >= 0 && p <= 1,
1417	1417	"p must be between 0, 1 for negative binomial distribution.");
1418	1418	if ( k == 0 ) return pow(p, cast(double) n);
…	…
1428	1428	/*Probability that k or more failures precede the nth success./
1429	1429	double negBinomCDFR(ulong k, ulong n, double p) {
1430		enforce(p >= 0 && p <= 1,
	1430	dstatsEnforce(p >= 0 && p <= 1,
1431	1431	"p must be between 0, 1 for negative binomial distribution.");
1432	1432
…	…
1443	1443	///
1444	1444	ulong invNegBinomCDF(double pVal, ulong n, double p) {
1445		enforce(p >= 0 && p <= 1,
	1445	dstatsEnforce(p >= 0 && p <= 1,
1446	1446	"p must be between 0, 1 for negative binomial distribution.");
1447		enforce(pVal >= 0 && pVal <= 1,
	1447	dstatsEnforce(pVal >= 0 && pVal <= 1,
1448	1448	"P-values must be between 0, 1.");
1449	1449
…	…
1540	1540	*/
1541	1541	double gammaCDF(double x, double a, double b) {
1542		enforce(x > 0, "x must be >0 in gamma distribution.");
1543		enforce(a > 0, "a must be >0 in gamma distribution.");
1544		enforce(b > 0, "b must be >0 in gamma distribution.");
	1542	dstatsEnforce(x > 0, "x must be >0 in gamma distribution.");
	1543	dstatsEnforce(a > 0, "a must be >0 in gamma distribution.");
	1544	dstatsEnforce(b > 0, "b must be >0 in gamma distribution.");
1545	1545
1546	1546	return gammaIncomplete(b, a*x);
…	…
1549	1549	/** ditto */
1550	1550	double gammaCDFR(double x, double a, double b) {
1551		enforce(x > 0, "x must be >0 in gamma distribution.");
1552		enforce(a > 0, "a must be >0 in gamma distribution.");
1553		enforce(b > 0, "b must be >0 in gamma distribution.");
	1551	dstatsEnforce(x > 0, "x must be >0 in gamma distribution.");
	1552	dstatsEnforce(a > 0, "a must be >0 in gamma distribution.");
	1553	dstatsEnforce(b > 0, "b must be >0 in gamma distribution.");
1554	1554
1555	1555	return gammaIncompleteCompl( b, a * x );
…	…
1558	1558	///
1559	1559	double invGammaCDFR(double p, double a, double b) {
1560		enforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
1561		enforce(a > 0, "a must be >0 in gamma distribution.");
1562		enforce(b > 0, "b must be >0 in gamma distribution.");
	1560	dstatsEnforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
	1561	dstatsEnforce(a > 0, "a must be >0 in gamma distribution.");
	1562	dstatsEnforce(b > 0, "b must be >0 in gamma distribution.");
1563	1563
1564	1564	double ax = gammaIncompleteComplInv(b, p);
…	…
1573	1573	///
1574	1574	double cauchyPDF(double X, double X0 = 0, double gamma = 1) {
1575		enforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
	1575	dstatsEnforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
1576	1576
1577	1577	double toSquare = (X - X0) / gamma;
…	…
1588	1588	///
1589	1589	double cauchyCDF(double X, double X0 = 0, double gamma = 1) {
1590		enforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
	1590	dstatsEnforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
1591	1591
1592	1592	return M_1_PI * atan((X - X0) / gamma) + 0.5L;
…	…
1601	1601	///
1602	1602	double cauchyCDFR(double X, double X0 = 0, double gamma = 1) {
1603		enforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
	1603	dstatsEnforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
1604	1604
1605	1605	return M_1_PI * atan((X0 - X) / gamma) + 0.5L;
…	…
1615	1615	///
1616	1616	double invCauchyCDF(double p, double X0 = 0, double gamma = 1) {
1617		enforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
1618		enforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
	1617	dstatsEnforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
	1618	dstatsEnforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
1619	1619
1620	1620	return X0 + gamma * tan(PI * (p - 0.5L));
…	…
1637	1637	///
1638	1638	double laplacePDF(double x, double mu = 0, double b = 1) {
1639		enforce(b > 0, "b must be > 0 for laplace distribution.");
	1639	dstatsEnforce(b > 0, "b must be > 0 for laplace distribution.");
1640	1640
1641	1641	return (exp(-abs(x - mu) / b)) / (2 * b);
…	…
1650	1650	///
1651	1651	double laplaceCDF(double X, double mu = 0, double b = 1) {
1652		enforce(b > 0, "b must be > 0 for laplace distribution.");
	1652	dstatsEnforce(b > 0, "b must be > 0 for laplace distribution.");
1653	1653
1654	1654	double diff = (X - mu);
…	…
1666	1666	///
1667	1667	double laplaceCDFR(double X, double mu = 0, double b = 1) {
1668		enforce(b > 0, "b must be > 0 for laplace distribution.");
	1668	dstatsEnforce(b > 0, "b must be > 0 for laplace distribution.");
1669	1669
1670	1670	double diff = (mu - X);
…	…
1683	1683	///
1684	1684	double invLaplaceCDF(double p, double mu = 0, double b = 1) {
1685		enforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
1686		enforce(b > 0, "b must be > 0 for laplace distribution.");
	1685	dstatsEnforce(p >= 0 && p <= 1, "P-values must be between 0, 1.");
	1686	dstatsEnforce(b > 0, "b must be > 0 for laplace distribution.");
1687	1687
1688	1688	double p05 = p - 0.5L;
…	…
1700	1700	/*Kolmogorov distribution. Used in Kolmogorov-Smirnov testing./
1701	1701	double kolmDist(double X) {
1702		enforce(X >= 0, "X must be >= 0 for Kolmogorov distribution.");
	1702	dstatsEnforce(X >= 0, "X must be >= 0 for Kolmogorov distribution.");
1703	1703
1704	1704	if(X == 0) {

trunk/infotheory.d

r159	r164
57	57	double entropyCounts(T)(T data)
58	58	if(isForwardRange!(T) && doubleInput!(T)) {
59		auto save = data;
	59	auto save = data.save();
60	60	return entropyCounts(save, sum!(T, double)(data));
61	61	}
…	…
262	262	} else static if(isArray!(T)) {
263	263	enum bool NeedsHeap = false;
264		} else static if(is(T == Joint!(typeof(T.init.tupleof)))
	264	} else static if(is(Joint!(typeof(T.init.tupleof)))
	265	&& is(T == Joint!(typeof(T.init.tupleof)))
265	266	&& allSatisfy!(isArray, typeof(T.init.tupleof))) {
266	267	enum bool NeedsHeap = false;
…	…
272	273	unittest {
273	274	auto foo = map!("a.dup")(cast(uint[][]) [[1]]);
274		auto bar = map!("a + 2")(~~cast(uint[]) [1,2,3~~]);
	275	auto bar = map!("a + 2")([1,2,3][]);
275	276	static assert(NeedsHeap!(typeof(foo)));
276	277	static assert(!NeedsHeap!(typeof(bar)));

trunk/random.d

r159	r164
122	122	public import std.random; //For uniform distrib.
123	123
124		import dstats.alloc;
	124	import dstats.alloc, dstats.base;
125	125
126	126	version(unittest) {
…	…
200	200	normData = *lastNormPtr;
201	201	*lastNormPtr = temp;
	202	}
	203
	204	typeof(this) save() @property {
	205	return this;
	206	}
	207
	208	// Doesn't do anything particularly useful. It's more a workaround for
	209	// bug 4305.
	210	typeof(randFun(args)) moveFront() {
	211	return frontElem;
202	212	}
203	213	}
…	…
263	273	///
264	274	double rNorm(RGen = Random)(double mean, double sd, ref RGen gen = rndGen) {
265		enforce(sd > 0, "Standard deviation must be > 0 for rNorm.");
	275	dstatsEnforce(sd > 0, "Standard deviation must be > 0 for rNorm.");
266	276
267	277	double lr = lastNorm;
…	…
298	308	///
299	309	double rCauchy(RGen = Random)(double X0, double gamma, ref RGen gen = rndGen) {
300		enforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
	310	dstatsEnforce(gamma > 0, "gamma must be > 0 for Cauchy distribution.");
301	311
302	312	return (rNorm(0, 1, gen) / rNorm(0, 1, gen)) * gamma + X0;
…	…
318	328	///
319	329	double rStudentT(RGen = Random)(double df, ref RGen gen = rndGen) {
320		enforce(df > 0, "Student's T distribution must have >0 degrees of freedom.");
	330	dstatsEnforce(df > 0, "Student's T distribution must have >0 degrees of freedom.");
321	331
322	332	double N = rNorm(0, 1, gen);
…	…
341	351	///
342	352	double rFisher(RGen = Random)(double df1, double df2, ref RGen gen = rndGen) {
343		enforce(df1 > 0 && df2 > 0,
	353	dstatsEnforce(df1 > 0 && df2 > 0,
344	354	"df1 and df2 must be >0 for the Fisher distribution.");
345	355
…	…
362	372	///
363	373	double rChiSquare(RGen = Random)(double df, ref RGen gen = rndGen) {
364		enforce(df > 0, "df must be > 0 for chiSquare distribution.");
	374	dstatsEnforce(df > 0, "df must be > 0 for chiSquare distribution.");
365	375
366	376	return 2.0 * stdGamma(df / 2.0L, gen);
…	…
384	394	///
385	395	int rPoisson(RGen = Random)(double lam, ref RGen gen = rndGen) {
386		enforce(lam > 0, "lambda must be >0 for Poisson distribution.");
	396	dstatsEnforce(lam > 0, "lambda must be >0 for Poisson distribution.");
387	397
388	398	static int poissonMult(ref RGen gen, double lam) {
…	…
466	476	///
467	477	int rBernoulli(RGen = Random)(double P = 0.5, ref RGen gen = rndGen) {
468		enforce(P >= 0 && P <= 1, "P must be between 0, 1 for Bernoulli distribution.");
	478	dstatsEnforce(P >= 0 && P <= 1, "P must be between 0, 1 for Bernoulli distribution.");
469	479
470	480	double pVal = uniform(0.0L, 1.0L, gen);
…	…
662	672	///
663	673	int rBinomial(RGen = Random)(int n, double p, ref RGen gen = rndGen) {
664		enforce(n >= 0, "n must be >= 0 for binomial distribution.");
665		enforce(p >= 0 && p <= 1, "p must be between 0, 1 for binomial distribution.");
	674	dstatsEnforce(n >= 0, "n must be >= 0 for binomial distribution.");
	675	dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 for binomial distribution.");
666	676
667	677	if (p <= 0.5) {
…	…
776	786	///
777	787	int rHypergeometric(RGen = Random)(int n1, int n2, int n, ref RGen gen = rndGen) {
778		enforce(n <= n1 + n2, "n must be <= n1 + n2 for hypergeometric distribution.");
779		enforce(n1 >= 0 && n2 >= 0 && n >= 0,
	788	dstatsEnforce(n <= n1 + n2, "n must be <= n1 + n2 for hypergeometric distribution.");
	789	dstatsEnforce(n1 >= 0 && n2 >= 0 && n >= 0,
780	790	"n, n1, n2 must be >= 0 for hypergeometric distribution.");
781	791
…	…
843	853
844	854	int rGeometric(RGen = Random)(double p, ref RGen gen = rndGen) {
845		enforce(p >= 0 && p <= 1, "p must be between 0, 1 for geometric distribution.");
	855	dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 for geometric distribution.");
846	856
847	857	if (p >= 0.333333333333333333333333) {
…	…
879	889	///
880	890	int rNegBinom(RGen = Random)(double n, double p, ref RGen gen = rndGen) {
881		enforce(n >= 0, "n must be >= 0 for negative binomial distribution.");
882		enforce(p >= 0 && p <= 1,
	891	dstatsEnforce(n >= 0, "n must be >= 0 for negative binomial distribution.");
	892	dstatsEnforce(p >= 0 && p <= 1,
883	893	"p must be between 0, 1 for negative binomial distribution.");
884	894
…	…
911	921	///
912	922	double rLaplace(RGen = Random)(double mu = 0, double b = 1, ref RGen gen = rndGen) {
913		enforce(b > 0, "b must be > 0 for Laplace distribution.");
	923	dstatsEnforce(b > 0, "b must be > 0 for Laplace distribution.");
914	924
915	925	double p = uniform(0.0L, 1.0L, gen);
…	…
935	945	///
936	946	double rExponential(RGen = Random)(double lambda, ref RGen gen = rndGen) {
937		enforce(lambda > 0, "lambda must be > 0 for exponential distribution.");
	947	dstatsEnforce(lambda > 0, "lambda must be > 0 for exponential distribution.");
938	948
939	949	double p = uniform(0.0L, 1.0L, gen);
…	…
997	1007	///
998	1008	double rGamma(RGen = Random)(double a, double b, ref RGen gen = rndGen) {
999		enforce(a > 0, "a must be > 0 for gamma distribution.");
1000		enforce(b > 0, "b must be > 0 for gamma distribution.");
	1009	dstatsEnforce(a > 0, "a must be > 0 for gamma distribution.");
	1010	dstatsEnforce(b > 0, "b must be > 0 for gamma distribution.");
1001	1011
1002	1012	return stdGamma(b, gen) / a;
…	…
1019	1029	///
1020	1030	double rBeta(RGen = Random)(double a, double b, ref RGen gen = rndGen) {
1021		enforce(a > 0, "a must be > 0 for beta distribution.");
1022		enforce(b > 0, "b must be > 0 for beta distribution.");
	1031	dstatsEnforce(a > 0, "a must be > 0 for beta distribution.");
	1032	dstatsEnforce(b > 0, "b must be > 0 for beta distribution.");
1023	1033
1024	1034	double Ga = void, Gb = void;
…	…
1091	1101	///
1092	1102	double rLogistic(RGen = Random)(double loc, double scale, ref RGen gen = rndGen) {
1093		enforce(scale > 0, "scale must be > 0 for logistic distribution.");
	1103	dstatsEnforce(scale > 0, "scale must be > 0 for logistic distribution.");
1094	1104
1095	1105	double U = uniform(0.0L, 1.0L, gen);
…	…
1113	1123	///
1114	1124	double rLogNorm(RGen = Random)(double mu, double sigma, ref RGen gen = rndGen) {
1115		enforce(sigma > 0, "sigma must be > 0 for log-normal distribution.");
	1125	dstatsEnforce(sigma > 0, "sigma must be > 0 for log-normal distribution.");
1116	1126
1117	1127	return exp(rNorm(mu, sigma, gen));
…	…
1135	1145	///
1136	1146	double rWeibull(RGen = Random)(double shape, double scale = 1, ref RGen gen = rndGen) {
1137		enforce(shape > 0, "shape must be > 0 for weibull distribution.");
1138		enforce(scale > 0, "scale must be > 0 for weibull distribution.");
	1147	dstatsEnforce(shape > 0, "shape must be > 0 for weibull distribution.");
	1148	dstatsEnforce(scale > 0, "scale must be > 0 for weibull distribution.");
1139	1149
1140	1150	return pow(rExponential(1, gen), 1. / shape) * scale;
…	…
1152	1162	///
1153	1163	double rWald(RGen = Random)(double mu, double lambda, ref RGen gen = rndGen) {
1154		enforce(mu > 0, "mu must be > 0 for Wald distribution.");
1155		enforce(lambda > 0, "lambda must be > 0 for Wald distribution.");
	1164	dstatsEnforce(mu > 0, "mu must be > 0 for Wald distribution.");
	1165	dstatsEnforce(lambda > 0, "lambda must be > 0 for Wald distribution.");
1156	1166
1157	1167	alias mu mean;
…	…
1187	1197	///
1188	1198	double rRayleigh(RGen = Random)(double mode, ref RGen gen = rndGen) {
1189		enforce(mode > 0, "mode must be > 0 for Rayleigh distribution.");
	1199	dstatsEnforce(mode > 0, "mode must be > 0 for Rayleigh distribution.");
1190	1200
1191	1201	return modesqrt(-2.0 log(1.0 - uniform(0.0L, 1.0L, gen)));

trunk/regress.d

r145	r164
36	36
37	37	import std.math, std.algorithm, std.traits, std.array, std.traits, std.contracts,
38		dstats.alloc, std.range, std.conv, dstats.distrib, dstats.cor, dstats.base;
	38	std.typetuple;
	39
	40	import dstats.alloc, std.range, std.conv, dstats.distrib, dstats.cor, dstats.base;
39	41
40	42	private void enforceConfidence(double conf) {
41		enforce(conf >= 0 && conf <= 1,
	43	dstatsEnforce(conf >= 0 && conf <= 1,
42	44	"Confidence intervals must be between 0 and 1.");
43	45	}
…	…
72	74	}
73	75
74		bool empty() {
	76	typeof(this) save() @property {
	77	return this;
	78	}
	79
	80	bool empty() @property {
75	81	return range.empty;
76	82	}
…	…
269	275	size_t i = 0;
270	276	foreach(elem; X) {
271		sum += cast(double) elem.front * betas[i];
	277	double frnt = elem.front();
	278	sum += frnt * betas[i];
272	279	i++;
273	280	}
…	…
277	284	this(F[] betas, U Y, R X) {
278	285	static if(is(typeof(X.length))) {
279		enforce(X.length == betas.length,
	286	dstatsEnforce(X.length == betas.length,
280	287	"Betas and X must have same length for residuals.");
281	288	} else {
282		enforce(walkLength(X) == betas.length,
	289	dstatsEnforce(walkLength(X) == betas.length,
283	290	"Betas and X must have same length for residuals.");
284	291	}
…	…
434	441	static if(isForwardRange!(T[0]) && isForwardRange!(typeof(XIn[0].front())) &&
435	442	T.length == 1) {
	443
	444	enum bool arrayX = true;
436	445	alias typeof(XIn[0].front) E;
437	446	E[] X = tempdup(XIn[0]);
438	447	} else static if(allSatisfy!(isForwardRange, T)) {
	448	enum bool arrayX = false;
439	449	alias XIn X;
440	450	} else {
…	…
445	455	double[][] xTx;
446	456	double[] xTy;
447		rangeMatrixMulTrans(xTy, xTx, Y, X);
	457
	458	typeof(X) xSaved;
	459	static if(arrayX) {
	460	xSaved = X.tempdup;
	461	foreach(ref elem; xSaved) {
	462	elem = elem.save;
	463	}
	464	} else {
	465	xSaved = saveAll(X).expand;
	466	}
	467
	468	rangeMatrixMulTrans(xTy, xTx, Y.save, X);
448	469	invert(xTx);
449	470	double[] betas = new double[X.length];

trunk/tests.d

r159	r164
32	32
33	33	import std.algorithm, std.functional, std.range, std.conv, std.math, std.traits,
34		std.contracts;
	34	std.contracts, std.typetuple;
35	35
36	36	import dstats.base, dstats.distrib, dstats.alloc, dstats.summary, dstats.sort,
…	…
43	43
44	44	private void enforceConfidence(double conf) {
45		enforce(conf >= 0 && conf <= 1,
	45	dstatsEnforce(conf >= 0 && conf <= 1,
46	46	"Confidence intervals must be between 0 and 1.");
47	47	}
…	…
149	149	if( (isSummary!T \|\| doubleIterable!T)) {
150	150	enforceConfidence(confLevel);
151		enforce(isFinite(testMean), "testMean must not be infinite or nan.");
	151	dstatsEnforce(isFinite(testMean), "testMean must not be infinite or nan.");
152	152
153	153	static if(isSummary!T) {
…	…
197	197	if( (doubleIterable!T \|\| isSummary!T) && (doubleIterable!U \|\| isSummary!U)) {
198	198	enforceConfidence(confLevel);
199		enforce(isFinite(testMean), "testMean must not be infinite or nan.");
	199	dstatsEnforce(isFinite(testMean), "testMean must not be infinite or nan.");
200	200
201	201	static if(isSummary!T) {
…	…
311	311	if( (isSummary!T \|\| doubleIterable!T) && (isSummary!U \|\| doubleIterable!U)) {
312	312	enforceConfidence(confLevel);
313		enforce(isFinite(testMean), "testMean cannot be infinite or nan.");
	313	dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or nan.");
314	314
315	315	static if(isSummary!T) {
…	…
428	428	if(doubleInput!(T) && doubleInput!(U) && isInputRange!T && isInputRange!U) {
429	429	enforceConfidence(confLevel);
430		enforce(isFinite(testMean), "testMean cannot be infinite or nan.");
	430	dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or nan.");
431	431
432	432	MeanSD msd;
…	…
466	466	if(isSummary!T) {
467	467	enforceConfidence(confLevel);
468		enforce(isFinite(testMean), "testMean cannot be infinite or nan.");
	468	dstatsEnforce(isFinite(testMean), "testMean cannot be infinite or nan.");
469	469
470	470	// Save typing.
…	…
640	640	auto res2 = fTest([thing1, thing2, thing3].dup);
641	641	assert(approxEqual(result.testStat, res2.testStat));
642		assert(~~result.p == res2.p~~);
	642	assert(approxEqual(result.p, res2.p));
643	643
644	644	thing1 = [2,7,1,8,2];
…	…
692	692	// The cast is to force conversions to double on alias this'd stuff
693	693	// like the Mean struct.
694		centers[i] = cast(double) central(category);
	694	centers[i] = cast(double) central(category.save);
695	695	}
696	696
…	…
1481	1481
1482	1482	static if(dstats.base.hasLength!T && dstats.base.hasLength!U) {
1483		enforce(before.length == after.length,
	1483	dstatsEnforce(before.length == after.length,
1484	1484	"Ranges must have same lengths for wilcoxonSignedRank.");
1485	1485
…	…
1827	1827	*/
1828	1828	double binomialTest(ulong k, ulong n, double p) {
1829		enforce(k <= n, "k must be <= n for binomial test.");
1830		enforce(p >= 0 && p <= 1, "p must be between 0, 1 for binomial test.");
	1829	dstatsEnforce(k <= n, "k must be <= n for binomial test.");
	1830	dstatsEnforce(p >= 0 && p <= 1, "p must be between 0, 1 for binomial test.");
1831	1831
1832	1832	enum epsilon = 1 - 1e-6; // Small but arbitrary constant to deal w/ rounding error.
…	…
2013	2013	if(doubleInput!(T) && doubleInput!(U)) {
2014	2014	if(countProp == Expected.PROPORTION) {
2015		enforce(isForwardRange!(U),
	2015	dstatsEnforce(isForwardRange!(U),
2016	2016	"Can't use expected proportions instead of counts with input ranges.");
2017	2017	}
…	…
2025	2025	double expectSum = 0;
2026	2026	multiplier = 0;
2027		auto obsCopy = observed;
2028		auto expCopy = expected;
	2027	auto obsCopy = observed.save;
	2028	auto expCopy = expected.save;
2029	2029	while(!obsCopy.empty && !expCopy.empty) {
2030	2030	multiplier += obsCopy.front;
…	…
2046	2046	// exist.
2047	2047	if(e == 0) {
2048		enforce(observed.front == 0,
	2048	dstatsEnforce(observed.front == 0,
2049	2049	"Can't have non-zero observed value w/ zero expected value.");
2050	2050	continue;
…	…
2297	2297	foreach(ri, range; ranges) {
2298	2298	size_t curLen = 0;
2299		foreach(elem; range) {
	2299	foreach(elem; range.save) {
2300	2300	colSums[ri] += cast(double) elem;
2301	2301	curLen++;
…	…
2389	2389	foreach(range; contingencyTable) {
2390	2390	foreach(elem; range) {
2391		enforce(elem >= 0,
	2391	dstatsEnforce(elem >= 0,
2392	2392	"Cannot have negative elements in a contingency table.");
2393	2393	}
…	…
2407	2407
2408	2408
2409		alias contingencyTable c;
	2409	alias contingencyTable c; // Save typing.
2410	2410	immutable oddsRatio = cast(double) c[0][0] * c[1][1] / c[0][1] / c[1][0];
2411	2411	if(alt == Alt.NONE) {
…	…
2514	2514	TestRes fisherExact(T)(const T[][] contingencyTable, Alt alt = Alt.TWOSIDE)
2515	2515	if(isIntegral!(T)) {
2516		enforce(contingencyTable.length == 2 &&
	2516	dstatsEnforce(contingencyTable.length == 2 &&
2517	2517	contingencyTable[0].length == 2 &&
2518	2518	contingencyTable[1].length == 2,
…	…
2888	2888	ConfInt pearsonCorTest(T)(double cor, T N, Alt alt = Alt.TWOSIDE, double confLevel = 0.95)
2889	2889	if(isNumeric!(T)) {
2890		enforce(N >= 0, "N must be >= 0 for pearsonCorTest.");
2891		enforce(cor > -1.0 \|\| approxEqual(cor, -1.0),
	2890	dstatsEnforce(N >= 0, "N must be >= 0 for pearsonCorTest.");
	2891	dstatsEnforce(cor > -1.0 \|\| approxEqual(cor, -1.0),
2892	2892	"Correlation must be between 0, 1.");
2893		enforce(cor < 1.0 \|\| approxEqual(cor, 1.0),
	2893	dstatsEnforce(cor < 1.0 \|\| approxEqual(cor, 1.0),
2894	2894	"Correlation must be between 0, 1.");
2895	2895	enforceConfidence(confLevel);
…	…
3367	3367	uint df = 0;
3368	3368	foreach(pVal; pVals) {
3369		enforce(pVal >= 0 && pVal <= 1,
	3369	dstatsEnforce(pVal >= 0 && pVal <= 1,
3370	3370	"P-values must be between 0, 1 for Fisher's Method.");
3371	3371	chiSq += log( cast(double) pVal);
…	…
3422	3422	qVals.length = pVals.length;
3423	3423	foreach(i, elem; pVals) {
3424		enforce(elem >= 0 && elem <= 1,
	3424	dstatsEnforce(elem >= 0 && elem <= 1,
3425	3425	"P-values must be between 0, 1 for falseDiscoveryRate.");
3426	3426	qVals[i] = cast(float) elem;
…	…
3534	3534
3535	3535	foreach(i, ref q; qVals) {
3536		enforce(q >= 0 && q <= 1,
	3536	dstatsEnforce(q >= 0 && q <= 1,
3537	3537	"P-values must be between 0, 1 for hochberg.");
3538	3538	q = min(1.0f, q * (cast(double) qVals.length - i));
…	…
3604	3604
3605	3605	foreach(i, ref q; qVals) {
3606		enforce(q >= 0 && q <= 1,
	3606	dstatsEnforce(q >= 0 && q <= 1,
3607	3607	"P-values must be between 0, 1 for holmBonferroni.");
3608	3608	q = min(1.0, q * (cast(double) qVals.length - i));

Download in other formats:

© Copyright 2004-2012 - The Sankaty Group, Inc. All rights reserved.

Make a Donation

Privacy Statement Terms of Service