749701_zh-CN

liliconcom

长远考虑：如何使用 numpy.polyfit 重现 Vec Quadratic() 结果？

对于我的过失表示歉意，具有讽刺意味的是，这是我在 JSL 中取得个人进步的直接结果。

我所在组织中一位非常乐于助人的计算机科学家正在尝试将最初用 JSL 编写的内容翻译成 Python。具体来说，挑战在于如何重现 JSL 函数“ Vec Quadratic() ”的功能？另请参阅随附的示例表。

背景:

要翻译的 JSL 应用程序会为置信区间和预测区间创建公式列。然后，这些公式用于计算用户指定的某个x值（又称独立变量）的置信区间和预测区间值。

问题:

• numpy.polyfit （ https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html ）是正确的选择吗？
• 如果是，有谁知道如何使用 numpy.polyfit 来重现附表中 Lower & Upper 95% Mean & Indiv 列的结果，或者可以为我们指明正确的方向？

屏幕截图:

• 我们想要在Python中实现的结果（红框内的结果）
  

• 上面显示的列是如何生成的
  

自动化和脚本基础数据分析与建模回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

您是否向 Python 社区询问过此事。就我个人而言，努力帮助某人离开 JMP 对 JMP 社区来说似乎适得其反。

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

一方面确实如此。但另一方面，Python 已集成到 JMP18 中，这意味着一定有人开始接受这一现实。我希望社区中有人能够解释一下，因为我认为 JSL 用户更有可能知道“ Vec Quadratic() ”的 Python 等价形式，而不是相反。（当我问 Python 社区时，他们会告诉我什么？）

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

虽然 jmp.DataTable 对象将使一些用户更容易地导出 JMP 数据，但改进的 Python 集成侧重于使 JMP 用户受益。

• 它为刚接触 JMP 但熟悉 Python 的用户可以立即投入工作提供一个环境。
• Python 集成的设计重点是提高 JMP 作为工作流中心的能力。减少离开 JMP 完成工作的需要。希望/愿望是 JMP 环境变得更加方便、更具吸引力，工作将被带入 JMP 而不是从 JMP 导出。JMP/Python 集成不仅仅是 JMP 和 Python 的简单聚合。例如：
  • 可以使用外部 Python 包、parquet、hdf5 等打开 JMP 尚不支持的文件格式……
  • 允许 JMP 用户通过 pyopencl 和 pyCUDA 等软件包使用桌面 GPU。
  • 其他统计功能，如层次密度聚类 (HDBSCAN)
  • 使用 matplotlib、plotly 等软件包实现额外的绘图和图形功能...
  • 通过 Python 接口与其他软件系统集成。SAS、MatLab、R、...

JMP 统计人员和测试人员投入了大量精力来创建可靠、准确且经过充分测试的代码。这就是价值主张。我们专注于统计，因此您可以专注于您的应用程序，而不必重新发明统计“轮子”。

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

我并不怀疑 JMP 统计人员为编写可靠的代码付出了很多努力。但我怀疑的是，他们希望这个论坛不出现与 Python 相关的问题。

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

谢谢，太客气了！我们可能找到了解决方案： https://numpy.org/doc/stable/reference/generated/numpy.cov.html

编辑：这篇文章的标题具有误导性。numpy.polyfit 可能不适合模拟“ Vec Quadratic() ”。

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

我们欢迎在 JMP 博客中提出有关 Python 的问题。我最近回答了很多问题。包括示例代码，以及如何在 JMP 中充分利用 Python 的讨论。

@txnelson 的评论是向 Python 社区提问的，就像询问统计博客一样。不过，您不是询问如何复制 JSL 的“Vec Quadratic()”，而是询问如何绘制您想要的结果，如何在 Python 中进行您想要的统计。

至于回答这个具体的问题，我帮不了你。我不知道“Vec Quadratic()”的内部工作原理。但我可以告诉你脚本索引说了什么。Vec Quadratic(S,X) 的脚本索引定义表明它计算为 Vec Diag(X*S*X)。脚本索引中的 Vec Diag(X) 描述是将方阵的对角线元素作为向量返回。除非我遗漏了一些看似简单的矩阵数学的东西。

