probabilidad - ¿Cómo se muestran al azar los datos lognormales en Python usando el CDF inverso y se especifican los percentiles de destino?

probabilidad y estadistica en python

Question

Estoy tratando de generar muestras aleatorias de una distribución lognormal en Python, la aplicación es para simular el tráfico de red. Me gustaría generar muestras tales que:

1. El resultado de muestra modal es 320 (~ 10 ^ 2.5)
2. El 80% de las muestras se encuentran dentro del rango de 100 a 1000 (10 ^ 2 a 10 ^ 3)

Mi estrategia es usar el CDF inverso (o la transformación Smirnov creo):

1. Use el PDF para una distribución normal centrada alrededor de 2.5 para calcular el PDF para 10 ^ x donde x ~ N (2.5, sigma).
2. Calcule el CDF para la distribución anterior.
3. Genere datos aleatorios uniformes a lo largo del intervalo de 0 a 1.
4. Use el CDF inverso para transformar los datos aleatorios uniformes en el rango requerido.

El problema es que cuando calculo el percentil 10 y 90 al final, tengo los números equivocados.

Aquí está mi código:

%matplotlib inline

import matplotlib
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.stats import norm

# find value of mu and sigma so that 80% of data lies within range 2 to 3
mu=2.505
sigma = 1/2.505
norm.ppf(0.1, loc=mu,scale=sigma),norm.ppf(0.9, loc=mu,scale=sigma)
# output: (1.9934025, 3.01659743)

# Generate normal distribution PDF
x = np.arange(16,128000, 16) # linearly spaced here, with extra range so that CDF is correctly scaled
x_log = np.log10(x)
mu=2.505
sigma = 1/2.505
y = norm.pdf(x_log,loc=mu,scale=sigma)
fig, ax = plt.subplots()
ax.plot(x_log, y, ''r-'', lw=5, alpha=0.6, label=''norm pdf'')

x2 = (10**x_log) # x2 should be linearly spaced, so that cumsum works (later)
fig, ax = plt.subplots()
ax.plot(x2, y, ''r-'', lw=5, alpha=0.6, label=''norm pdf'')
ax.set_xlim(0,2000)

# Calculate CDF
y_CDF = np.cumsum(y) / np.cumsum(y).max()
fig, ax = plt.subplots()
ax.plot(x2, y_CDF, ''r-'', lw=2, alpha=0.6, label=''norm pdf'')
ax.set_xlim(0,8000)

# Generate random uniform data
input = np.random.uniform(size=10000)

# Use CDF as lookup table
traffic = x2[np.abs(np.subtract.outer(y_CDF, input)).argmin(0)]

# Discard highs and lows
traffic = traffic[(traffic >= 32) & (traffic <= 8000)]

# Check percentiles
np.percentile(traffic,10),np.percentile(traffic,90)

Que produce la salida:

(223.99999999999997, 2480.0000000000009)

... y no el (100, 1000) que me gustaría ver. Cualquier consejo apreciado!

Answer 1

Primero, no estoy seguro acerca de Use the PDF for a normal distribution centred around 2.5 . Después de todo, log-normal es sobre logaritmo e base (también conocido como registro natural), lo que significa 320 = 10 ^2.5 = e ^5.77 .

En segundo lugar, abordaría el problema de una manera diferente. Necesita m y s para muestrear desde Log-Normal .

Si miras el artículo de wiki anterior, podrías ver que se trata de una distribución de dos parámetros. Y tienes exactamente dos condiciones:

Mode = exp(m - s*s) = 320
80% samples in [100,1000] => CDF(1000,m,s) - CDF(100,m,s) = 0.8

donde CDF se expresa mediante la función de error (que es una función bastante común que se encuentra en cualquier biblioteca)

Entonces dos ecuaciones no lineales para dos parámetros. Resuélvalos, encuentre m y s y póngalo en cualquier muestreo log-normal estándar

Answer 2

El enfoque de Severin es mucho más delicado que mi intento original de usar la transformación de Smirnov. Este es el código que funcionó para mí (usando fsolve para encontrar s, aunque es bastante trivial hacerlo manualmente):

# Find lognormal distribution, with mode at 320 and 80% of probability mass between 100 and 1000
# Use fsolve to find the roots of the non-linear equation

%matplotlib inline

import matplotlib
import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import fsolve
from scipy.stats import lognorm
import math

target_modal_value = 320

# Define function to find roots of
def equation(s):

    # From Wikipedia: Mode = exp(m - s*s) = 320
    m = math.log(target_modal_value) + s**2

    # Get probability mass from CDF at 100 and 1000, should equal to 0.8.
    # Rearange equation so that =0, to find root (value of s)
    return (lognorm.cdf(1000,s=s, scale=math.exp(m)) - lognorm.cdf(100,s=s, scale=math.exp(m)) -0.8)

# Solve non-linear equation to find s
s_initial_guess = 1
s =  fsolve(equation, s_initial_guess)

# From s, find m
m = math.log(target_modal_value) + s**2
print(''m=''+str(m)+'', s=''+str(s)) #(m,s))

# Plot
x = np.arange(0,2000,1)
y = lognorm.pdf(x,s=s, scale=math.exp(m))
fig, ax = plt.subplots()
ax.plot(x, y, ''r-'', lw=5, alpha=0.6, label=''norm pdf'')
plt.plot((100,100), (0,1), ''k--'')
plt.plot((320,320), (0,1), ''k-.'')
plt.plot((1000,1000), (0,1), ''k--'')
plt.ylim(0,0.0014)
plt.savefig(''lognormal_100_320_1000.png'')