## fitting – Finding astroids – Mathematica Stack Exchange

An astroid is a particular mathematical curve: a hypocycloid with four cusps.

Given a lot of data points I’m trying to fit such a curve to the data.
As an astroid is also an algebraic curve of grade 6 I used LinearModelFit:

points=Import("https://pastebin.com/raw/Bi3DPvNj","Data");
lin ={ #1 + #2, #1^2 + #2^2, #1^3 + #2^3, #1^4 + #2^4, #1^5 + #2^5, #1^6 + #2^6, #1 #2, #1^2 #2 + #1 #2^2, #1^3 #2 + #1 #2^3, #1^4 #2 + #1 #2^4, #1^5 #2 + #1 #2^5, #1^2 #2^2, #1^3 #2^2 + #1^2 #2^3, #1^4 #2^2+ #1^2 #2^4, #1^3 #2^3} & @@@ points;
lm = LinearModelFit(lin,Prepend(Table(Subscript(a,i),{i,1,14}),1),Table(Subscript(a,i),{i,1,14}));
Quiet@lm("ANOVATable")


Which results in

w(x_, y_) :=lm("BestFitParameters").{1, x+y,  x^2+y^2, x^3+ y^3, x^4+ y^4, x^5+ y^5, x^6+y^6, x y, x y^2+ x^2 y, x y^3+ x^3 y, x y^4+x^4 y, x y^5+x^5 y, x^2 y^2, x^2 y^3+x^3 y^2, x^2 y^4+ x^4 y^2}-x^3 y^3;
ContourPlot(w(x, y) == 0, {x, 0.8, 1.0}, {y, 0.8, 1.0}, Prolog -> Point(points),ContourStyle->{Red},PlotPoints->109,ImageSize->Large)


Actually, the fit is okay, except near 2 of the cusps (green). I’m pretty sure the fit could be improved, but how?

## fitting – Using LinearModelFit or Predict to attempt a multiple linear regression model in Mathematica

I’m using Mathematica for a project, in which I would like to fit a multiple linear regression model to my data (which is made up of both numerical and categorical variables).

I have been trying to use the function LinearModelFit to do this, but I have had no luck. Based on the MMA documentation available, the examples they’ve given only use simple numerical data such as {{0, 1}, {1, 0}, {3, 2}, {5, 4}}. In my case, I have 3 numerical variables, 4 categorical variables and my response variable, so I haven’t been able to achieve this with LinearModelFit.

I then tried to use Predict(list1 -> list 2, Method->”LinearRegression”), with all my variable inputs in a list which I inserted in place of list1, and then my list of response variable values in a list which I inserted in place of list 2. This resulted in an error “Incompatible variable type (!(“Numerical”)) and variable value”.

I’m wondering does anyone know if there is a different function I should be using or if I’m perhaps not using the 2 functions mentioned above correctly?

Thanks.

## fitting – How to perform a LinearModelFit on a data-set with units of measurement?

Thanks for contributing an answer to Mathematica Stack Exchange!

But avoid

• Making statements based on opinion; back them up with references or personal experience.

Use MathJax to format equations. MathJax reference.

## fitting – non linear Fit with a parametric value

i have 4 sets of data of JCS as a function of T at 4 different B values
the equation is
JCS = JCST*Exp(-(B/BSC));

JCST = JCS0*(1 – T/78.4)^alpha2;

BSC = BSC0*Exp(-(T/TSC));

where JCS0 BSC and TCS have to be the same for all the dataset while B changes for each data but they are known.
i use to do a fit whit a single data set setting B every fit, but I get different best values for JCS0 BSC and TCS.
how can I perform a single Fit of JCS vs T with all 4 sets of data including the fact that B is not a fitting parameter but is a know value?

thanx

## fitting – How to test if NonlinearModelFit parameters are different from each other

Thanks for contributing an answer to Mathematica Stack Exchange!

But avoid

• Making statements based on opinion; back them up with references or personal experience.

Use MathJax to format equations. MathJax reference.

## plotting – Fitting two convoluted peaks

If I have the following data:

https://pastebin.com/ya03z9bP

which plotted as:

ListLinePlot(datawithnoliquidline,
PlotStyle -> Directive(Thick, Black),
PlotRange -> {{40, 110}, {-0.02, All}}, Frame -> True,
FrameStyle -> 14, Axes -> False, GridLines -> Automatic,
GridLinesStyle -> Lighter(Gray, .8),
FrameTicks -> {Automatic, Automatic},
FrameLabel -> (Style(#, 20, Bold) & /@ {"T ((Degree)C)",
Row({"!(*SubscriptBox((C), (P)))", " (", " J/gK)"})}),
LabelStyle -> {Black, Bold, 14})


gives two peaks as in the picture:

Questions:
1)How can I correctly fit these two peaks?
2) How can I calculate the areas of the two fitted peaks?

My approach, using gaussian fit (which doesn’t seem correct) is as follows:

ff2(x_, areaa1_, areaa3_, siga1_, siga3_, meda1_, meda3_) :=
areaa1 PDF(NormalDistribution(meda1, siga1), x) +
areaa3 PDF(NormalDistribution(meda3, siga3), x) ;

ma1guess = 6;
ma3guess = 1.3;
siga1guess = 4;
siga3guess = 3;
meda1guess = 82;
meda3guess = 97;

averagenematicarea = 1.2;(*Average nematic area from all three DSC
runs*)
STDnematicarea = 0.2; (*Standard deviation of nematic
area*)(*Aceptable shift above and below the nematic area where the
fits will be constraint (e.g if shift is 0.2 and area=1.2, then the
fits would be constraint between 1 J/g and 1.4 J/g)*)

averagesmecticarea = 1;(*Average smectic area from all three DSC runs*)

STDsmecticarea = 0.2; (*Standard deviation of smectic
area*)(*Aceptable shift above and below the nematic area where the
fits will be constraint (e.g if shift is 0.2 and area=1.2, then the
fits would be constraint between 1 J/g and 1.4 J/g)*)

averagenematiconset = 88.2;(*Average nematic onset from all three DSC
runs*)
STDnematiconset = 0.2; (*Standard deviation of nematic
onset*)(*Aceptable shift above and below the nematic onset where the
fits will be constraint (e.g if shift is 0.5 and onset=88, then the
fits would be constraint between 87.5C and 88.5C)*)

nlm3 = NonlinearModelFit(
datawithnoliquidline, {ff2(x, areaa1, areaa3, siga1, siga3, meda1,
meda3), areaa3 >= 0, areaa1 >= 0, 0 <= siga1 <= 20,
0 <= siga3 <= 20, 60 < meda1 < 85,
meda3 - 2*siga3 > 88.6}, {{areaa1, ma1guess}, {areaa3,
ma3guess}, {siga1, siga1guess}, {siga3, siga3guess}, {meda1,
meda1guess}, {meda3, meda3guess}}, x);

fp = nlm3("BestFitParameters");

p1 =(*Original data*)
ListLinePlot(datawithnoliquidline,
PlotStyle -> Directive(Thick, Black),
PlotRange -> {{40, 110}, {-0.02, All}}, Frame -> True,
FrameStyle -> 14, Axes -> False, GridLines -> Automatic,
GridLinesStyle -> Lighter(Gray, .8),
FrameTicks -> {Automatic, Automatic},
FrameLabel -> (Style(#, 20, Bold) & /@ {"T ((Degree)C)",
Row({"!(*SubscriptBox((C), (P)))", " (", " J/gK)"})}),
LabelStyle -> {Black, Bold, 14});
p2 = Plot({nlm3(x),
areaa3 PDF(NormalDistribution(meda3, siga3), x) /.
fp, +areaa1 PDF(NormalDistribution(meda1, siga1), x) /. fp}, {x,
40, 110},
PlotStyle -> {Directive(Red, Dashing({0.02, 0.04}),
AbsoluteThickness(5)), Directive(Green, AbsoluteThickness(2)),
Directive(Orange, AbsoluteThickness(2)),
Directive(Blue, AbsoluteThickness(2)), Directive(Pink, Dashed),
Directive(Cyan, Dashed)}, PlotRange -> All);

Show(p1, p2)


which gives the following:

So, my problem is really fitting correctly the second peak, which I am not sure how to do.

## Fitting data to an expression

I am trying to fit the following data to an expression of the form shown in the code below but it doesn’t seem to fit that well and I don’t know why it isn’t as this kind of data is routinely fit to such expressions. I would greatly appreciate any help or pointer on what I am doing wrong. Thanks.

data = {{2., 15858.5}, {2.1, 10855.9}, {2.2, 7346.87}, {2.3, 4908.57}, {2.4, 3229.68}, {2.5, 2084.4}, {2.6, 1310.79}, {2.7, 793.984}, {2.8, 453.257}, {2.9, 232.31}, {3., 92.1445}, {3.1, 5.8997}, {3.2, -44.8126}, {3.3, -72.4833}, {3.4, -85.5145}, {3.5, -89.4885}, {3.6, -88.0608}, {3.7, -83.587}, {3.8, -77.5604}, {3.9, -70.9073}, {4., -64.1871}, {4.1, -57.7239}, {4.2, -51.6913}, {4.3, -46.1686}, {4.4, -41.1779}, {4.5, -36.7066}, {4.6, -32.7236}, {4.7, -29.1886}, {4.8, -26.059}, {4.9, -23.2921}, {5., -20.8469}, {5.1, -18.6868}, {5.2, -16.7775}, {5.3, -15.0885}, {5.4, -13.5933}, {5.5, -12.2677}, {5.6, -11.0904}, {5.7, -10.0433}, {5.8, -9.11012}, {5.9, -8.27735}, {6., -7.53272}, {6.1, -6.86538}, {6.2, -6.26625}, {6.3, -5.72797}, {6.4, -5.24286}, {6.5, -4.80566}, {6.6, -4.41044}, {6.7, -4.05333}, {6.8, -3.72946}, {6.9, -3.43566}, {7., -3.1688}, {7.1, -2.92572}, {7.2, -2.70467}, {7.3, -2.50286}, {7.4, -2.31854}, {7.5, -2.14996}, {7.6, -1.99571}, {7.7, -1.85441}, {7.8, -1.72465}, {7.9, -1.60573}, {8., -1.49591}, {8.1, -1.39518}, {8.2, -1.30215}, {8.3, -1.21645}, {8.4, -1.13741}, {8.5, -1.06396}, {8.6, -0.996458}, {8.7, -0.933501}, {8.8, -0.875442}, {8.9, -0.821579}, {9., -0.771564}, {9.1, -0.724697}, {9.2, -0.681676}, {9.3, -0.641105}, {9.4, -0.603681}, {9.5, -0.568705}, {9.6, -0.535828}, {9.7, -0.505399}, {9.8, -0.476719}, {9.9, -0.450137}, {10., -0.425305}, {10.1, -0.401871}, {10.2, -0.380186}, {10.3, -0.35955}, {10.4, -0.340314}, {10.5, -0.322126}, {10.6, -0.305338}, {10.7, -0.289249}, {10.8, -0.274209}, {10.9, -0.260219}, {11., -0.246928}, {11.1, -0.234337}, {11.2, -0.222445}, {11.3, -0.211603}, {11.4, -0.20111}, {11.5, -0.190967}, {11.6, -0.181874}, {11.7, -0.17278}, {11.8, -0.164386}, {11.9, -0.156691}, {12., -0.149346}, {12.1, -0.142001}, {12.2, -0.135356}, {12.3, -0.12906}, {12.4, -0.123114}, {12.5, -0.117518}, {12.6, -0.111922}, {12.7, -0.106676}, {12.8, -0.102129}, {12.9, -0.0972324}, {13., -0.0930354}, {13.1, -0.0888383}, {13.2, -0.0846412}, {13.3, -0.0811436}, {13.4, -0.0772963}, {13.5, -0.0741485}, {13.6, -0.0706509}, {13.7, -0.0675031}, {13.8, -0.064705}, {13.9, -0.061907}, {14., -0.0591089}, {14.1, -0.0566606}, {14.2, -0.0542123}, {14.3, -0.051764}, {14.4, -0.0496655}, {14.5, -0.047567}, {14.6, -0.0454684}, {14.7, -0.0433699}, {14.8, -0.0416211}, {14.9, -0.0398723}, {15., -0.0381235}, {15.1, -0.0363747}, {15.2, -0.0349757}, {15.3, -0.0335767}, {15.4, -0.0321776}, {15.5, -0.0307786}, {15.6, -0.0293796}, {15.7, -0.0279806}, {15.8, -0.0269313}, {15.9, -0.025882}, {16., -0.024483}, {16.1, -0.0234337}, {16.2, -0.0223844}, {16.3, -0.0216849}, {16.4, -0.0206357}, {16.5, -0.0195864}, {16.6, -0.0188869}, {16.7, -0.0181874}, {16.8, -0.0171381}, {16.9, -0.0164386}, {17., -0.0157391}, {17.1, -0.0150396}, {17.2, -0.01434}, {17.3, -0.0136405}, {17.4, -0.012941}, {17.5, -0.0125913}, {17.6, -0.0118917}, {17.7, -0.0111922}, {17.8, -0.0108425}, {17.9, -0.010143}, {18., -0.0097932}, {18.1, -0.00909368}, {18.2, -0.00874393}, {18.3, -0.00839417}, {18.4, -0.00769465}, {18.5, -0.0073449}, {18.6, -0.00699514}, {18.7, -0.00664538}, {18.8, -0.00629563}, {18.9, -0.00594587}, {19., -0.00559611}, {19.1, -0.00524636}, {19.2, -0.0048966}, {19.3, -0.00454684}, {19.4, -0.00419708}, {19.5, -0.00384733}, {19.6, -0.00349757}, {19.7, -0.00349757}, {19.8, -0.00314781}, {19.9, -0.00279806}, {20., -0.0024483}, {20.1, -0.0024483}, {20.2, -0.00209854}, {20.3, -0.00174879}, {20.4, -0.00174879}, {20.5, -0.00139903}, {20.6, -0.00139903}, {20.7, -0.00104927}, {20.8, -0.000699514}, {20.9, -0.000699514}, {21., -0.000349757}, {21.1, -0.000349757}, {21.2, 0.}, {21.3, 0.}, {21.4, 0.000349757}, {21.5, 0.000349757}, {21.6, 0.000349757}, {21.7, 0.000699514}, {21.8, 0.000699514}, {21.9, 0.00104927}, {22., 0.00104927}, {22.1, 0.00104927}, {22.2, 0.00139903}, {22.3, 0.00139903}, {22.4, 0.00174879}, {22.5, 0.00174879}, {22.6, 0.00174879}, {22.7, 0.00174879}, {22.8, 0.00209854}, {22.9, 0.00209854}, {23., 0.00209854}, {23.1, 0.0024483}, {23.2, 0.0024483}, {23.3, 0.0024483}, {23.4, 0.0024483}, {23.5, 0.00279806}, {23.6, 0.00279806}, {23.7, 0.00279806}, {23.8, 0.00279806}, {23.9, 0.00314781}, {24., 0.00314781}, {24.1, 0.00314781}, {24.2, 0.00314781}, {24.3, 0.00314781}, {24.4, 0.00349757}, {24.5, 0.00349757}, {24.6, 0.00349757}, {24.7, 0.00349757}, {24.8, 0.00349757}, {24.9, 0.00349757}}

