topic Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results? in Discussions

Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Ressel — Thu, 25 Apr 2024 13:41:15 GMT

Apologies for the transgression, which is; ironically; a direct result of my personal progress in JSL.

A very helpful computer scientist in my organization is trying to translate something initially written in JSL to Python. Specifically, the challenge is how to reproduce what the JSL function 'Vec Quadratic()' is doing? See also the attached example table .

Background:

The JSL application to be translated creates formula columns for the confidence and prediction intervals. The formulas are then used to calculate the values of the confidence and prediction intervals for some user specified value of x, aka independent variable.

Questions:

Is numpy.polyfit (https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html) the right choice for this?
If it is, anyone who knows how to use numpy.polyfit to reproduce the results for Lower & Upper 95% Mean & Indiv columns in the attached table or can point us into the right direction?

Screenshots for orientation:

What we want to achieve in Python (results inside the red frame)

How the columns shown above were generated

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

txnelson — Thu, 25 Apr 2024 14:00:49 GMT

Have you ask the Python community about this. Personally, it seems counter productive to the JMP Community to put effort into assisting in moving someone away from JMP.

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Ressel — Thu, 25 Apr 2024 19:10:39 GMT

True on the one hand. On the other hand, though, Python is integrated into JMP18, which suggests there must be people getting to terms with this reality. I was hoping for someone in the community being able to explain, since I consider it infinitely more likely that a JSL user would know the Python equivalent of 'Vec Quadratic()', rather than vice versa. (And what will the Python community tell me when I ask them?)

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Paul_Nelson — Fri, 26 Apr 2024 00:11:10 GMT

While the jmp.DataTable object will enable some user's an easier path to export JMP data, the improved Python integration is focused on benefiting the JMP user.

It provides an environment that users new to JMP, but familiar with Python can be productive immediately.
The design of the Python integration focuses on improving JMP's capability at being the hub of the workflow. At reducing the need to leave JMP to get the job done. The hope / desire is that the JMP environment becomes more convenient, more enticing, that work will be brought into JMP rather than exported away from JMP. That JMP / Python integration becomes greater than simply the aggregation of JMP and Python. For example:
- File formats not yet supported by JMP can be opened using external Python packages, parquet, hdf5, ...
- Allow JMP users to utilize desktop GPUs through packages such as pyopencl, and pyCUDA.
- Additional statistical capabilities such as hierarchical density clustering (HDBSCAN)
- Additional plotting and graphing capabilities with packages such as matplotlib, plotly, ...
- Integration with other software systems through the Python interface. SAS, MatLab, R, ...

The JMP statisticians and testers put a lot of effort into creating code that is reliable, accurate and well tested. That is the value proposition. We focus on the statistics so you can focus on your application, rather than having to reinvent the statistical 'wheel'.

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Ressel — Fri, 26 Apr 2024 08:44:46 GMT

I don't doubt that JMP statisticians put a lot of effort into creating reliable code. What I doubt is the expectation to keep this forum free from Python related questions.

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Craige_Hales — Fri, 26 Apr 2024 11:08:40 GMT

Maybe the

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

idea from will point you in the right direction. I do pictures, not statistics, so I might be way off. But to make the ellipse picture I had to answer a similar question.

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Ressel — Fri, 26 Apr 2024 13:56:37 GMT

Thanks, that's very kind! It is possible we found a solution: https://numpy.org/doc/stable/reference/generated/numpy.cov.html

Edit: The title of this posting is misleading. numpy.polyfit is probably no suitable for emulating 'Vec Quadratic()'.

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Paul_Nelson — Fri, 26 Apr 2024 14:17:51 GMT

We welcome Python questions in the JMP blogs. I have been answering a lot of them lately. Including sample code, and discussions on how best to use Python in JMP to its best advantage.

@txnelson comment on asking the Python community is a valid one, as would asking statistical blogs. Though instead of asking how to replicate JSL's 'Vec Quadratic()' you phrase it in terms of asking how to graph the results you want, how to do the statistics you want in Python.

As for answering this specific ask, I can't help you. I have no knowledge of how 'Vec Quadratic()' works internally. But I can tell you what the Scripting Index says. The scripting index definition of Vec Quadratic( S, X ) states it evaluates as Vec Diag ( X * S * X ). Vec Diag( X ) description in the Scripting index is Returns the diagonal elements of the square matrix as a vector. Unless I'm missing something that appears to be straightforward matrix math.