至于如何从在 JMP 中运行的 Python 程序中调用“Vec Quadriatic()”以便进一步处理结果……这是一种很有趣的问题。

您可以使用 numpy 数组表示 S 和 X，并且可以执行 Python Get(S)、Python Get(X) 将 numpy 数组值作为 JSL 矩阵放入 JSL 环境中。您可以调用 jmp.run_jsl( 'Vec Quadratic( S, X)' )。然后，从 JMP 中的 Python，您可以直接从 Python 中的数据表中获取列结果。请参阅 jmp.DataTable 和 jmp.DataTable.Column 对象的脚本索引。

源自脚本索引中的 JSL“Vec Quadratic()”示例。

import jmp
import numpy as np

s = np.array( [[1,3,5],[3,2,6], [5,6,1] ])
x = np.array( [[1,3,5],[2,4,6] ])
jmp.run_jsl('''
S = Python Get(s);
X = Python Get(x);
show(Vec Quadratic(S, X));
''')

给出结果

/*:

//:*/

S = Python Get(s);
X = Python Get(x);
show(Vec Quadratic(S, X));

/*:

Vec Quadratic(S, X) = [292, 528];

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

@Paul_Nelson ，谢谢。在发帖之前，我已经知道 JSL 函数“Vec Quadratic()”的作用了。抱歉，我在这篇文章中没有充分表达我的问题。这个论坛里并不是每个人都是英语母语人士。为了完整起见：

你的评论：

尽管不是询问如何复制 JSL 的“Vec Quadratic()”......

我原来的问题：

如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

---

“再现”一词的使用是否令人反感？

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

...你注意到它在矩阵操作方面存在一些问题:)

调整括号后，可得到：

import numpy as np

s = np.array( [[1.,3.,5.],[3.,2.,6.], [5.,6.,1.] ])
x = np.array( [[1.,3.,5.],[2.,4.,6.] ])

res= x @ (s @ x.T)
print(res)
print(res.diagonal())

... 适合

Names Default To Here( 1 );
exS = [1 3 5, 3 2 6, 5 6 1];
exX = [1 3 5, 2 4 6];

Print(r= exX * exS * exX`);

tmp=as list(r);
print(tmp[1][1],tmp[2][2])

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

@Ressel我认为您的请求中令人困惑的部分是 Vec Quadratic() 只是执行脚本索引中描述的矩阵数学运算。Vec Quadratic() 的结果只是矩阵运算的输出。

给出 95% 值的公式本身是由 Fit Y by X 平台生成的。'Fit y By X' 也是负责创建图表的。所以据我理解，这个问题更准确的说法是：如何在 Python 或 numpy 中对 y By x 进行双变量拟合，并设置置信度？

我认为Vec Quadratic() 的焦点是错误的，因为它没有创建公式或常量，它只是生成的公式中的数学运算。

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

@Paul_Nelson ，这是正确的，我后来才明白过来（并且厚颜无耻地利用这一点接受了我自己的答案作为解决方案）。关于这一点以及问题的确切表述，老实说，我当时真的不知道。学得越多，收获越多……非常感谢您的分享，祝您周末愉快。
（附言：作为一家成熟的统计软件公司，你几乎可以指望工程师没有计算才能成为商业模式的一部分。无论好坏。）

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

我赞赏 Jmp for Python 的开放性:)
正如@Paul_Nelson提到的 - 现在有了 Jmp18，如果 JMP 中缺少某个功能，每个用户都可以轻松地从 Python 中获取它 - 只需 2-3 行代码。非常有用！

另一方面，维基百科上关于矩阵操作、DOE、回归、神经网络等的一些知识不必担心被“窃取”。

凭借来自最先进、最新的 Python 软件包的“功能”，JMP 的“软技能”、无错误的出现、直观的界面和卓越的用户体验将成为其优势并使其在竞争中脱颖而出：

回复：愿望清单 - 新标签：用户体验

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

也许

NormalContour = sqrt(ChiSquare Quantile( 0.95, 2 ));

Ellipse的想法会给你指明正确的方向。我做的是图片，不是统计，所以我可能离题太远。但为了制作椭圆图，我必须回答一个类似的问题。

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

如果你问 Copilot：

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

此外，我认为通过 JMP 从事数据分析和自动化工作的用户不太可能熟练掌握 Python 或任何其他编程和脚本语言，这是很自然的。JMP 正在积极尝试招募这类用户，并使使用 JSL 编写脚本变得容易（我认为）。然而，一旦用户发展并学习，就可以引发对额外知识的渴望，至少在某些情况下是这样，在克服了理解一些基本计算概念的初始障碍之后。当然，这类用户的第一直觉是向社区寻求帮助，特别是在 JMP 18 中宣传 Python 集成之后。最后，专业数据科学家的任务是复制外行统计人员在 JMP 中开发的内容以用于实际商业环境，这一事实不言而喻。虽然，我可能全都错了。

回复：长远来看：如何使用 numpy.polyfit 来重现 Vec Quadratic() 结果？

如果有人感兴趣，下面的 Python 代码看起来将准确重现 JMP 在计算回归置信区间方面所做的工作。此代码块可以直接粘贴到 Jupyter 笔记本中。除了计算线性回归中每个 X 的置信区间的 Y 坐标外，它还会创建一个图，计算回归系数、截距、相关不确定性并制表结果。

不要相信这是我自己写的，但对于学习仍然有用。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# the data as array
x_coordinate = np.array([0, 1, 2, 3])
y_coordinate = np.array([0.0, 0.5, 1.2, 5.1])

# linear model
func = lambda x, a, b: a * x + b

# model fit
best_fit_ab, covar = curve_fit(func, x_coordinate, y_coordinate)
sigma_ab = np.sqrt(np.diagonal(covar))

# extract fitted parameters and uncertainties
a, b = best_fit_ab
a_err, b_err = sigma_ab
text_res = f"Best fit parameters:\na = {a:.10f} ± {a_err:.10f}\nb = {b:.10f} ± {b_err:.10f}"
print(text_res)

# calculate the standard error of the regression
residuals = y_coordinate - func(x_coordinate, a, b)
s_err = np.sqrt(np.sum(residuals**2) / (len(x_coordinate) - 2))

# define a function to calculate the confidence intervals
def predict_ci(x, func, popt, pcov, s_err, alpha=0.05): # adjust alpha here
    from scipy.stats import t
    df = len(x_coordinate) - 2
    t_val = t.ppf(1 - alpha / 2, df)  # t-value for the given confidence level and degrees of freedom
    y = func(x, *popt)
    ci = t_val * s_err * np.sqrt(1/len(x_coordinate) + (x - np.mean(x_coordinate))**2 / np.sum((x_coordinate - np.mean(x_coordinate))**2))
    return y - ci, y + ci

# generate points for high-resolution x-axis
hires_x = np.linspace(x_coordinate.min(), x_coordinate.max(), 100)

# calculate confidence intervals
ci_lower_orig, ci_upper_orig = predict_ci(x_coordinate, func, best_fit_ab, covar, s_err)

# tabulate the original data and confidence intervals
table_orig = np.column_stack((x_coordinate, y_coordinate, ci_lower_orig, ci_upper_orig))
print("\nTabulated values (original data):")
print("X Coordinate\tY Coordinate\t\tCI Lower\t\tCI Upper")
for row in table_orig:
    print(f"{row[0]:.2f}\t\t{row[1]:.10f}\t\t{row[2]:.10f}\t{row[3]:.10f}") # adjust decimals for tabulation

# plotting the data
plt.scatter(x_coordinate, y_coordinate, label='Data', color='blue')

# plotting the fitted model
plt.plot(hires_x, func(hires_x, *best_fit_ab), 'black', label='Fitted model')

# plotting the confidence intervals
ci_lower, ci_upper = predict_ci(hires_x, func, best_fit_ab, covar, s_err)
plt.fill_between(hires_x, ci_lower, ci_upper, color='black', alpha=0.15, label='Confidence interval')

# adding text with parameters
plt.text(x_coordinate.min() + 0.5, y_coordinate.max() - 0.5, text_res)

# customize plot
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('Y Coordinate vs. X Coordinate')
plt.legend()
plt.show()