<jupyter_start><jupyter_text>Métodos Montecarlo Fing 2022 - Entrega 3 **Autor**: <NAME>, marzo-abril 2022. **Email**: <EMAIL>, <EMAIL> Ejercicio 6.1se idealiza una montaña como un cono inscrito en una region cuadrada de lado 1 km. La base de la montaña es circular, con centro en (0.5, 0.5) y radio r = 0.4km, y la altura es H = 8km. La altura de cada punto (x, y) de la montaña está dada por la función:f (x, y) = H − H/r × sqrt((x − 0.5)2 + (y − 0.5)2), en la zona definida por el círculo, y 0 fuera del círculo.El volumen total de la montaña (en km cúbicos) puede verse como la integral de la función altura en la región.**Parte a:**Escribir un programa para calcular el volúmen por Monte Carlo. Realizar 10^6 replicaciones y estimar el valor de ζ y el error cometido (con nivel de confianza 0.95), utilizando como criterio la aproximación normal.<jupyter_code>import random
import math
import tabulate
import time
from IPython.core.display import HTML
random.seed()

import cm2c.fing.mmc.integral as mmci
import cm2c.fing.mmc.utils as mmcutils
reloj_ppal = mmcutils.timeit()
mmci.version()
# Validacion: integro f(x) = x**2 en (0,1)

import math

r = mmci.integracionMonteCarlo(lambda x: x[0]**2, 1, 10**5, lambda x: random.uniform(0,1))

print("El resultado debe aproximarse a 1/3: ", r[0])<jupyter_output>El resultado debe aproximarse a 1/3:  0.3332985227530386<jupyter_text>Ahora, para resolver el problema del volumen de la montaña defino una función en R2 que me devuelva la altura de la montaña abstracta:<jupyter_code>H = 8.0 # altura en km
r = 0.4 # radio de la base en km
n = 10**6 # cantidad de muestras para la parte (a)
delta = 0.05

import math

def Montana(x):
    """
    x es un vector de dos elementos
    devuelve la altura estimada
    """
    
    # calculo distancia al centro
    d = math.sqrt( (x[0]-0.5)**2 + (x[1]-0.5)**2 )
    
    if d > 0.4:
        return 0.0
    else:
        return H - (H/r)*d
## end def Montana

(estimZ, estimV, _, _) = mmci.integracionMonteCarlo(Montana, 2, n, mmcutils.sortearPuntoRN)

(icn0, icn1) = mmci.intConfianzaAproxNormal(estimZ, estimV, n, delta)

epsilon_est = estimZ-icn0

print(" ")
print("Volumen estimado por MMC {:.5f} km3".format(estimZ))
print("Varianza estimada : {:.5e}".format(estimV))
print("Intervalo de confianza para delta {} : ({:.5f}, {:.5f}) ".format(delta, icn0, icn1))
print("Error estimado: {:.5e}".format(epsilon_est))<jupyter_output>Volumen estimado por MMC 1.33848 km3
Varianza estimada : 3.55871e-06
Intervalo de confianza para delta 0.05 : (1.33478, 1.34218) 
Error estimado: 3.69738e-03<jupyter_text>**Parte b:**En base al valor estimado en la parte a, calcular el número dereplicaciones necesario para obtener un error absoluto menor a 10^-3(con nivel de confianza 0.95).- Paso 1: Realizar un número n' de pruebas preliminares y estimar la varianza.    En este caso ya tenemos una estimación para n' = 10^6<jupyter_code>print("Volumen estimado por MMC {:.5f} km3".format(estimZ))
print("Varianza estimada : {:.5e}".format(estimV))<jupyter_output>Volumen estimado por MMC 1.34109 km3
Varianza estimada : 3.56573e-06<jupyter_text>- Paso 2  Calcular el tamaño de la muestra requerido según la aproximación normal, el cual se define por:    ```npuntoN(epsilon, delta) = norm.ppf(1-delta/2)**2 * estimV / (epsilon**2)```<jupyter_code>from scipy.stats import norm

def npuntoN(delta_, epsilon_, estimV_, n_):
    return (norm.ppf(1-delta_/2)**2)*estimV_*n_/(epsilon_**2)

npuntoN_est = math.ceil ( npuntoN(0.05, 0.001, estimV, n) )

print(f'{estimV} , {n}')
print(f'{npuntoN_est:,}')<jupyter_output>3.56573492794628e-06 , 1000000
13,697,624<jupyter_text>- Paso 3  Repetir las simulaciones con N > npuntoN_estimado, asegurando que las semillas son diferentes (para garantizar la independencia de los experimentos)<jupyter_code># tryN = [ 14*10**6, 16*10**6, 20*10**6, 30*10**6]
tryN = [ 14*10**4, 16*10**4, 20*10**4, 30*10**4]

table1 = [ ['N', 'Volumen est. MMC', 'Varianza est.', 'Int. Confianza', 'Error est.', 'Tiempo ejecución'] ]

reloj = mmcutils.timeit()

for n in tryN:
    (estimZ, estimV, _, _) = mmci.integracionMonteCarloParalelo(Montana, 2, n, 8)
    (icn0, icn1) = mmci.intConfianzaAproxNormal(estimZ, estimV, n, delta)
    epsilon_est = estimZ-icn0
    table1.append([f'{n:,}', f'{estimZ:.5f}', f'{estimV:.5e}', f'({icn0:.5f}, {icn1:.5f})', f'{epsilon_est:.5e}', f'{reloj.lap():.3f}'])

