Tarea 2 estocásticos

import math import plotly.express as px import numpy as np import pandas as pd from decimal import Decimal, InvalidOperation

precalced = dict() def dynamic_factorial(val): key = ('factorial', val) try: if key not in precalced: precalced[key] = math.factorial(val.as_integer_ratio()[0]) return precalced[key] except InvalidOperation: print(val) def poisson(lmbda, t, k): key = ('poisson', lmbda, t, k) if key not in precalced: mult = Decimal(lmbda * t) true_k = Decimal(k) exp = (-mult).exp() multiplication = Decimal(1) if not mult and not true_k else (mult**true_k) nominator = multiplication*exp denominator = dynamic_factorial(true_k) precalced[key] = nominator / denominator return precalced[key] def poisson_sum(lmbda, t, start=0, end=1): key = ('poisson_sum', lmbda, t, start, end) if key not in precalced: precalced[key] = sum(poisson(lmbda, t, k) for k in range(start, end)) return precalced[key] def probability_func(n, t, lmbda, p, q, r): key = ('probability_func', n, t, lmbda, p,q, r) if key not in precalced: precalced[key] = sum( poisson_sum(lmbda*q,t,start=0, end=i) * poisson_sum(lmbda*r,t,start=0, end=i) * poisson(lmbda * p, t, i) for i in range(1,n+1) ) return precalced[key]

data_points = [] time_points = 30 n_points = list(range(1, 150, 20)) hyper_lambda = 16 hyper_p = 0.55 hyper_q = 0.30 hyper_r = 0.15 for t in range(time_points): for n in n_points: data_points.append({"t": t, "N": n, "Probabilidad": probability_func(n, t, hyper_lambda, hyper_p, hyper_q, hyper_r)}) df = pd.DataFrame(data_points) fig = px.line(df, x='t', y='Probabilidad', color='N', title=r"$\text{Grafico de } P(\mathbb{N}_1(t) \geq \mathbb{N}_j(t) ~\forall j \in \{2,3\})\text{ para distintos } n \text{ tal que } n \rightarrow \infty $", line_shape="spline") fig

data_points = [] time_points = 30 n_points = list(range(1, 150, 20)) hyper_lambda = 16 hyper_p = 0.55 hyper_q = 0.30 hyper_r = 0.15 for t in range(time_points): for n in n_points: data_points.append({"t": t, "N": n, "Probabilidad": probability_func(n, t, hyper_lambda, hyper_q, hyper_p, hyper_r)}) df = pd.DataFrame(data_points) fig = px.line(df, x='t', y='Probabilidad', color='N', title=r"$\text{Grafico de } P(\mathbb{N}_2(t) \geq \mathbb{N}_j(t) ~\forall j \in \{1,3\})\text{ para distintos } n \text{ tal que } n \rightarrow \infty $", line_shape="spline") fig

data_points = [] time_points = 30 n_points = list(range(1, 150, 20)) hyper_lambda = 16 hyper_p = 0.55 hyper_q = 0.30 hyper_r = 0.15 for t in range(time_points): for n in n_points: data_points.append({"t": t, "N": n, "Probabilidad": probability_func(n, t, hyper_lambda, hyper_r, hyper_p, hyper_q)}) df = pd.DataFrame(data_points) fig = px.line(df, x='t', y='Probabilidad', color='N', title=r"$\text{Grafico de } P(\mathbb{N}_3(t) \geq \mathbb{N}_j(t) ~\forall j \in \{1,2\})\text{ para distintos } n \text{ tal que } n \rightarrow \infty $", line_shape="spline") fig

def question_b(lmbda, p, q, r, k, t): return ( poisson(lmbda * p, t, k) + poisson(lmbda * q, t, k) + poisson(lmbda * r, t, k) )

data_points = [] time_points = 100 n_points = list(range(1, 150, 20)) hyper_lambda = 16 hyper_p = 0.55 hyper_q = 0.30 hyper_r = 0.15 for t in range(time_points): for n in n_points: data_points.append({"t": t, "k": n, "Probabilidad": question_b(hyper_lambda, hyper_p, hyper_q, hyper_r, n, t)}) df = pd.DataFrame(data_points) print(data_points) fig = px.line(df, x='t', y='Probabilidad', color='k', title=r"$\text{Grafico de } \mathbb{N}(t) = \mathbb{N}_1(t) + \mathbb{N}_2(t) + \mathbb{N}_3(t) \text{ para distintos } k $", line_shape="spline") fig

import random import statistics random.seed(16624823)

staying = 0.45 moving = 0.35 falling = 0.2

def simulate(prob_staying, prob_moving, prob_falling, boxes, iterations, max_seconds=math.inf): outcomes = [] for _ in range(iterations): on_tempered = True steps_taken = 0 seconds_passed = 0 while on_tempered and steps_taken < boxes and seconds_passed < max_seconds: choice, *_ = random.choices(["stay", "move", "fall"],[prob_staying, prob_moving, prob_falling], k=1) if choice == "move": steps_taken += 1 elif choice == "fall": on_tempered = False steps_taken += 1 seconds_passed += 1 outcomes.append(steps_taken) return outcomes simulation = simulate(staying, moving, falling, 25, 10000) print("Promedio de casilla", statistics.mean(simulation)) print("Proporción personas que llegan a una casilla k>=13", (len(list(filter(lambda steps: steps >= 13, simulation))) / len(simulation))*100) simulation2 = simulate(staying, moving, falling, 25, 10000, max_seconds=60) print("Casilla mas probable despues de 60 segundos transcurridos:", sorted(simulation2, key=lambda steps: simulation2.count(steps), reverse=True)[0]) sorted(simulation2, key=lambda steps: simulation2.count(steps), reverse=True)

def simulate_dynamic(prob_staying, prob_moving, prob_falling, boxes, iterations, lmbda, max_seconds=math.inf, hard_break_step=None): outcomes = [] for _ in range(iterations): on_tempered = True steps_taken = 0 seconds_passed = 0 if not hard_break_step: hard_break_step = boxes while on_tempered and steps_taken < hard_break_step and seconds_passed < max_seconds: lambda_k = lmbda * steps_taken dyn_staying = prob_staying + lambda_k dyn_moving = prob_moving - lambda_k/2 dyn_falling = prob_falling - lambda_k/2 choice, *_ = random.choices(["stay", "move", "fall"],[dyn_staying, dyn_moving, dyn_falling], k=1) if choice == "move": steps_taken += 1 elif choice == "fall": on_tempered = False steps_taken += 1 seconds_passed += 1 outcomes.append(steps_taken) return outcomes

from functools import partial def to_prob_of_reaching(simulation, k): hits = sum(filter(lambda outcome: outcome >= k, simulation)) total = len(simulation) return hits/total values = [] partial_sim = partial(simulate_dynamic, staying, moving, falling, 25, 10000, hard_break_step=13) for pre_lmbda in range(0, int(0.1/0.01)): calculated_lmbda = pre_lmbda*0.01 for max_seconds in range(0, 300): values.append([calculated_lmbda, max_seconds, to_prob_of_reaching(partial_sim(calculated_lmbda, max_seconds=max_seconds), k=13)])

df = pd.DataFrame(np.array(values), columns=['Lambda','t', 'Probabilidad']) fig = px.line(df, x='t', y='Probabilidad',color='Lambda', title=r"$\text{Grafico de probabilidad de que lleguen a la casilla } k=13$", line_shape="linear") fig