Question: I am trying to estimate quantiles for some snow data using the log pearson type 3 distribution in Python and comparing with R. I do this by reading in the data, log transforming it, fitting Pearson type 3, estimating quantiles, then transforming back from log space. In python: import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats import lmoments3 as lm import lmoments3.distr as ld data=np.array([[1079], [ 889], [ 996], [1476], [1567], [ 897], [ 991], [1222], [1372], [1450], [1077], [1354], [1029], [1699], [ 965], [1133], [1951], [1621], [1069], [ 930], [1039], [1839]]) return_periods = np.array([2,3,5,10,20,50,100,200,1000]) log_data = np.log(data) params = stats.pearson3.fit(log_data) #Max likelihood estimation method quantiles = np.exp(stats.pearson3.ppf(1 - 1 / return_periods, *params)) paramsmm=ld.pe3.lmom_fit(log_data) #lmoments estimation method paramsmm2=(paramsmm["skew"], paramsmm['loc'], paramsmm['scale'][0]) quantilesmm = np.exp(ld.pe3.ppf(1 - 1 / return_periods, *paramsmm2)) print(quantiles) print(quantilesmm)  in R: library(lmom) library(lmomco) library(FAdist) swe_data <- c(1079, 889, 996, 1476, 1567, 897, 991, 1222, 1372, 1450, 1077, 1354, 1029, 1699, 965, 1133, 1951, 1621, 1069, 930, 1039, 1839) return_periods <- c(2, 3, 5, 10, 20, 50, 100, 200, 1000) exceedance_probabilities <- 1 / return_periods # P = 1/T nonexceedance_probabilities <- 1 - exceedance_probabilities # P_nonexceedance = 1 - P_exceedance log_swe <- log(swe_data) loglmoments <- lmom.ub(log_swe_data) fit_LP3 <- parpe3(loglmoments) #pars estimated using lmoments LP3_est=exp(quape3(nonexceedance_probabilities, fit_LP3)) print(LP3_est)  The quantiles estimated are the following: MLE/scipy stats: params=(2.0246357656236125, 7.10812763271725, 0.32194785836668816) #skew, loc scale quantiles=[1105.86050592 1259.46110488 1484.67412496 1857.18767881 2324.18036925 3127.68767927 3916.2007443 4904.15011095 8271.24322709]  Lmoments/python: params=(-2.2194418726874434, 7.1069179914286424, 0.07535915093401913) #skew, loc scale quantiles=[1251.30865382 1276.35189073 1291.29995882 1300.06624583 1303.59129662 1305.31725745 1305.78638777 1305.98555852 1306.11275037]  Lmoments/R: params= (7.1069180 0.2566677 0.9365001) #mu, sigma, gamma quantiles=[1173.116 1313.849 1485.109 1721.131 1969.817 2326.812 2623.112 2945.728 3814.692]  I would expect the latter two methods, both using lmoments, to produce the same result. Based on comparisons with other distributions, it seems like R is giving the most realistic result. Any explanation for the large differences? How might I get a similar result in Python?

Answer: Okey to my understanding basically, The main difference between the Python and R results stems from the estimation methods used: Python's scipy.stats uses Maximum Likelihood Estimation (MLE) by default. R's lmom package uses L-moments estimation. These methods can produce different parameter estimates, especially for small sample sizes or highly skewed distributions like the Log-Pearson Type III.To get a similar or closer result you can use the lmoments3 library in Python, which implements the L-moments method. Ensure correct parameter extraction from the lmoments3 results. To get them to agree, use the lmoments3 library in Python – it uses the same L-moments method as R. Just double-check how you're pulling out the parameters, especially the 'scale' value, from the results. Here's how you can get python result much closer to R results. import numpy as np import scipy.stats as stats import lmoments3 as lm import lmoments3.distr as ld # Input data data = np.array([1079, 889, 996, 1476, 1567, 897, 991, 1222, 1372, 1450, 1077, 1354, 1029, 1699, 965, 1133, 1951, 1621, 1069, 930, 1039, 1839]) # Return periods and probabilities return_periods = np.array([2, 3, 5, 10, 20, 50, 100, 200, 1000]) non_exceedance_probabilities = 1 - 1 / return_periods # Log-transform the data log_data = np.log(data) # L-moments method (using lmoments3) paramsmm = ld.pe3.lmom_fit(log_data) quantiles_lmom = np.exp(ld.pe3.ppf(non_exceedance_probabilities, paramsmm["skew"], paramsmm['loc'], paramsmm['scale'])) print("L-moments method results:") print(quantiles_lmom) # Maximum Likelihood Estimation method (using scipy.stats) params_mle = stats.pearson3.fit(log_data) quantiles_mle = np.exp(stats.pearson3.ppf(non_exceedance_probabilities, *params_mle)) print("\nMaximum Likelihood Estimation results:") print(quantiles_mle) # R results for comparison r_quantiles = np.array([1173.116, 1313.849, 1485.109, 1721.131, 1969.817, 2326.812, 2623.112, 2945.728, 3814.692]) print("\nR results:") print(r_quantiles) # Compare results print("\nDifference between Python L-moments and R results:") print(quantiles_lmom - r_quantiles)