拟合分布、拟合优度、p 值。是不是可以使用 Scipy (Python) 做到这一点？

Posted 2023-03-12

【中文标题】拟合分布、拟合优度、p 值。是不是可以使用 Scipy (Python) 做到这一点？

【英文标题】：Fitting distributions, goodness of fit, p-value. Is it possible to do this with Scipy (Python)?拟合分布、拟合优度、p 值。是否可以使用 Scipy (Python) 做到这一点？

【发布时间】：2011-09-30 17:23:22

【问题描述】：

简介：我是一名生物信息学家。在我对所有人类基因（大约 20 000 个）进行的分析中，我搜索了一个特定的短序列基序，以检查该基序在每个基因中出现的次数。

基因以四个字母（A、T、G、C）的线性序列“写入”。例如：CGTAGGGGGTTTAC...这是遗传密码的四个字母表，就像每个细胞的秘密语言，是DNA实际存储信息的方式。

我怀疑某些基因中特定短基序序列 (AGTGGAC) 的频繁重复对于细胞中特定的生化过程至关重要。由于基序本身非常短，因此使用计算工具很难区分基因中真正的功能示例和那些偶然看起来相似的示例。为了避免这个问题，我得到所有基因的序列并连接成一个字符串并打乱。存储了每个原始基因的长度。然后对于每个原始序列长度，通过从级联序列中重复随机选择A或T或G或C并将其转移到随机序列来构建随机序列。这样，生成的随机序列集具有相同的长度分布，以及整体 A、T、G、C 组成。然后我在这些随机序列中搜索主题。我将这个程序执行了 1000 次，然后取平均值。

15000 个不包含给定基序的基因 5000 个基因包含 1 个基序 3000 个基因，包含 2 个基序 1000 个基因，包含 3 个基序 ... 1个基因包含6个基序

因此，即使对真正的遗传密码进行 1000 次随机化，也没有任何基因具有超过 6 个基序。但在真正的遗传密码中，有少数基因包含超过 20 次出现的基序，这表明这些重复可能是功能性的，不可能完全偶然地找到如此丰富的它们。

问题：

我想知道在我的分布中找到一个基序出现 20 次的基因的概率。所以我想知道偶然发现这样一个基因的概率。我想在 Python 中实现这个，但我不知道如何。

我可以用 Python 做这样的分析吗？

任何帮助将不胜感激。

【问题讨论】：

请注意，您已大幅修改您的问题。是否可以将此问题恢复为您的原始问题以及所有新细节的明确“更新”部分？或者也许只是一个新问题？谢谢

你可以考虑在Biostar上提问

我问一个新问题：***.com/questions/6620471/…

【参考方案1】：

In SciPy documentation 您将找到所有已实现的连续分布函数的列表。每一个都有a fit() method，返回对应的shape参数。

即使您不知道要使用哪种分布，您也可以同时尝试多种分布并选择更适合您的数据的分布，如下面的代码所示。请注意，如果您不了解分布，则可能难以拟合样本。

import matplotlib.pyplot as plt
import scipy
import scipy.stats
size = 20000
x = scipy.arange(size)
# creating the dummy sample (using beta distribution)
y = scipy.int_(scipy.round_(scipy.stats.beta.rvs(6,2,size=size)*47))
# creating the histogram
h = plt.hist(y, bins=range(48))

dist_names = ['alpha', 'beta', 'arcsine',
              'weibull_min', 'weibull_max', 'rayleigh']

for dist_name in dist_names:
    dist = getattr(scipy.stats, dist_name)
    params = dist.fit(y)
    arg = params[:-2]
    loc = params[-2]
    scale = params[-1]
    if arg:
        pdf_fitted = dist.pdf(x, *arg, loc=loc, scale=scale) * size
    else:
        pdf_fitted = dist.pdf(x, loc=loc, scale=loc) * size
    plt.plot(pdf_fitted, label=dist_name)
    plt.xlim(0,47)
plt.legend(loc='upper left')
plt.show()

参考资料：

- Distribution fitting with Scipy

- Fitting empirical distribution to theoretical ones with Scipy (Python)?

【讨论】：

上面的代码向我抛出以下消息错误：“AttributeError: 'module' object has no attribute 'arcsineweibull_min'” on the "dist = getattr(scipy.stats, dist_name)" 语句。我的版本是：scipy 是 0.13.3，numpy 是 1.8.0，matplotlib 是 1.3.1。

@srodriguex 谢谢！有一个小错字，我刚刚修复了它

@SaulloCastro 我如何将这个fit() 函数应用于3D 表面拟合，我刚刚使用scipy.linalg.lstsq 实现了？我如何确认我对数据的拟合度很好。谢谢。

这个例子不能正常工作了，你能更新一下吗？

@Khris 我会试着找点时间来处理它。如果您已经知道与新 SciPy 包不兼容的内容，请随意编辑...

【参考方案2】：

