Γλώσσα Προγραμματισμού ΙΙ - Εαρινό εξάμηνο 2015

Διάλεξη 26ης Μαρτίου 2015

Συνεχίζουμε τη συζήτηση για την βιβλιοθήκη NumPy με περισσότερες λεπτομέρεις για τις συναρτήσεις της NumPy που παράγουν τυχαίους αριθμούς.

Τυχαίοι αριθμοί

Αρκετή από την αρχική δουλειά στην θεωρία των πιθανοτήτων αφιερώθηκε στην ανάλυση τυχερών παιχνιδιών, κυρίως παιχνιδιών με ζάρια. Το ενδιαφέρoν ενός εκ των πρωτοπόρων της θεωρίας των πιθανοτήτων, του Γάλλου μαθηματικού Blaise Pascal, κινήθηκε από την ερώτηση ενός φίλου κατά πόσο θα ήταν επικερδές να στοιχηματίσει ότι δεν θα φέρει εξάρες αν ρίξει ένας ζεύγος από ζάρια 24 φορές. Σήμερα γνωρίζουμε ότι η πιθανότητα να μην εμφανιστούν εξάρες είναι (35/36)²⁴, δηλαδή περίπου 0.51. Μπορούμε να γράψουμε ένα πικρό πρόγραμμα για να προσομειώσουμε το συγκεκριμένο παιχνίδι και να προσεγγίσουμε πειραματικά την παραπάνω πιθανότητα.


def rollDie():
	return np.random.random_integers(1,6)

def playGame(numTrials):
	yes = 0.0
	for i in range(numTrials):
		for j in range(24):
			d1 = rollDie()
			d2 = rollDie()
			if d1 == 6 and d2 == 6:
				yes += 1
				break
	print 'Probability of losing = ' + str(1.0 - yes / numTrials)

Αν καλέσουμε τη συνάρτηση playGame με όρισμα, 100000, 200000, 400000 και 800000 παίρνουμε, αντίστοιχα, τα αποτελέσματα


Probability of losing = 0.50897
Probability of losing = 0.50998
Probability of losing = 0.508605
Probability of losing = 0.5088375

ενώ η τιμή της έκφρασης (35/36)²⁴ είναι κατά προσέγγιση 0.5085961238690966. Δοκιμάστε και εσείς να προσομοιώσετε το παραπάνω παιχνίδι. Σημειώνουμε ότι η κλήση random_integers(a,b) επιστρέφει ένα τυχαίο αριθμό στο διάστημα [a,b]. Άλλες χρήσιμες συναρτήσεις από το πακέτο random της βιβλιοθήκης NumPy είναι:

rand(d0,d1, ..., dn): Επιστρέφει τυχαίους αριθμούς (στο διάστημα [0, 1)) με το δεδομένο σχήμα (d0, d1, ..., dn). Για παράδειγμα, η κλήση np.random.rand(2,2) επιστρέφει ένα 2 x 2 πίνακα με τυχαία στοιχεία.
randn(d0,d1, ..., dn): Επιστρέφει ένα δείγμα σε σχήμα (d0, d1, ..., dn) από την τυπική κανονική κατανομή.
randint(low, high=None, size=None): Επιστρέφει τυχαίους ακέραιους στο διάστημα [low, high). Αν η παράμετρος high δωθεί ως None τότε οι αριθμοί επιλέγονται από το διάστημα [0, low). H τιμή της παραμέτρου size είναι είτε ένας ακέραιος είτε ένα tuple ακεραίων. Για παράδειγμα, η κλήση random.randint(2, size=10) επιστρέφει 10 τυχαίος ακέραιους στο διάστημα [0, 2) ενώ η κληση random.randint(5, size=(2, 4)) επιστρέφει ένα 2 x 4 πίνακα τα στοιχεία του οποίου είναι τυχαίοι ακέραιοι στο διάστημα [0, 5).
random(size=None): Επιστρέφει τυχαίους αριθμούς κινητής υποδιαστολής στο διάστημα [0.0, 1.0). Το όρισμα size μπορεί να είναι είτε ένας θετικός ακέραιος είτε ένα tuple. Σε πρίπτωση που το όρισμα size δωθεί, η συνάρτηση επιστρέφει έναν μόνο τυχαίο αριθμό.

Ως ένα τελευταίο παράδειγμα, θα υπολογίσουμε προσεγγιστικά την πιθανότητα να έχουμε το λεγόμενο "χρώμα", στα αγγλικά flush, δηλαδή πέντε φύλλα του ιδίου χρώματος, σε ένα παιχνίδι poker. Μια τράπουλα αποτελείται από 52 χαρτιά χωρισμένα σε τέσσερα χρώματα, τις κούπες (hearts), τα μπαστούνια (clubs), τα καρώ (diamonds) και σπαθιά (spades). Μπορούμε επομένως να φτιάξουμε εύκολα μια τράπουλα με την παρακάτω συνάρτηση:


def make_deck():
	suits = ['C', 'D', 'H', 'S']
	ranks = ['2', '3', '4', '5', '6', '7', '8', '9', '10', 'J', 'Q', 'K', 'A']
	deck = [s+r for s in suits for r in ranks]
	np.random.shuffle(deck)                        # ανακάτεμα της τράπουλας
	return deck

Αν θέλουμε να μοιράσουμε πέντε φύλλα από την τράπουλα, αρκεί να δώσουμε σε κάθε παίκτη πέντε φύλλα, για παρέδειγμα με την εντολή hand = deck[:5]. Θα χρειαστούμε επίσης και ένα τρόπο να διαπιστώσουμε αν τα χαρτιά τα οποία πήραμε είναι όλα του ιδίου χρώματος:


def same_suit(hand):
	suits = [card[0] for card in hand]
	count = {}
	for s in suits:
		n = suits.count(s)
		if n > 1:
			count[s] = n
	return count

Για να υπολογίσουμε τώρα την πιθανότητα να έχουμε "χρώμα" μπορούμε να ανακατέψουμε την τράπουλα, να μοιράσουμε πέντε φύλλα και να μετρήσουμε πόσα χέρια περιέχουν "χρώμα". Ο λόγος αυτού του αριθμού πρός τον συνολικό αριθμό που μοιράσαμε τα πέντε φύλλα θα είναι κοντά στην πραγματική πιθανότητα. Δεν είναι δύσκολο να υπολογίσει κανείς ότι η πιθανότητα χρώματος είναι ακριβώς 5148/2598960, δηλαδή περίπου 0.0019807923169267707.


# Shuffle the deck N times

N = 900000

flush = 0

for i in range(N):
	deck = make_deck()
	hand = deck[:5]

#	Count hands with five cards of the same suit (flush)

	count = same_suit(hand)
	if 5 in count.values():
		flush += 1

# Approximate probabilities

print 'Probability of a flush:', flush/float(N), 'Exact probability:', 5148/2598960.0

Πολυώνυμα

Η βιβλιοθήκη NumPy παρέχει το πακέτο numpy.polynomial για την αναπαράσταση πολυωνύμων, όπως και πλήθος συναρτήσεων που δρούν πάνω σε αυτά. Η αναπαράσταη ενός πολυωνύμου γίνεται με την αναγραφή των συντελεστών του από το σταθερό όρο προς τον μεγιστοβάθμιο όρο σε μια λίστα:


from numpy.polynomial import Polynomial as P
>>> p = P([1,2,3])
>>> p
Polynomial([ 1.,  2.,  3.], [-1.,  1.], [-1.,  1.])
>>> p.coef
array([ 1.,  2.,  3.])

Σημειώνουμε ότι η πρώτη γραμμή κώδικα εισάγει το όνομα P ως συντόμευση του ονόματος Polynomial. Διαφορετικά, η 2η γραμμή που κατασκευάζει το πολυώνυμο p θα έπρεπε να γραφεί ως p = np.polynomial.Polynomial([1, 2, 3]). Το πολυώνυμο στο οποίο αναφερόμαστε εδώ είναι φυσικά το $p(x) = 1 + 2x + 3x^2$. Προς το παρόν δεν θα ασχοληθούμε με τις λίστες [-1., 1.] και [-1., 1.] που επιστρέφει η συνάρτηση Polynomial. Πράξεις μεταξύ πολυωνύμων μπορούν να γίνουν με τον αναμενόμενο τρόπο:


>>> p + p
Polynomial([ 2.,  4.,  6.], [-1.,  1.], [-1.,  1.])
>>> p - p
Polynomial([ 0.], [-1.,  1.], [-1.,  1.])
>>> p * p
Polynomial([  1.,   4.,  10.,  12.,   9.], [-1.,  1.], [-1.,  1.])
>>> p**2
Polynomial([  1.,   4.,  10.,  12.,   9.], [-1.,  1.], [-1.,  1.])
>>> p / 2
Polynomial([ 0.5,  1. ,  1.5], [-1.,  1.], [-1.,  1.])
>>> p // P([-1, 1])
Polynomial([ 5.,  3.], [-1.,  1.], [-1.,  1.])
>>> p % P([-1, 1])
Polynomial([ 6.], [-1.,  1.], [-1.,  1.])

Σημειώνουμε ότι για την διαίρεση πολυωνύμων θα μπορούσαμε να είχαμε γράψει:


>>> quo, rem = divmod(p, P([-1, 1]))
>>> quo
Polynomial([ 5.,  3.], [-1.,  1.], [-1.,  1.])
>>> rem
Polynomial([ 6.], [-1.,  1.], [-1.,  1.])

Για να υπολογίσουμε την τιμή ενός πολυωνύμου σε ένα ή περισσότερα σημεία αρκεί να δώσουμε το ή τα σημεία μέσα σε παρενθέσεις μετά το όνομα του πολυωνύμου. Οι ρίζες του πολυωνύμου μπορούν να βρεθούν με τη μέθοδο roots():


>>> x = np.arange(5)
>>> p(x)
array([  1.,   6.,  17.,  34.,  57.])
>>> p.roots()
array([-0.33333333-0.47140452j, -0.33333333+0.47140452j])

Τόσο ολοκλήρωση όσο και παραγώγιση πολυωνύμων υποστηρίζεται από συναρτήσεις του πακέτου numpy.polynomial:


>>> p = P([2, 6])
>>> p.integ()
Polynomial([ 0.,  2.,  3.], [-1.,  1.], [-1.,  1.])
>>> p.integ(2)
Polynomial([ 0.,  0.,  1.,  1.], [-1.,  1.], [-1.,  1.])

Παρατηρήστε ότι στο πρώτο παράδειγμα η σταθερά της ολοκλήρωσης είναι μηδέν ενώ στο δεύτερο παράδειγμα το πολυώνυμο $p(x) = 2 + 6x^2$ ολοκληρώνεται δυό φορές (με τη σταθερά ολοκλήρωσης να εξακολουθεί να είναι μηδέν). Η μέθοδος deriv() για την παραγώγιση ενός πολυωνύμου λειτουργεί με ανάλογο τρόπο:


>>> p = P([1, 2, 3])
>>> p.deriv(1)
Polynomial([ 2.,  6.], [-1.,  1.], [-1.,  1.])
>>> p.deriv(2)
Polynomial([ 6.], [-1.,  1.], [-1.,  1.])