tabulate.tabulate(table1, tablefmt='html')<jupyter_output><empty_output><jupyter_text>**Conclusión** La formula de estimación de cantidad de muestras nos da una referencia, en este caso unos 13.6x10^6 muestras. Probando con cantidades de muestras levemente mayores ya logramos bajar el error de 10^-3 que era lo pedido.   Ejercicio 6.2 Problema: Se desea estimar la integral de la función ```F5(X) = x1*x2^2*x3^3*x4^4*x5^5``` sobre elhipercubo J^m de dimensión m = 5.**Parte a:**Revisar los códigos preparados para el ejercicio 6.1, elegir unode ellos como punto de partida.Sobre esa base, modificarlo para realizar cálculo por Monte Carlo de la integral planteada en el ejercicio6.2. Realizar 10^6 replicaciones y estimar el valor de ζ. Calcular analíticamente el valor exacto de la integral.<jupyter_code># Definimos la función

def F5(x):
    """
    x es un vector en R^5
    """
    return x[0] * (x[1]**2) *  (x[2]**3) * (x[3]**4) * (x[4]**5)
# end def

# Calculo la integral con 10^6 replicas
n = 10**6

(estimZ, estimV, _, _) = mmci.integracionMonteCarlo(F5, 5, n)

HTML(f'<h4>estimZ: {estimZ:.7f} - estimV: {estimV:.5e}</h4>')<jupyter_output><empty_output><jupyter_text>Cálculo analítico de la integral de la función F5 en J^5:![Calculo analítico de la integral de F5](integral_analitica.gif "Cálculo analítico del valor de la integral en J^5 de F5") **Parte b:**En base al valor estimado en la parte a, calcular el número dereplicaciones necesario para obtener un error menor a 10−4 (con nivelde confianza 0.95).<jupyter_code>npuntoN_est = math.ceil ( npuntoN(0.05, 0.0001, estimV, n) )

print(f'Estimación de cantidad de muestras: {npuntoN_est:,}')<jupyter_output>Estimación de cantidad de muestras: 37,176<jupyter_text>**Parte c:**Decimos que un intervalo de confianza cubre el valor exactocuando este último pertenece al intervalo.Realizar L = 500 experimentos con semillas diferentes, cada unoconsistente en estimar por Monte Carlo con el nro. de replicaciones dela parte b el valor de la integral, así como intervalos de confianza denivel 0.9, 0.95 y 0.99.Para cada nivel de confianza, calcular el nivel decobertura empírico (en que porcentaje de los 500 experimentos elintervalo de confianza calculado cubrió el valor exacto).Discutir los resultados, comparando la cobertura empírica con laespecificada.<jupyter_code>table2 = [ ['L', 'n', '1-delta', '% Experimentos exitosos', 'Tiempo ejecución'] ]
L = 500
delta = 0.10
ZAnalitico = 1.0/720.0
npuntoN_est = math.ceil ( npuntoN(delta, 0.0001, estimV, n) )
reloj = mmcutils.timeit()

def Lexperimentos():
    cobertura = []
    for j in range(0,L):
        random.seed()
        (estimZ_, estimV_, _, _) = mmci.integracionMonteCarloParalelo(F5, 5, npuntoN_est, 8)
        (icn0, icn1) = mmci.intConfianzaAproxNormal(estimZ_, estimV_, npuntoN_est, delta)
        epsilon = ZAnalitico-icn0
        if ZAnalitico>=icn0 and ZAnalitico<=icn1:
            cobertura.append(1)
        else:
            cobertura.append(0)
    # end for
    # La cantidad de unos en el array cobertura es la cantidad de experimentos donde el valor analitico esta contenido
    # en el intervalo de confianza
    exitos = cobertura.count(1)

    table2.append( [L, f'{npuntoN_est:,}', f'{1-delta:.2f}', f'{exitos/L*100:.2f}%', f'{reloj.lap():.3f}' ] )
# end def

Lexperimentos()
print(reloj.lap())

delta = 0.05
npuntoN_est = math.ceil ( npuntoN(delta, 0.0001, estimV, n) )
Lexperimentos()
delta = 0.01
npuntoN_est = math.ceil ( npuntoN(delta, 0.0001, estimV, n) )
Lexperimentos()
tabulate.tabulate(table2, tablefmt='html')<jupyter_output><empty_output><jupyter_text>**Discusión de los resultados de cobertura**Me resulta interesante ver de que la cantidad de muestras que obtenemos de la formula basada en la normal partiendo de la estimación de la varianza obtenida para 10^6 muestras es significativamente más pequeña que 10^6 (del orden de 3.6 * 10^4).Cuando realizamos los L=500 experimentos partiendo de estas cantidades de muestras calculadas de esta forma (bastante más pequeñas) obtenemos porcentajes de cobertura que parecen bastante buenos.Al disminuir el tamaño del intervalo de confianza, naturalmente la cantidad de muestras crece. Es interesante ver también que el porcentaje de experimentos exitosos también crece.  Datos adicionales y referencias Información acerca del software y hardware utilizados**Software:** - Python 3.8.10 corriendo en Windows WSL2 (Windows Subsystem for Linux) - Jupyter Notebook  **Librerias:**  - scipy norm - pathos multiprocessing (para parelilizar ejecuciones) **Hardware:** - PC Windows 11, con WSL2 - CPU Intel Core i5 10400F (6 cores) - 16 GB de RAM<jupyter_code>print(f"%% FIN - tiempo total de ejecución {reloj_ppal.lap():.3f}s")<jupyter_output>%% FIN - tiempo total de ejecución 38.740s
{'fraction_non_alphanumeric': 0.11801919720767888, 'fraction_numerical': 0.049738219895287955, 'mean_word_length': 4.649414664202095, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '106094', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 3744, 'n_tokens_neox': 3464, 'n_words': 1072, 'path': 'e03_integrales_mmc.ipynb', 'repo_name': 'carlosm3011/fing-montercarlo-2022', 'size': '9203'}