> |
restart:
Digits:= 15:
St:= Statistics
:
|
> |
Zscore:= proc(V::Vector)
uses St=Statistics;
local mu:= St:-Mean(V), sigma:= St:-StandardDeviation(V);
mu, sigma, (V -~ mu)/sigma
end proc
:
|
> |
#Dependent variable:
Hvals:= <210.84, 218.19, 219.17, 222.27, 226.35, 228.45, 230.44, 234.65, 239.75>:
H__avg, H__sd, H__z:= Zscore(Hvals)
:
|
> |
#Independent variable:
Qvals:= <1800, 1750, 1730, 1700, 1650, 1635, 1600, 1570, 1500>:
Q__avg, Q__sd, Q__z:= Zscore(Qvals)
:
|
> |
#Proposed model, with (apparently) A known beforehand and Q__0, k__f, k__1, and Q__s to be
#determined.
model:=
subs(
Q= Q__N*Q__sigma + Q__mu,
H__N = (A*(1 - Q/Q__0) - k__f*Q^2 - k__1*(Q - Q__s)^2 - H__mu)/H__sigma
)
:
|
> |
pv:= St:-LinearFit(
[1,Q,Q^2], Q__z, H__z, Q, summarize,
output= [parametervector, residualsumofsquares]
);
|
Summary:
----------------
Model: .52025574e-1-1.0044218*Q-.58528770e-1*Q^2
----------------
Coefficients:
Estimate Std. Error t-value P(>|t|)
Parameter 1 0.0520 0.0414 1.2565 0.2556
Parameter 2 -1.0044 0.0319 -31.5038 0.0000
Parameter 3 -0.0585 0.0325 -1.8025 0.1215
----------------
R-squared: 0.9941, Adjusted R-squared: 0.9921
![[Vector(3, {(1) = 0.5202557353350474e-1, (2) = -1.0044218089795456, (3) = -0.5852877022519349e-1}), 0.47590683831025216e-1]](/view.aspx?sf=297321_Answer/a1cf2c739f5b7f214c2af08bb3e9bddf.gif)
Solve for the symbolic paarameters in the model by equating coefficients of the powers of Q. We have 3 coefficients and 4 parameters, so I'll just set the value of k__1 (arbitrary choice) to the MathCAD value, 1e-5.
> |
Known_params:=
[A, k__1, Q__mu, Q__sigma, H__mu, H__sigma]
=~ [350.156, 1e-5, Q__avg, Q__sd, H__avg, H__sd]
;
PV:= solve(
eval(
{seq(pv[1]=~ coeff~(expand(rhs(model)), Q__N, [0,1,2]))},
Known_params
)
);
|
![[A = 350.156, k__1 = 0.1e-4, Q__mu = HFloat(1659.4444444444443), Q__sigma = HFloat(95.01461875826149), H__mu = HFloat(225.5677777777778), H__sigma = HFloat(8.940258354457344)]](/view.aspx?sf=297321_Answer/bc7ce0f6ae987ac0d60b95e2a7e443c1.gif)

Compute the sum of the squared residuals. It should be the same as returned by the LinearFit command, allowing for some round off error.
> |
add(
eval(
eval(H__N - rhs(model), [PV[2][], Known_params[]])^2,
[H__N= H__z[k], Q__N= Q__z[k]]
),
k= 1..n
);
|

Compute the sum of the squared residuals using MathCAD's parameter values.
> |
add(
eval(
eval(
H__N - rhs(model),
[
([Q__0, k__f, Q__s]=~ [6962.32773, 2e-5, 2198.86103])[],
Known_params[]
]
)^2,
[H__N= H__z[k], Q__N= Q__z[k]]
),
k= 1..n
);
|


|