import numpy as np
import pandas as pd
import scipy.stats as st
import statsmodels.api as sm
import matplotlib as mpl
import matplotlib.pyplot as plt
import math
import random

mpl.style.use("ggplot")

def danoes_formula(data):
    """
    DANOE'S FORMULA
    https://en.wikipedia.org/wiki/Histogram#Doane's_formula
    """
    N = len(data)
    skewness = st.skew(data)
    sigma_g1 = math.sqrt((6*(N-2))/((N+1)*(N+3)))
    num_bins = 1 + math.log(N,2) + math.log(1+abs(skewness)/sigma_g1,2)
    num_bins = round(num_bins)
    return num_bins

def plot_histogram(data, results, n):
    ## n first distribution of the ranking
    N_DISTRIBUTIONS = k: results[k] for k in list(results)[:n]

    ## Histogram of data
    plt.figure(figsize=(10, 5))
    plt.hist(data, density=True, ec='white', color=(63/235, 149/235, 170/235))
    plt.title('HISTOGRAM')
    plt.xlabel('Values')
    plt.ylabel('Frequencies')

    ## Plot n distributions
    for distribution, result in N_DISTRIBUTIONS.items():
        # print(i, distribution)
        sse = result[0]
        arg = result[1]
        loc = result[2]
        scale = result[3]
        x_plot = np.linspace(min(data), max(data), 1000)
        y_plot = distribution.pdf(x_plot, loc=loc, scale=scale, *arg)
        plt.plot(x_plot, y_plot, label=str(distribution)[32:-34] + ": " + str(sse)[0:6], color=(random.uniform(0, 1), random.uniform(0, 1), random.uniform(0, 1)))
    
    plt.legend(title='DISTRIBUTIONS', bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.show()

def fit_data(data):
    ## st.frechet_r,st.frechet_l: are disbled in current SciPy version
    ## st.levy_stable: a lot of time of estimation parameters
    ALL_DISTRIBUTIONS = [        
        st.alpha,st.anglit,st.arcsine,st.beta,st.betaprime,st.bradford,st.burr,st.cauchy,st.chi,st.chi2,st.cosine,
        st.dgamma,st.dweibull,st.erlang,st.expon,st.exponnorm,st.exponweib,st.exponpow,st.f,st.fatiguelife,st.fisk,
        st.foldcauchy,st.foldnorm, st.genlogistic,st.genpareto,st.gennorm,st.genexpon,
        st.genextreme,st.gausshyper,st.gamma,st.gengamma,st.genhalflogistic,st.gilbrat,st.gompertz,st.gumbel_r,
        st.gumbel_l,st.halfcauchy,st.halflogistic,st.halfnorm,st.halfgennorm,st.hypsecant,st.invgamma,st.invgauss,
        st.invweibull,st.johnsonsb,st.johnsonsu,st.ksone,st.kstwobign,st.laplace,st.levy,st.levy_l,
        st.logistic,st.loggamma,st.loglaplace,st.lognorm,st.lomax,st.maxwell,st.mielke,st.nakagami,st.ncx2,st.ncf,
        st.nct,st.norm,st.pareto,st.pearson3,st.powerlaw,st.powerlognorm,st.powernorm,st.rdist,st.reciprocal,
        st.rayleigh,st.rice,st.recipinvgauss,st.semicircular,st.t,st.triang,st.truncexpon,st.truncnorm,st.tukeylambda,
        st.uniform,st.vonmises,st.vonmises_line,st.wald,st.weibull_min,st.weibull_max,st.wrapcauchy
    ]
    
    MY_DISTRIBUTIONS = [st.beta, st.expon, st.norm, st.uniform, st.johnsonsb, st.gennorm, st.gausshyper]

    ## Calculae Histogram
    num_bins = danoes_formula(data)
    frequencies, bin_edges = np.histogram(data, num_bins, density=True)
    central_values = [(bin_edges[i] + bin_edges[i+1])/2 for i in range(len(bin_edges)-1)]

    results = 
    for distribution in MY_DISTRIBUTIONS:
        ## Get parameters of distribution
        params = distribution.fit(data)
        
        ## Separate parts of parameters
        arg = params[:-2]
        loc = params[-2]
        scale = params[-1]
    
        ## Calculate fitted PDF and error with fit in distribution
        pdf_values = [distribution.pdf(c, loc=loc, scale=scale, *arg) for c in central_values]
        
        ## Calculate SSE (sum of squared estimate of errors)
        sse = np.sum(np.power(frequencies - pdf_values, 2.0))
        
        ## Build results and sort by sse
        results[distribution] = [sse, arg, loc, scale]
        
    results = k: results[k] for k in sorted(results, key=results.get)
    return results
        
def main():
    ## Import data
    data = pd.Series(sm.datasets.elnino.load_pandas().data.set_index('YEAR').values.ravel())
    results = fit_data(data)
    plot_histogram(data, results, 5)

if __name__ == "__main__":
    main()