The model expression is:

vdc(RR_) :=
Module({R = RR, (Alpha)0 = 16.36606, (Alpha)1 =
0.70172, (Beta)0 = 17.19338, (Beta)1 = 0.09574},
A(num1_Integer) := (Alpha)0 num1^-(Alpha)1;
B(num2_Integer) := (Beta)0 Exp(-(Beta)1*num2);
(Rho) = 5.5 + 1.25 R0;
R0 = 6.46;
(Chi)(R_,
n_Integer) := (1 - Exp(-(A(n)*R/(Rho)) - (B(n)*R^2/(Rho)^2)))^n;
v(R_) = Sum((Chi)(R, n)*Subscript(CC, n)/R^n, {n, 6, 10, 2})
);
vtotal(R_) := (a* Exp(-b * R)) - vdc(R)


and here is my trial code

dataplot = ListPlot(data, PlotRange -> {{2., 8.0}, {-400.0, 700}},
ImageSize -> Large);
fit = FindFit(data,
vtotal(R), {Subscript(CC, 6), Subscript(CC, 8), Subscript(CC, 10),
a, b}, R, Method -> NMinimize, AccuracyGoal -> Infinity);

Show({pt1, Plot(vtotal(R) /. fit, {R, 2., 8.0})},
AxesLabel -> {"R",
"Potential Energy,!(*SuperscriptBox((cm), (-1)))"})
$$$$


