Zufällige Auswahl mit Speichenrad-Verfahren: Qualität des Codes



  • Hi,

    ich habe mal versucht das Speichenrad-Verfahren umzusetzen, weil ich gesehen habe, dass es nicht für jede Sprache eine Implementierung gibt.
    Das Speichenradverfahren hat eine Laufzeit von O(n) und ist sehr stabil und hat eine gute Diversität bei gleichverteilten Wahrscheinlichkeiten.
    Hier kann man auf S. 41 mehr dazu lesen:

    http://www-home.htwg-konstanz.de/~bittel/msi_robo/Vorlesung/08_MonteCarloLokalisierung.pdf

    Hier ist meine Umsetzung:

    import java.util.*;
    
    public class RandomChoice {
        public static List<Integer> choice(final Vector<Double> weights) {
            // Save accumulated sum of the weights in order to calculate the intervals
            // between accumulated weights
            final double accumulatedSum[] = new double[weights.size()];
    
            accumulatedSum[0] = weights.get(0);
    
            for (int i = 1; i < weights.size(); ++i) {
                accumulatedSum[i] = accumulatedSum[i-1] + weights.get(i);
            }
    
            final double sum = accumulatedSum[weights.size() - 1];
            final double spokeGap = sum / weights.size();
    
            // nextDouble() is an element of [0.0, 1.0]
            final double spin = new Random().nextDouble() * spokeGap;
    
            final List<Integer> choices = new ArrayList<>(weights.size());
            int i = 0;
            double intervalBegin = 0.0;
            double intervalEnd = accumulatedSum[i];
    
            double currentSpoke = spin;
    
            for (int j = 0; j < weights.size(); ++j) {
                currentSpoke += spokeGap;
                currentSpoke %= sum;
    
                while (intervalEnd <= currentSpoke || intervalBegin > currentSpoke) {
                    ++i;
    
                    i %= weights.size();
    
                    intervalBegin = accumulatedSum[i] - weights.get(i);
                    intervalEnd = accumulatedSum[i];
                }
    
                choices.add(i);
    
                if (choices.size() == weights.size()) break;
            }
    
            return choices;
        }
    
        public static void main(String args[]) {
            Vector<Double> test = new Vector<>();
    
            test.add(1.0);
            test.add(2.0);
            test.add(4.0);
            test.add(5.0);
            test.add(10.0);
    
            List<Integer> choices = RandomChoice.choice(test);
    
            for (int i : choices) {
                System.out.print(i + ", ");
            }
    
            System.out.println("\nLength: " + choices.size());
    
        }
    }
    

    Ich habe darauf geachtet so wenig wie möglich zu alloziieren. Ich musste an einer Stelle auch der Versuchung widerstehen zu sortieren, da das die Laufzeit auf O(n * log n) verschlechtern würde. Nun würde ich gerne eure Meinung dazu hören, bevor ich den Code auf andere Programmiersprachen portiere.

    L. G.,
    ShadowClone



  • Ich habe deine Implementierung noch nicht angeschaut, versuche erst den Algorithmus zu verstehen.
    Gibt es irgendeinen Grund weshalb die Anzahl der verschiedenen Werte/Gewichte (M) gleich der Anzahl der Ziehungen sein muss?
    Ich meine das ginge doch auch im verallgemeinerten Fall mit M Gewichten und N Ziehungen, ja? Mit d=W/N.



  • scrontch schrieb:

    Ich meine das ginge doch auch im verallgemeinerten Fall mit M Gewichten und N Ziehungen, ja? Mit d=W/N.

    Das stimmt. Da spricht nichts dagegen. Das wäre noch allgemeiner.



  • Ok, der Grund, weshalb ich nichts über das Speichenrad-Verfahren im Netzt gefunden habe ist, weil mein Prof. den Low Variance Resampling Algorithmus als Speichenrad veranschaulicht hat.
    Wenn ich wieder Zeit habe, werde ich mal andere Implementierungen anschauen, falls ich welche finde. – Eigentlich ist gerade Klausurzeit. 😃



  • Ich vermute fast, daß der Geist des Speichenradverfahrens der ist, daß auf dem Prozessor selbst die Zahlen modulo 2^(eingebaute Zahlenbreite) eine Speichenrad sind.

    W sei 2^32=4294967296.

    Oder W sei 2^8=256.

    Man müßte die Gewichte skalieren, damit die Gewichtesumme eine Zweierpozenz ist. Das bringt Ungenauigkeit, wenn W zu klein ist. Aber es ist granatenschnell.



  • FÜr das Skalieren muss ich Umrechnen und Casten. Zu dem kommt – wenn ich das richtig im Kopf habe – nochmals eine Schleife mit den Kosten von O(n) hinzu. Das ist zwar insgesamt immer noch O(n), ist aber konträr zu der ursprünglichen Absicht der Optimierung.
    Lohnt sich das?



  • Code nochmals leicht optimiert:

    public static List<Integer> choice(final Vector<Double> weights) {
        // Save accumulated sum of the weights in order to calculate the intervals
        // between accumulated weights
        final double accumulatedSum[] = new double[weights.size()];
    
        accumulatedSum[0] = weights.get(0);
    
        for (int i = 1; i < weights.size(); ++i) {
            accumulatedSum[i] = accumulatedSum[i-1] + weights.get(i);
        }
    
        final double sum = accumulatedSum[weights.size() - 1];
        final double spokeGap = sum / weights.size();
    
        // nextDouble() is an element of [0.0, 1.0]
        final double spin = new Random().nextDouble() * spokeGap;
    
        final List<Integer> choices = new ArrayList<>(weights.size());
        int i = 0;
    
        double currentSpoke = spin;
    
        for (int j = 0; j < weights.size(); ++j) {
            currentSpoke += spokeGap;
            currentSpoke %= sum;
    
            // accumulatedSum[i] is the interval end
            // (accumulatedSum[i] - weight.get(i)) is the interval begin
            // Following condition must be true: interval begin <= currentSpoke < interval end
            // If the condition is not true, then increment i in order to get the next greater intervals.
            while (accumulatedSum[i] <= currentSpoke || (accumulatedSum[i] - weights.get(i)) > currentSpoke) {
                ++i;
                i %= weights.size();
            }
    
            choices.add(i);
        }
    
        return choices;
    }
    


  • Der Plan von Seite 41 ist perfekt. Sorry für mein Posting vorhin. Ich hatte es nicht verstanden.

    Ich bin noch nicht sicher, ob ich Deinen Code verstehe. Ich habs mal ganz anders angedacht. Bestimmt falsch, aber vielleicht hilft es.

    List resampling(List pm){
        List result;
    
        //Aus der Partikelmenge {x[1], x[2], .., x[M]}
        //mit Gewichten w1,..., wM sollen M Partikeln
        //zufällig entsprechend ihrem Gewicht
        //gezogen werden.
        ///Ich mag annehmen, daß in pm Objekte mit .x und .w sind. 
    
        //Ordne Gewichte w1,..., wM auf einer
        //kreisförmigen Skala an. 
        //Es sei W = w1+...+wM die Gewichtsumme.
        ///Ich ordne sie nur auf dem Papier an. 
        ///Aber dei Gewichtesumme W brauch ich. 
        double W=0;
        foreach(i in x)
            W+=i.w;
    
        //Bilde ein Speichenrad mit M Speichen
        //und einem Speichenabstand von d = W/M.
        ///Auch das baue ich nur auf dem Papier. 
    
        //Drehe 1. Speiche auf einen zufälligen Wert
        //r ∈ [0, d]. 
        ///Ok, das mache ich mal…
        double r=Math.Rand()*W/pm.getSize();
    
        //Wähle jeden Partikel x[i] aus, auf dessen
        //Gewicht w[i] eine Speiche zeigt.
        //Zeigen mehrere Speichen auf ein Partikel,
        //dann wird der Partikel mehrfach ausgewählt. 
        int opferIndex=0;
        double gepacktGewicht=0;
        for(double stichGewicht=r;stichGewicht<M;stichGewicht+=d){
            while(pm[opferIndex].w<stichGewicht){
                ++opferIndex;
            }
            result.add(pm[opferIndex])
        }
        return result; 
    }
    

    Uups, hab ja gar nicht weniger Code als Du.



  • Ich werde mir da nochmal Gedanken machen. Vll. erst in eine Woche. Ich sollte eigentlich auf meine letzte Klausur lernen...



  • /// Chooses n samples by their weights. The greater their weights the more likely they get chosen.
    ///
    /// @invariant sum of weights must not overflow.
    /// @param samples The to be selected samples
    /// @param weights Weights that get chosen by their weight/probability. One weight can be greater 1.
    /// @param n Number of randomly chosen indices by weight.
    /// @return randomly selected samples by their weights
    pub fn random_choice<'a, T>(samples: &'a Vec<T>, weights: &Vec<f64>, n: usize) -> Vec<&'a T>{
        // TODO Check, if weight.len() > 0 and n > 0
    
        let sum:f64 = weights.iter().fold(0.0, |acc, &i| acc + i);
        let spoke_gap: f64 = sum / n as f64;
    
        // next_f64() ∈ [0.0, 1.0)
        let spin = thread_rng().next_f64() * spoke_gap;
    
        let mut i: usize = 0;
        let mut accumulated_weights = weights[0];
        let mut choices: Vec<&T> = Vec::with_capacity(n);
        let mut current_spoke: f64 = spin;
    
        for _ in 0 .. n {
            while accumulated_weights < current_spoke {
                i += 1;
                accumulated_weights += weights[i];
            }
            choices.push(&samples[i]);
            current_spoke += spoke_gap;
        }
    
        choices
    }
    
    pub fn random_choice_in_place<T: Copy>(samples: &mut Vec<T>, weights: &Vec<f64>) {
        // TODO Check, if weight.len() > 0 and n > 0
        let sum:f64 = weights.iter().fold(0.0, |acc, &i| acc + i);
        let n: usize = weights.len();
        let spoke_gap: f64 = sum / n as f64;
    
        // next_f64() ∈ [0.0, 1.0)
        let spin = thread_rng().next_f64() * spoke_gap;
    
        let mut j: usize = 0;
        let mut accumulated_weights = weights[0];
        let mut current_spoke: f64 = spin;
    
        for i in 0 .. n {
            while accumulated_weights < current_spoke {
                j += 1;
                accumulated_weights += weights[j];
            }
            samples[i] = samples[j];
            current_spoke += spoke_gap;
        }
    }
    

    Ich glaube, besser geht es nicht. Ich habe sogar eine in place Variante in der es 0 Heap-Alloziierungen gibt. 🕶


Log in to reply