In [2]:
%pylab
%matplotlib inline
import scipy.stats as st
In de bovenstaande afbeelding zie je een mier in zijn natuurlijke leefomgeving: een rij hokjes. De mier zal gaan wandelen, op zoek naar voedsel of andere eerste levensbehoeften. We modelleren het wandelgedrag van de mier als volgt: iedere tijdseenheid (stap) loopt de mier met een kans van 50% 1 vakje naar links en zo niet, dan loopt hij 1 vakje naar rechts. Het doel voor de mier is om de suiker te bereiken in een van de uiterste vakjes.
Beantwoord de volgende vraag met een simulatie (experiment 10000 keer uitvoeren): hoeveel tijdseenheden (stapjes) duurt het gemiddeld tot de mier in het meest linkse vakje of het meest rechtse vakje aankomt?
In [3]:
n = 10000
steps_to_exit = []
for i in range(n):
x = 0
steps = 0
while -7 < x < 7:
x += np.random.choice([-1, 1]) # step left or right
steps += 1
steps_to_exit.append(steps)
print("Gemiddeld aantal stappen tot suiker: {:.3f}".format(mean(steps_to_exit)))
Het model in oefening 1 beschrijft een random walk in één dimensie. In deze oefening doe je een random walk in 2D.
Een man stapt om middernacht stomdronken uit een café midden in de stad. Hij wil naar huis, maar zijn richtingsgevoel en coördinatie laten hem enigszins in de steek: iedere stap die hij zet is in een willekeurige richting. Hij loopt steeds exact 1 stap in een willekeurige richting, noord, zuid, oost of west. Hoe groot is de kans dat de man in hooguit 10000 stappen terug komt bij het café?
In [9]:
n = 100
successes = 0
steps_used = []
for i in range(n):
x, y = 0, 0
steps = 0
for j in range(10000):
direction = np.random.choice(list("NESW"))
if direction == "N":
y += 1
elif direction == "E":
x += 1
elif direction == "S":
y -= 1
else:
x -= 1
steps += 1
if (x, y) == (0, 0):
# we zijn terug in het cafe (0, 0)
successes += 1
steps_used.append(steps)
break
else:
# steps_used.append(10000)
pass
print("Kans op terugkeren in café in maximaal 10000 stappen: {:.3f}".format(successes / n))
print(np.mean(steps_used))
print(len(steps_used))
In [10]:
n = 10000
successes = 0
for i in range(n):
sides = list(np.random.rand(3))
sides.sort()
if sides[2] < sum(sides[:2]):
successes += 1
print("Kans op mogelijkheid driehoek bij benadering {:.3f}".format(successes / n))
In [11]:
%pylab
%matplotlib inline
x = linspace(0, 1, 400)
y = sqrt(1 - x**4)
plt.fill_between(x, y)
plt.axis('equal')
plt.ylim([0, 1.1])
plt.show()
Met behulp van Monte Carlo simulatie kun je de integraal wel vrij eenvoudig benaderen. Je kunt een rechthoek kiezen waar je integraal binnen valt en dan voor willekeurige punten in de rechthoek bepalen of ze onder of boven de grafiek van de functie vallen.
Benader de gegeven integraal met een Monte Carlo simulatie.
In [18]:
n = 500000
under_graph = 0
for i in range(n):
x = np.random.rand()
y = np.random.rand()
fx = sqrt(1 - x**4)
if y < fx:
under_graph += 1
print("Schatting kans willekeurig punt onder grafiek = {:.3f}".format(under_graph / n))
print("Schatting integraal = {:.5f}".format((1 * 1) * (under_graph / n)))
In [ ]: