In [1]:
import numpy as np
from scipy.stats import lomax
from scipy.optimize import fsolve
import pandas as pd
In [2]:
def lomax_moments_test(mean ,var): 
    # Solve the system of equations for the shape and scale parameters
    def equations(params):
        shape, scale = params
        eq1 = mean - scale / (shape - 1)
        eq2 = var - shape * scale**2 / ((shape - 1)**2 * (shape - 2))
        return [eq1, eq2]
    
    shape, scale = fsolve(equations, (3, 3))
    
    return shape, scale

def lomax_moments(data):
    mean = np.mean(data)
    var = np.var(data)
    
    # Solve the system of equations for the shape and scale parameters
    def equations(params):
        shape, scale = params
        eq1 = mean - scale / (shape - 1)
        eq2 = var - shape * scale**2 / ((shape - 1)**2 * (shape - 2))
        return [eq1, eq2]
    
    shape, scale = fsolve(equations, (3, 3))
    
    return shape, scale
In [3]:
c= 1.3
scale = 10

scale / (c - 1) #mean = 
c* scale**2 / ((c- 1)**2 * (c - 2)) #var = 

sample_means = []
sample_variances = []

for _ in range(50000):
    # Generate 1000 samples from the Lomax distribution
    samples = lomax.rvs(c, scale=scale, size=50)
    sample_mean = np.mean(samples)
    sample_means.append(sample_mean )    
    sample_variance = np.var(samples)
    sample_variances.append(sample_variance)    

print("mean of simulated means = " + str(np.mean(sample_means)))
print("median of simulated means = " + str(np.median(sample_means)))

print("mean of simulated variances = " + str(np.mean(sample_variances)))
print("median of simulated variances = " + str(np.median(sample_variances)))

mean of simulated means = 35.32232097822492
median of simulated means = 22.199454165491623
mean of simulated variances = 20133683.651214894
median of simulated variances = 1999.904985311406
In [ ]: