Programm läuft nicht korrekt

Ramanujan

Ich habe den Fehler (bzw. vielmehr eine Löusng) gefunden. Verstehen tue ich ihn nicht.

Wenn ich die Eingabematrix so erzeuge, läuft alles durch (es gilt d=2):

Matrix<double> a_pre(n, d);
for (int i = 0; i < n; i++) {
	if (i < n/2) {
		a_pre(i, 0) = useful::random_LL(-50, 49) + useful::rnd(0,1);
		a_pre(i, 1) = useful::rnd(-.001, .001);
	} else {
		a_pre(i, 1) = useful::random_LL(-50, 49) + useful::rnd(0,1);
		a_pre(i, 0) = useful::rnd(-.001, .001);
	}
}

Wenn ich sie so erzeuge, stimmen die Ergebnisse nicht mehr:

Matrix<double> a_pre(n, d);
for (int i = 0; i < n; i++) {
	if (i < n/2) {
		a_pre(i, 0) = useful::rnd(-50, 50);
		a_pre(i, 1) = useful::rnd(-.001, .001);
	} else {
		a_pre(i, 1) = useful::rnd(-50, 50);
		a_pre(i, 0) = useful::rnd(-.001, .001);
	}
}

Hinter den useful-random-Funktionen steckt das folgende:

double useful::rnd() {
	return ((double) random_LL() / RANDOM_LL_MAX);
}

double useful::rnd(double min, double max) {
	return rnd() * (max - min) + min;
}

const unsigned long long useful::RANDOM_LL_MAX = -1ULL;
unsigned long long useful::random_LL() {
	long long rnd = 0;
	// 64bit-Zufallszahl erzeugen
	const unsigned lng = 8 * sizeof(long long);
	for (unsigned j = 0; j < lng; j += 15) {
		long long t = rand() & 32767; // Mindestgröße von RAND_MAX = 32767
		rnd |= t << j;
	}
	return rnd;
}

long long useful::random_LL(long long min, long long max) {
	unsigned long long size = max - min;
	unsigned long long rnd;
	do {
		rnd = random_LL();
	} while (rnd > RANDOM_LL_MAX - RANDOM_LL_MAX % size);
	return min + rnd % size;
}

In beiden Fällen liegen alle Werte der Matrix zwischen -50 und 50 und es gibt kein NaN-Werte.

Eigentlich sehe ich nur noch eine Möglichkeit für das Verhalten:
Bei dem Programm geht es darum, eine sogenannte Einbettung zu berechnen. Dazu wird eine implizit gespeicherte (s x n)-Matrix S erstellt, die als Einträge -1 und 1 hat. Die Matrix S*A/sqrt(s) hat nun Ähnliche Eigenschaften wie A. Unter anderem gilt, dass (A^T)*A (hier steht A^T für die Transponierte von A) etwa gleich ((SA)^T)*SA ist. Wählt man z.B. s = 34 sind so Abweichungen von höchstens 20% zu erwarten. Für die Matrix S hab ich einen eigenen pseudo-Zufallsgenerator geschrieben (der vor allem das Element Sij in O(1) Zeit berechnen kann). Jetzt könnte es natürlich trotzdem sein, dass irgendwelche stochastischen Abhängigkeiten zwischen der rand()-basierten rnd()-Funktion und der Matrix S bestehen, die dazu führen, dass das Ergebnis nicht mehr korrekt ist. Das ist extremst unwahrscheinlich, aber ich kann mir das anders nicht erklären. Zumindest würde das auch erklären, warum es bei Linux geht: Andere rand()-Implementierung => keine Abhängigkeit mehr => alles ok.

Oder sieht jemand noch einen anderen möglichen Grund (der möglicherweise weniger an den Haaren herbei gezogen ist)?

volkard

Kannste mit wenig Aufwand nen anderen Zufallszahlengenerator unter Windows testen? Haste boost oder sowas rumfliegen?

SeppJ

Zwei Bemerkungen:
1. Dir ist aber schon klar, dass die Summe zweier Gleichverteilungen keine Gleichverteilung mehr ist, oder? Das heißt deine Codes 1 und 2 erzeugen unterschiedlich verteilte Zahlen.
2. Bezüglich deines Verdachts: Darum nimmt man rand() auch nicht für kritische Anwendungen. Nimm Generatoren von geprüfter Qualität. In Numerikbibliotheken, Boost und der C++11-Standardbibliothek wirst du fündig. Zur allergrößten Not implementier selber einen bekannten Generator - ist nicht schwierig. Dann entfällt auch der ganze Aufwand, den du derzeit treibst, der Nachbehandlung der rand()-Zahlen. Und ja: Typische rand()-Implementierungen haben Korrelationen drin.

