如何计算scipy中曲线拟合的可能性?

2024-04-20

我有一个非线性模型拟合,如下所示:

深色实线是模型拟合,灰色部分是原始数据。

问题的简短版本:如何获得该模型拟合的可能性,以便我可以执行对数似然比测试?假设残差服从正态分布。

我对统计学比较陌生,我目前的想法是:

  1. 从曲线拟合中得到残差,并计算残差的方差;

  2. Use this equation Log-likelihood for normal distributions And plug in the variance of residual into sigma-squared, x_i as experiment and mu as model fit;

  3. 计算对数似然比。

谁能帮我解决这两个完整版问题?

  1. 我的方法正确吗? (我想是的,但如果能够确定的话那就太好了!)

  2. python/scipy/statsmodels 中是否有任何现成的函数可以为我执行此操作?


你的似然函数

它只是高斯分布的概率密度函数的对数之和。

是的可能性为您的残差拟合 mu 和 sigma,而不是可能性给定你的数据的模型。简而言之,你的方法是wrong.

正弦你正在做非线性最小二乘,按照@usethedeathstar已经提到的,你应该直接去F-test。 。考虑以下示例,修改自http://www.walkingrandomly.com/?p=5254 http://www.walkingrandomly.com/?p=5254,我们进行F-test using R。我们将讨论如何将其转化为python到底。

# construct the data vectors using c()
> xdata = c(-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9)
> ydata = c(0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001)
# some starting values
> p1 = 1
> p2 = 0.2
> p3 = 0.01

# do the fit
> fit1 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), start=list(p1=p1,p2=p2))
> fit2 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata)+p3*xdata, start=list(p1=p1,p2=p2,p3=p3))

# summarise
> summary(fit1)

Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
p1 1.881851   0.027430   68.61 2.27e-12 ***
p2 0.700230   0.009153   76.51 9.50e-13 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.08202 on 8 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 2.189e-06

> summary(fit2)

Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata) + p3 * xdata

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
p1  1.90108    0.03520  54.002 1.96e-10 ***
p2  0.70657    0.01167  60.528 8.82e-11 ***
p3  0.02029    0.02166   0.937     0.38    
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.08243 on 7 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: 2.476e-06

> anova(fit2, fit1)
Analysis of Variance Table

Model 1: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata) + p3 * xdata
Model 2: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1      7   0.047565                             
2      8   0.053813 -1 -0.0062473  0.9194 0.3696

这里我们有两个模型,fit1有2个参数,因此残差有8个自由度;fit2有一个附加参数,并且残差有 7 个自由度。模型 2 明显更好吗?否,F 值为 0.9194,上(1,7)自由度并不重要。

要获得方差分析表:残差 DF 很容易。残差平方和:0.08202*0.08202*8=0.05381 and 0.08243*0.08243*7=0.04756293(注意:“剩余标准误差:7 个自由度上的 0.08243”, ETC)。在python,你可以通过(y_observed-y_fitted)**2, since scipy.optimize.curve_fit()不返回残留物。

The F-ratio is 0.0062473/0.047565*7并获得 P 值:1-scipy.stats.f.cdf(0.9194, 1, 7).

把它们放在一起我们就有了python相等的:

In [1]:

import scipy.optimize as so
import scipy.stats as ss
xdata = np.array([-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9])
ydata = np.array([0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001])
def model0(x,p1,p2):
    return p1*np.cos(p2*x) + p2*np.sin(p1*x)
def model1(x,p1,p2,p3):
    return p1*np.cos(p2*x) + p2*np.sin(p1*x)+p3*x
p1, p2, p3 = 1, 0.2, 0.01
fit0=so.curve_fit(model0, xdata, ydata, p0=(p1,p2))[0]
fit1=so.curve_fit(model1, xdata, ydata, p0=(p1,p2,p3))[0]
yfit0=model0(xdata, fit0[0], fit0[1])
yfit1=model1(xdata, fit1[0], fit1[1], fit1[2])
ssq0=((yfit0-ydata)**2).sum()
ssq1=((yfit1-ydata)**2).sum()
df=len(xdata)-3
f_ratio=(ssq0-ssq1)/(ssq1/df)
p=1-ss.f.cdf(f_ratio, 1, df)
In [2]:

print f_ratio, p
0.919387419515 0.369574503394

正如@usethedeathstar 指出的:当你的残差呈正态分布时,非线性最小二乘IS最大可能性。因此F检验和似然比检验是等价的。因为,F-ratio 是似然比 λ 的单调变换.

或者以描述性的方式,请参阅:http://www.stata.com/support/faqs/statistics/chi-squared-and-f-distributions/ http://www.stata.com/support/faqs/statistics/chi-squared-and-f-distributions/

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

如何计算scipy中曲线拟合的可能性? 的相关文章