## fitting – Finding a set of line segments to fit noisy data

Looking for a methodology to choose line segments that are a rough fit to a given set of data. In this example, the data are {x,y} pairs. For example, if the data looked like what is shown on the left, then would like to find a few line segments that go through the data, as shown on the right.

For this application

• line segments are required – curves will not work with other parts of the system
• line segments are continuous, so that the end of one line segment is the beginning of the next.
• the number of line segments is arbitrary – chosen either by the user or by an improved algorithm

A methodology that works is shown below. Any recommendations for other methods that might be more general or more efficient would be appreciated.

The methodology below uses FixedPoint and FindMinimum. At the inner level, it uses FindMinimum to determine new y-values for pairs of points, starting with the points 1 and 2, proceeding to points 2 and 3, and ending with points n-1 and n. At the outer level, the methodology below uses FixedPoint to repeat this process or stop after the maximum number of iterations is reached.
The methodology below pushes the following responsibilities to the user:

• number of points to use for the line segments
• x-value for each point
• range of x and y values (though this could easily be automated)

Seeking suggestions about other approaches or improvements to what is shown below.
Thanks!

(*problem definition*)
ptsData = {N@#,
N@((-3.5 #^2 + 3 #) Exp(3 #) ) (1 +
RandomReal({-0.075, +0.075}))} & /@  RandomReal({0, 1}, 500);
xyStart = {#, 0} & /@ {0, 0.2, 0.5, 0.6, 0.75, 0.85, 0.95, 1.0};
xRange = {0, 1};
yRange = {-20, 10};
(*analysis*)
xyNew = findNewYvaluesFromData(ptsData, xRange, yRange, xyStart, 10)
(*results*)
ListPlot( ptsData, PlotRange -> { Automatic, {-5, 5} },
Epilog -> {Orange, AbsoluteThickness(2), AbsolutePointSize(5),
Line(xyNew) , Red, Point(xyNew)})


And below is the methodology implemented thus far

Clear(findNewYvaluesFromData)
(*repeatdly improve y values in the list xyIn, until convergence or
maximum number of iterations, nIts*)
findNewYvaluesFromData(
xyData_, {xminIn_, xmaxIn_}, {yminIn_, ymaxIn_}, xyIn_, nIts_) :=
FixedPoint(
findNewYvaluesFromData(
xyData, {xminIn, xmaxIn}, {yminIn, ymaxIn}, #) &, xyIn, nIts)

(*improve y values in the list xyIn, by minimizing the deviation
between xyData and a linear interpolation of the list xyIn*)
findNewYvaluesFromData(
xyData_, {xminIn_, xmaxIn_}, {yminIn_, ymaxIn_}, xyIn_) :=
Fold(update2YvaluesFromData(
xyData, {xminIn, xmaxIn}, {yminIn, ymaxIn},  #1, #2 ) &, xyIn,
makePairsij(Range@Length@xyIn) )

Clear(update2YvaluesFromData)
(*improve y values at postions i,j in the list xyIn  *)
(*y values are improved by comparing a linear interpolation of the
list xyIn with xyData *)
(*FindMinimum is used to determine the improved y values.*)
update2YvaluesFromData(
xyData_, {xminIn_, xmaxIn_}, {yminIn_, ymaxIn_}, xyIn_, {i_, j_}) :=
Module({xyNew, r, yi, yj},
r = FindMinimum(
avgErr2YvaluesFromData(xyData, {xminIn, xmaxIn}, xyIn, {i, j},
yi, yj), {yi, xyIn((i, 2)), yminIn, ymaxIn}, {yj, xyIn((j, 2)),
yminIn, ymaxIn}, AccuracyGoal -> 2 , PrecisionGoal -> 2);

xyNew = xyIn;
xyNew((i, 2)) = yi /. r((2));
xyNew((j, 2)) = yj /. r((2));
xyNew
)

Clear(avgErr2YvaluesFromData)
(*compare xyData with a linear interpolation function  over the range
(xmin, xmax) *)
(*linear interpolation function uses xyIn with y values replaced at
positions i and j *)
avgErr2YvaluesFromData(xyData_, {xminIn_, xmaxIn_}, xyIn_, {i_, j_},
yi_?NumericQ, yj_?NumericQ) := Module({xyNew, fLin, sum, x},
xyNew = xyPairsUpdate(xyIn,  {xminIn, xmaxIn}, {i, j}, yi, yj);
fLin = Interpolation(xyNew, InterpolationOrder -> 1);
Fold(#1 + Abs(Last@#2 - fLin(First@#2 ) ) &, 0, xyData)  /
Max(1, Length@ xyData)
)

Clear(makePairsij)
(*choose adjacent pairs from a list *)
(*makePairsij(list_) := {list((#)), list((#+1))} & /@
Range(Length@list - 1)*)
makePairsij(list_) :=
ListConvolve({1, 1}, list, {-1, 1}, {}, #2 &, List)

Clear(xyPairsUpdate)
(*prepare xyV list for Interpolation function*)
(*1) ensure that there is a point at xmin and xmax*)
(*2) remove duplicates*)
xyPairsUpdate(xyV_, {xminIn_, xmaxIn_}, {i_, j_}, yi_, yj_) :=
Module({xyNew},
(*to do: remove duplicate values*)
xyNew = Sort(xyV);
xyNew = DeleteDuplicates(xyNew, Abs(First@#1 - First@#2) < 0.0001 &);
xyNew((i, 2)) = yi;
xyNew((j, 2)) = yj;
xyNew =
If(xminIn < xyNew((1, 1)),
Prepend(xyNew, {xminIn, xyNew((1, 2))}), xyNew);
xyNew =
If(xmaxIn > xyNew((-1, 1)),
Append(xyNew, {xmaxIn, xyNew((-1, 2))}), xyNew);
xyNew
)

Clear(xyPairsCheck)
(*prepare xyV list for Interpolation function*)
(*1) ensure that there is a point at xmin and xmax*)
(*2) remove duplicates*)
xyPairsCheck(xyV_, {xminIn_, xmaxIn_}, {i_, j_}) := Module({xyNew},
(*to do: remove duplicate values*)
xyNew = Sort(xyV);
xyNew = DeleteDuplicates(xyNew, Abs(First@#1 - First@#2) < 0.0001 &);
xyNew
)



## fitting – Non linear model fit for complicated function (involving complex number and definite integral)

I am trying to fit a complex function which is basically a function for anelastic relaxation. I have 3 columns in my dataset Data , first column is the temperature (which I have denoted by x here), 2nd column is the real part and 3rd column is the imaginary part:

Data = Import("C:\Users\Sonu\Downloads\datat.asc");
real = Data((All, {1, 2}));
imag = Data((All, {1, 3}));
Model = A*(1 -NIntegrate((((1/Sqrt(2*Pi*s*H))*Exp((-Log((H/d))^2/(2*s^2))))/((1 + (I*w*t*
Exp(((H)/(1.38*10^-23*x))))))), {H, 0, Infinity}));

fit = ResourceFunction("MultiNonlinearModelFit")(Rationalize({real, imag}, 0), ComplexExpand(ReIm@Model), Rationalize({{A, 1.0*10^-4}, {tw, 1.0*10^-4}, {d, 10}, {s, 0.25}},  0), {x});
fit("ParameterConfidenceIntervalTable")
`

I am getting this error

NIntegrate::inumr: The integrand E^(-(Log(H/d)^2/(2 s^2)))/(Sqrt(2 (Pi)) Sqrt(H s) (1+I E^((7.24638*10^22 H)/x) t w)) has evaluated to non-numerical values for all sampling points in the region with boundaries {{(Infinity),0.}}.

Here A is the elastic constant, d is energy, t is relaxation time, w is the resonance frequency, s is dimensionless parameter. Intial guess of A,d,t,w,s are {1e-4, 7, 1e-12, 1e8, 0.25}

## python – Empirical probability density function and fitting the distribution with an approximating polynomial and via a Bernstein approach.?

Thanks for contributing an answer to Stack Overflow!