Ramanujan

@1.) Die Funktion useful::random_LL(-50, 49) gibt mir ganze Zahlen von -50 bis 49. useful::rnd(0,1) addiert Werte zwischen 0 und 1 hinzu. Somit ergibt sich eine Gleichverteilung von -50 bis 50. Ok, wenn man es genau nimmt, sind die ganzen Werte leicht häufiger, da meine Implementierung von rnd() auch 0 und 1 als Werte zulässt. Aber das fällt nicht ins Gewicht.

@2.) Ja, ok, ich werd mal gucken, was C++11 so bietet.

Edit:
Wenn ich sowas hier benutze:

Matrix<double> a_pre(n, d); 
for (int i = 0; i < n; i++) { 
     if (i < n/2) { 
         a_pre(i, 0) = useful::rnd(-50, 50); 
         for (int k = 0; k < rand() % 100; k++)
             rand();
         a_pre(i, 1) = useful::rnd(-.001, .001); 
     } else { 
         a_pre(i, 1) = useful::rnd(-50, 50); 
         for (int k = 0; k < rand() % 100; k++)
             rand();
         a_pre(i, 0) = useful::rnd(-.001, .001); 
     } 
}

Also quasi eine zufällige Anzahl von rand()-Aufrufen dazwischen schiebe. Läuft wieder alles korrekt. Naja, da lag es wohl wirklich an der Abhängigkeit. Möglicherweise hatten die Vorzeichen von S und A ja vorher ein irgendwie gemeinsames Muster.

SeppJ

Ramanujan schrieb:

@1.) Die Funktion useful::random_LL(-50, 49) gibt mir ganze Zahlen von -50 bis 49. useful::rnd(0,1) addiert Werte zwischen 0 und 1 hinzu. Somit ergibt sich eine Gleichverteilung von -50 bis 50.

Nein. Wie schon gesagt, ist die Summe von zwei Gleichverteilungen keine Gleichverteilung. Das ist ganz elementar zu wissen, wenn du mit Verteilungen arbeitest, wie die Summe zweier Zufallsvariablen verteilt ist^*. Schnapp dir unbedingt mal ein Buch über Stochastik! Das musst du wissen!

Beweis durch dein eigenes Gegenbeispiel: Die Zahl -50 hat die Wahrscheinlichkeit von random_LL==-50, also 1/100, mal die Wahrscheinlichkeit von rnd==0, also 1/2, also insgesamt 1/200. Gleiches gilt für die Zahl 50. Es sollten aber 1/101 sein, wenn es eine Gleichverteilung wäre. Ebenso ist die Wahrscheinlichkeit für eine Zahl N zwischen -50 und 50 nicht mehr 1/101, sondern die Wahrscheinlichkeit
P(N) = P(random_LL==N-1) * P(rnd==1) + P(random_LL==N)*P(rnd==0) = 1/200 + 1/200 = 1/100 != 1/101

qed.

^*: Ganz interessant und wohl eine der wichtigsten Erkenntnisse der Mathematik (riesige Teile der Naturwissenschaft beruhen darauf) ist übrigens, was die Summe unendlich vieler Zufallszahlen ist, (fast) egal wie die konkrete Verteilung dieser Zahlen aussieht. Das ist nämlich im Limes eine Gaussverteilung.

Ramanujan

rnd() liefert reelle Werte zwischen 0 und 1 (jeweils einschließlich). Die Wahrscheinlichkeit für den Wert für 0 ist 1/2^64 (wenn sizeof(long long) == 8 gilt). Da allerdings rnd() auch 1 als Ergebnis liefern kann, sind ganze Werte um 1/2^64 wahrscheinlicher (also theoretisch keine Gleichverteilung). Wenn rnd() von 0 (einschließlich) bis 1 (ausschließlich) liefern würde, wäre aber alles ok. Da double ziemlich ungenau ist (im Vergleich zu long long), spielt das aber keine Rolle, man kann also praktisch die Verteilung von einer Gleichverteilung nicht unterscheiden.