随机推荐

  • Nginx 站点已启用,站点可用:无法在 Ubuntu 12.04 中的配置文件之间创建软链接

    我正在尝试在 etc nginx 中的启用站点和可用站点目录中包含服务器块的配置文件之间创建软链接 我正在使用的命令是 sudo ln s sites available foo conf sites enabled 当我执行时 ls l
  • EXTJS 中选项卡面板的单击侦听器

    我在 extjs 中使用选项卡面板 我想在单击选项卡时显示警报 但我不知道如何 这就是我现在所做的 xtype tabpanel activeTab 0 region center items xtype panel title All i
  • 从 kubectl 输出显示失败的 pod

    我想写一个包装kubectl仅显示失败的 Pod 这意味着它应该只显示 Ready 列值不相同的项目 即0 1 0 2 1 2 2 3 etc kubectl get pods all namespaces NAMESPACE NAME R
  • 无法识别的配置节 system.serviceModel

    当我将网站发布到我的 Plesk 服务器时 出现以下错误 无法识别的配置节 system serviceModel 您的主机支持哪个版本的 NET 框架 The
  • 使用 Visual Studio 安装项目设置 InstallPath 注册表项

    我正在使用使用 Visual Studio 安装项目设计的 msi 安装程序来部署我的应用程序 如何将注册表项设置为应用程序的安装路径 实际上 当我在寻找同样的东西时 还提到了以下解决方案 在注册表项中使用 TARGETDIR
  • 如何从java类调用python脚本[重复]

    这个问题在这里已经有答案了 我有一个 java web 应用程序 我需要使用一个简单的网络爬虫从网页中读取 html 我在java中找不到任何简单的解决方案 但得到了一个非常简单的 python 脚本来解决我的问题 现在如何从我的 java
  • 在 Python 中使用 XLRD 迭代行和列

    我正在使用 python xlrd 模块来解析 Excel 文件 Excel 文件如下所示 Title A B C attribute 1 1 2 3 attribute 2 4 5 6 attribute 3 7 8 9 我想要以下格式的
  • 如何根据列值是否位于 Spark DataFrame 中的一组字符串中来过滤行

    是否有一种更优雅的方法根据字符串集中的值进行过滤 def myFilter actions Set String myDF DataFrame DataFrame val containsAction udf action String g
  • 如果当前行的宽度太窄,则将子级的溢出移至下一行

    EDIT 我正在构建一个简单的活动日历 使用 HTML CSS 目前正在处理多日活动 我是 HTML CSS 的初学者 有一个非常简单的问题 但我似乎找不到答案 如果没有 如何使子 div 溢出到下一行 div 等 屏幕上 或 div 行
  • 安全地运行 docker

    我知道 docker 守护进程需要以 root 身份运行 https docs docker com articles security 所以我被告知这可能会导致一些安全隐患 例如如果容器遭到破坏 攻击者可以更改主机的系统文件 发生攻击时
  • 我如何使其解密而不是加密?

    想知道如何从加密代码中获取此代码并使用相同的代码来创建解密 我知道这意味着我必须反转一些指令并重新排序 但我无法弄清楚哪些指令需要重新排序 哪些不需要 编辑 这是完整的函数 可以让事情变得更清晰一些 对堆栈溢出非常陌生 因此对于任何混淆表示
  • 同一台服务器的不同端口是否算跨域? (Ajax 方面)

    XMLHttpRequest 可以发送请求到http mydomain example 81 from http mydomain example 对于被视为具有相同来源的两个文档 协议 http https 域和端口 默认 80 或 xx
  • 以编程方式启动和停止 IIS Express

    我正在尝试用 C 构建一个小型应用程序 它应该启动 停止 IIS Express 工作进程 为此 我想使用 MSDN 上记录的官方 IIS Express API http msdn microsoft com en us library
  • React-native:如何在不单击react-native-maps中的标记的情况下显示工具提示

    我正在使用react native maps模块 我已经给出了纬度和经度值 并且我使用MapView Marker在单击标记时显示标记和工具提示 但是 现在我想在地图最初加载时显示工具提示 而无需单击标记 这是我的代码
  • 使用 eval 加载模块

    我在 Perl 和内置函数方面遇到了一些麻烦eval http perldoc perl org functions eval html 我浏览了网络 但找不到任何答案或示例代码 我想动态加载模块 在执行时间之前我不知道它们 module
  • 嵌套 RibbonApplicationMenuItem 时出错

    我想建立一个RibbonApplicationMenu 其中应嵌套一个RibbonApplicationMenuItem or RibbonApplicationSplitMenuItem 例如喜欢这个
  • Azure Web 角色如何在没有入口点的情况下运行?

    出于好奇 我打开了我的 Azure Web 角色项目 导航到包含以下内容的文件 RoleEntryPoint后代阶级和完全删除了该类定义 然后我打包了该角色并将其部署在 Azure 中 该角色启动时没有任何错误指示 这怎么可能行得通 除了
  • 如何找出 WPF 应用程序中的焦点在哪里?

    我的 WPF 应用程序中有一个搜索屏幕 该屏幕作为 TabControl 的 TabItem 中的 UserControl 实现 当用户切换到 搜索 选项卡时 我希望焦点进入一个特定字段 因此 我向 Xaml 中的 UserControl
  • AKSequencer 计数一或两个小节

    在当前序列开始播放之前需要播放 1 或 2 个小节进行计数 只需点击一下即可计入 能够做类似的事情会很酷 player sequencer setTime MusicTimeStamp 4 将时间设置为0 不起作用 使用 AKSequenc
  • 如何计算scipy中曲线拟合的可能性?

    我有一个非线性模型拟合 如下所示 深色实线是模型拟合 灰色部分是原始数据 问题的简短版本 如何获得该模型拟合的可能性 以便我可以执行对数似然比测试 假设残差服从正态分布 我对统计学比较陌生 我目前的想法是 从曲线拟合中得到残差 并计算残差的