As for how do I call 'Vec Quadriatic()' from my Python program running in JMP so I can further process the results... That is the kind of question that would be interesting.

You can use a numpy array for S and X and you can do the Python Get(S), Python Get(X) to get the numpy array values into JSL environment as a JSL Matrix. And you can call jmp.run_jsl( 'Vec Quadratic( S, X)' ). Then from Python in JMP you can get the column results directly from the data table in Python. See the Scripting Index on jmp.DataTable and jmp.DataTable.Column objects.

Derived from the JSL "Vec Quadratic()' example in the Scripting Index.

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)); ''')

giving results

/*:

//:*/

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

/*:

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

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Ressel — Fri, 26 Apr 2024 17:32:59 GMT

@Paul_Nelson, thank you. I already knew what the JSL function 'Vec Quadratic()' did before posting. Apologies for inadequately phrasing my question in this post. Not everyone is this forum is an English native speaker. And just for completeness:

Your comment:

Though instead of asking how to replicate JSL's 'Vec Quadratic()' ...

My original question:

How to use numpy.polyfit to reproduce Vec Quadratic() results?

Is it the use of "reproduce" that is offensive?

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

hogi — Fri, 26 Apr 2024 17:44:03 GMT

if you ask Copilot:

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

hogi — Fri, 26 Apr 2024 17:48:21 GMT

... you notice that it has some issues with matrix manipulations :)

after adjusting the brackets, you get:

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())

... which fits to

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])

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Paul_Nelson — Fri, 26 Apr 2024 17:51:27 GMT

@Ressel I think the confusing part of your request is that Vec Quadratic() is simply performing matrix math as described by the Scripting Index. The results of Vec Quadratic() is simply the output of matrix operations.

The formulas themselves, which give the 95% values, are generated by the Fit Y by X platform. 'Fit y By X' is also what is responsible for creating the plot. So the question as I understand it, is more accurately; How do I do a Bivariate Fit of y By x with confidence levels in Python or numpy?

I think the Vec Quadratic() is the wrong focus since it didn't create the formula or the constants, it's just a mathematical operation in the generated formula.

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Ressel — Fri, 26 Apr 2024 19:09:35 GMT

@Paul_Nelson, this is correct, which I only figured out later (and shamelessly exploited to accept my own answer as solution). With regards to this as well as the exact phrasing of the question, I honestly didn't have it in me. The more one learns ... Thank you very much for sharing & have a good weekend.
(P.s.: As a sophisticated, statistical software company you'd almost count on engineers without computational talent being part of the business model. For better or worse.)

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

hogi — Fri, 26 Apr 2024 18:49:17 GMT

I applaud the openness of Jmp for Python :)
As @Paul_Nelson mentions - now with Jmp18, if a functionality is missing in JMP, every user can easily grab it from Python - with 2-3 lines of code. amazingly useful!

On the other hand, some Wikipedia knowledge about matrix manipulation, DOEs, regression, neural networks, etc. is nothing to worry about being "stolen".

With the "functionality" coming from state-of-the-art, up-to-date Python packages, it will be the "soft skills", the absence of bugs, the intuitive interface, the superior user experience that will speak for JMP and separate it from the competition:

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Ressel — Fri, 26 Apr 2024 19:05:55 GMT

Moreover, I think it's only natural that users getting into data analytical & automation work via JMP are unlikely to have a firm grasp of Python, or any other programing & scripting language for that matter. JMP is proactively trying to recruit this type of user and makes it (I think) easy to step into scripting in JSL. Once the user develops and learns, though, this can be expected to trigger the wish for additional knowledge, at least in some cases and after the initial barrier of understanding a few basic, computational concepts is overcome. Of course the first intuition for this kind of user will be to ask the community for help, especially after the advertised Python integration in JMP 18. Finally, the fact that a professional data scientist is tasked with copying what a lay statistician has developed in JMP for use in an actual business setting speaks for itself. Although, I might be getting it all wrong.

Re: Long shot: How to use numpy.polyfit to reproduce Vec Quadratic() results?

Ressel — Thu, 23 May 2024 18:51:22 GMT

If anyone's interested, the below Python code looks like it will exactly reproduce what JMP is doing in terms of calculating the confidence interval around a regression. This code block can be pasted into a Jupyter notebook directly. Besides calculating the Y- coordinates of the confidence interval for every X in the linear regression it does also create a plot, calculates the regression coefficient, intercept, associated uncertainties and tabulates the results.

Don't believe I wrote this myself, but still useful for learning.

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()