Im Allgemeinen hast du natürlich Recht: Wenn X und Y gleichverteilt, ist X+Y nicht zwingend ebenfalls gleichverteilt.

SeppJ

Ramanujan schrieb:

rnd() liefert reelle Werte zwischen 0 und 1 (jeweils einschließlich). Die Wahrscheinlichkeit für den Wert für 0 ist 1/2^64 (wenn sizeof(long long) == 8 gilt). Da allerdings rnd() auch 1 als Ergebnis liefern kann, sind ganze Werte um 1/2^64 wahrscheinlicher. Wenn rnd() von 0 (einschließlich) bis 1 (ausschließlich) liefern würde, wäre aber alles ok. Da double ziemlich ungenau ist (im Vergleich zu long long), spielt das aber keine Rolle.

Deine Argumentation ist nicht überzeugend, der Datentyp spielt keine Rolle. Aber wenn du nicht zur Kenntnis nehmen möchtest, wenn man dich auf Fehler hinweist, dann will ich dich auch nicht zu deinem Glück zwingen.

hustbaer

[Unsinn]
~~Freund Ramanujan...
Wenn du mit zwei Würfeln würfelst, und das Ergebnis zusammenzählst, wie ist dann die Verteilung? Ist das dann gleichverteilt?

Nein, ist es nicht. Genau so wenig wie die Werte in deinem Beispiel gleichverteilt sind.~~[/Unsinn]

Davon abgesehen ist die Implementierung von random_LL furchtbar. Sie verlässt sich auf rand() , und da rand() meist ein ganz einfacher LCG ist, hagelt es nur so kurze Serien in den unteren Bits.
Die Werte die random_LL produziert sind also voll von diversen statistischen Auffälligkeiten. Das selbe gilt dann natürlich auch für rnd() .

Und wie diese "statistischen Auffälligkeiten" aussehen ändert sich natürlich auch, wenn du 1x mehr oder weniger random_LL aufrufst.

SeppJ schrieb:

Ramanujan schrieb:

@1.) Die Funktion useful::random_LL(-50, 49) gibt mir ganze Zahlen von -50 bis 49. useful::rnd(0,1) addiert Werte zwischen 0 und 1 hinzu. Somit ergibt sich eine Gleichverteilung von -50 bis 50.

Nein. Wie schon gesagt, ist die Summe von zwei Gleichverteilungen keine Gleichverteilung.

In diesem Fall schon.

hustbaer schrieb:

Wenn du mit zwei Würfeln würfelst, und das Ergebnis zusammenzählst, wie ist dann die Verteilung? Ist das dann gleichverteilt?

Du bringst ein gutes Beispiel. Der Trick beruht ja darauf, dass sich z.B. 7 auf mehrere Arten schreiben lässt und diese nachher nicht mehr unterscheidbar sind.

Bei Rammi ist es aber anders. Er hat irgendeine Ganzzahl. Und er addiert eine Kommazahl zwischen 0 und 1 (ich hoffe mal, exklusive 1) hinzu. Vor- und Nachkommastellen quasi.

Wenn du die Zahl 14.62341 bekommst, gibt es genau 1 Möglichkeit für die Ganzzahl und genauso nur 1 für die Kommazahl. Jede erhaltene Zahl lässt sich so eindeutig zerlegen.

Folglich ist das hier ein Spezialfall, in dem die Summe zweier Gleichverteilungen gleichverteilt ist.

hustbaer

*self-facepalm*

Ja, haste recht. Sorry.
Der Rest meines Beitrags ist trotzdem nicht ganz verkehrt

Danke für die sehr sachliche und freundliche Korrektur!

distributor schrieb:

addiert eine Kommazahl zwischen 0 und 1 (ich hoffe mal, exklusive 1) hinzu

Ne, inklusive 1, aber dass das zu nem minimalen Fehler führt hat er selbst schon geschrieben. Ich hab nur erst jetzt gecheckt wieso er meint dass das relevant sein könnte...

Ramanujan

Endlich einer, der mich versteht. Danke distributor.

Ganz nebenbei sollte die Einbettung immer funktionieren, auch wenn die Daten der Einbettungsmatrix nicht gleichverteilt sind.

Auf rand() werde ich demnächst wohl verzichten...

Edit:
C++11 bietet ja recht viele Zufallsgeneratoren an. Welchen soll ich benutzen?