# Hysteresis

## An R package for fitting rate-dependent hysteretic loops

Hysteresis loops occur when an output variable can have multiple possible values at one input value depending on the history of the system and the direction of change in the input. This package contains functions to simulate, fit, and obtain parameter values along with their standard errors (Yang and Parkhurst) from hysteresis loops of the form $$x_{t}=b_{x}cos(2t\pi/T+\phi)+c_{x}+e_{x,t}$$ $$y_{t}=b_{y}*cos(2t\pi/T+\phi)^{n}+R*sin(2t\pi/T+\phi)^{m}+c_{y}+e_{y,t}$$

where e is a random error term. These generalized transcendental equations (Lapshin) form a hysteresis loop for a given frequency or period and set of time points t=1,2…n.

The plot below uses the function mloop which simulates hysteresis loops to show the effects of choosing various odd values for n and m.

library(knitr)
library(hysteresis)
par(mfrow=c(3,3),mai=c(0,0.2,0.2,0),ann=FALSE,xaxt="n",yaxt="n",oma=c(0,0,3,0))

for (i in c(1,3,15)){
for (j in c(1,3,15)){
obj<-mloop(m=i,n=j,n.points=100,period=99)
plot(floop(obj$x,obj$y,m=i,n=j,period=99),xlim=c(-0.8,0.8),ylim=c(-0.8,0.8))
if (i==1) title(paste("n=",j,sep=""))
if (j==1) title(ylab=paste("m=",i,sep=""),line=0,cex.sub=2)
}
}
title("Hysteresis Loops for Odd Values of m and n",outer=TRUE)


It is also possible to use even values for n.

par(mfrow=c(3,3),mai=c(0,0.2,0.2,0),ann=FALSE,xaxt="n",yaxt="n",oma=c(0,0,3,0))

for (i in c(1,3,15)){
for (j in c(2,4,16)){
obj<-mloop(m=i,n=j,n.points=100,period=99)
plot(floop(obj$x,obj$y,m=i,n=j,period=99),xlim=c(-0.8,0.8),ylim=c(-0.8,0.8))
if (i==1) title(paste("n=",j,sep=""))
if (j==2) title(ylab=paste("m=",i,sep=""),line=0,cex.sub=2)
}
}
title("Hysteresis Loops for Odd Values of m and Even Values of n",outer=TRUE)


A special case is when n=1 and m=1, this makes the hysteresis loop an ellipse. The centroid of the hysteresis loop is given by cx and cy as shown in the plot below of ellipses.

obj<-mloop(cx=0,cy=0,n.points=100,period=99)
obj2<-mloop(cx=1.5,cy=0,n.points=100,period=99)
obj3<-mloop(cx=0,cy=1.5,n.points=100,period=99)
plot(obj$x,obj$y,type="l",xlim=c(-2,3),ylim=c(-2,3),xlab="x",ylab="y",col="#6600CC",main="Centroid Given by cx and cy")
points(0,0,pch=19,col="#6600CC")
text(x=0,y=0.15,"(cx=0,cy=0)",col="#6600CC")
lines(obj2$x,obj2$y,col="#00FF66")
points(1.5,0,pch=19,col="#00FF66")
text(x=1.5,y=0.15,"(cx=1.5,cy=0)",col="#00FF66")
lines(obj3$x,obj3$y,col="#FF6600")
points(0,1.5,pch=19,col="#FF6600")
text(x=0,y=1.65,"(cx=0,cy=1.5)",col="#FF6600")


The saturation points describe where the direction of the input changes sign. The distances from the central point to the saturation points are given by b.x and b.y (saturation points at ($c_{x} \pm b_{x},c_{y} \pm b_{y}$))

for (i in c(1,2,4)){
obj<-mloop(b.x=i,n.points=100,period=99)
plot(obj$x,obj$y,xlim=c(-5,10),ylim=c(-1.4,1.4),type="l",main=paste("b.x=",i,sep=""),xlab="x",ylab="y")
points(i,0.8,pch=19)
legend(i,1,legend=c("Saturation Point","x=cx+b.x","y=cy+b.y"),bty="n")
}


for (i in c(0.8,2,4)){
obj<-mloop(b.y=i,n.points=100,period=99)
plot(obj$x,obj$y,xlim=c(-1,2),ylim=c(-5,5),type="l",main=paste("b.y=",i,sep=""),xlab="x",ylab="y")
points(0.6,i,pch=19)
legend(0.6,i,legend=c("Saturation Point","x=cx+b.x","y=cy+b.y"),bty="n")
}


Retention, or the output split point R, is the vertical distance from center to upper loop trajectory.

for (i in c(1,2,4)){
obj<-mloop(retention=i,n.points=100,period=99)
plot(obj$x,obj$y,xlim=c(-1,1),ylim=c(-5,5),type="l",main=paste("retention=",i,sep=""),xlab="x",ylab="y")
segments(0,0,0,i)
text(0.3,0.5,"Retention")
}


Finally the phase.angle, $$\phi$$, changes the location of points along the loop, but does not change the form of the loop itself. When phase.angle is zero, the loop starting point is also the saturation point.

obj<-mloop(retention=0.5,n.points=100,period=99)
plot(obj$x,obj$y,type="l",xlab="x",ylab="y",main="Starting Points for Different Values of phase.angle",xlim=c(-0.6,0.8))
for (i in c(0,90,180,260)){
obj2<-mloop(phase.angle=i,retention=0.5,n.points=1,period=99)
points(obj2$x,obj2$y,pch=19,col="gold",cex=2)
points(obj2$x,obj2$y,col="gold",cex=4)
text(obj2$x+.08,obj2$y,round(i,2))
}


## Fitting Ellipses

### The Process

Hysteresis contains one method for fitting hysteresis loops given any n and m in the function floop. In the special case of an ellipse where n=1 and m=1, four methods are available in the function fel. The two-step simple harmonic regression (harmonic2) method, the default, generally produces estimates that are less biased and have lower variances than those produced by the other methods. Since the focus is on rate-dependent hysteresis, knowledge of time for the observations is required (or if unknown, times may be assumed to be equally spaced). On the other hand, if the objective is solely to fit an ellipse, observation times are not needed for the other three methods.

set.seed(24)
ellipse1 <- mel(method=2,retention=0.4,b.x=0.6,b.y=0.8,cx=0,cy=0,sd.x=0.1,sd.y=0.1,phase.angle=0,period=24,n.points=24)
#The function **mel** can be used as an alternative to **mloop** for simulating ellipses, and it is useful because it offers four different ellipse parameterizations.
model <- fel(ellipse1$x,ellipse1$y,method="harmonic2",period=24,times="equal")
#period=24 and times="equal" are used to say that 24 equally spaced points make up an ellipse.
model

## Call:
## fel(x = ellipse1$x, y = ellipse1$y, method = "harmonic2", period = 24,
##     times = "equal")
##
## Delta Method Standard Errors and 95% C.I.'s:
##             Estimates    S.E.      low     high
## b.x          0.611746 0.02348  0.56242  0.66108
## b.y          0.829397 0.03389  0.75819  0.90061
## phase.angle  0.009077 0.03838 -0.07156  0.08971
## cx          -0.013770 0.01660 -0.04865  0.02111
## cy          -0.027886 0.02397 -0.07824  0.02247
## retention    0.423527 0.03389  0.35232  0.49474
## coercion     0.278211 0.02252  0.23090  0.32553
## area         0.813959 0.07224  0.66218  0.96574
## lag          1.803394 0.13902  1.51132  2.09546
## split.angle 53.588291 0.02678 53.53202 53.64456
## ampx         0.611746 0.02348  0.56242  0.66108
## ampy         0.931276 0.03389  0.86007  1.00248
## rote.deg    57.956870 1.53618 54.72947 61.18427


In addition to the fundamental values of the model, fel also calculates a wide variety of derived parameters. Definitions for these parameters can be found using help(loop.parameters).

model$Estimates  ## b.x b.y phase.angle cx cy retention ## 0.611745822 0.829396933 0.009076523 -0.013769807 -0.027886124 0.423527491 ## coercion area lag split.angle hysteresis.x hysteresis.y ## 0.278210971 0.813958926 1.803393783 53.588291060 0.454781971 0.510645114 ## ampx ampy rote.deg semi.major semi.minor focus.x ## 0.611745822 0.931275904 57.956870279 1.088509257 0.238023858 0.563540232 ## focus.y eccentricity ## 0.900344074 0.975798963  A wide variety of functions have S3 methods for objects of class ellipsefit produced by fel. The most important of these is summary.ellipsefit which can be used to bootstrap and summarize models produced by fel. summary(model,N=10000,studentize=TRUE)  ## Summary Call: ## summary.ellipsefit(object = model, N = 10000, studentize = TRUE) ## Call for Original Fit: ## fel(x = ellipse1$x, y = ellipse1$y, method = "harmonic2", period = 24, ## times = "equal") ## Ellipse Fitting Method: ## [1] "harmonic2" ## [1] "Two step simple harmonic least squares" ## ## Bootstrapped Estimates: ## Boot.Estimate Bias Std.Error B.q0.025 B.q0.975 ## b.x 0.61101 7.370e-04 0.025318 0.56113 0.66161 ## b.y 0.82986 -4.613e-04 0.040530 0.74803 0.90862 ## phase.angle 0.00923 -1.533e-04 0.041201 -0.07001 0.09025 ## cx -0.01376 -1.166e-05 0.017812 -0.04791 0.02183 ## cy -0.02785 -3.806e-05 0.025572 -0.07379 0.02578 ## retention 0.42397 -4.384e-04 0.050119 0.32551 0.52197 ## coercion 0.27844 -2.286e-04 0.033084 0.21389 0.34285 ## area 0.81384 1.154e-04 0.101911 0.61466 1.01507 ## lag 1.80413 -7.357e-04 0.218102 1.38333 2.22993 ## split.angle 53.66241 -7.412e-02 1.756372 50.15239 57.02858 ## hysteresis.x 0.45569 -9.123e-04 0.050774 0.35576 0.55270 ## hysteresis.y 0.50877 1.876e-03 0.072502 0.37467 0.65630 ## ampx 0.61101 7.370e-04 0.025318 0.56113 0.66161 ## ampy 0.93037 9.057e-04 0.036473 0.85875 1.00120 ## rote.deg 57.96201 -5.144e-03 1.655765 54.64828 61.21009 ## rote.rad 1.01163 -8.977e-05 0.028899 0.95379 1.06832 ## semi.major 1.08727 1.243e-03 0.033664 1.02207 1.15375 ## semi.minor 0.23826 -2.352e-04 0.028867 0.18219 0.29472 ## focus.x 0.56347 6.796e-05 0.027301 0.51032 0.61747 ## focus.y 0.89986 4.798e-04 0.037850 0.82377 0.97388 ## eccentricity 0.97615 -3.483e-04 0.006198 0.96278 0.98690  Another important S3 method is for the function plot. plot(model,main="2-step Simple Harmonic Regression Ellipse Example")  In addition, S3 methods exist for fitted, print, and residuals. ### Comparison of Ellipse Estimation Methods The two most useful ellipse estimation methods implemented by fel are the ‘harmonic2’ and ‘direct’ methods. The ‘direct’ method (Flusser and Halir) fits an ellipse without requiring time information and is more stable than the other two methods in fel, ‘lm’ and ‘nls’, which are based on linear least squares and ellipse-specific nonlinear regression respectively. The ‘direct’ method does not yet have delta method standard errors available. modeldirect <- fel(ellipse1$x,ellipse1$y,method="direct",period=24,times="equal") summodel<-summary(modeldirect,N=10000,studentize=TRUE) summodel  ## Summary Call: ## summary.ellipsefit(object = modeldirect, N = 10000, studentize = TRUE) ## Call for Original Fit: ## fel(x = ellipse1$x, y = ellipse1$y, method = "direct", period = 24, ## times = "equal") ## Ellipse Fitting Method: ## [1] "direct" ## [1] "Direct specific least squares" ## ## Bootstrapped Estimates: ## Boot.Estimate Bias Std.Error B.q0.025 B.q0.975 ## b.x 0.6455 -0.013318 0.027425 0.5943 0.70208 ## b.y 0.8143 -0.044875 0.046436 0.7259 0.91089 ## cx -0.1004 0.008620 0.026720 -0.1508 -0.04656 ## cy -0.1536 0.013048 0.032489 -0.2216 -0.09479 ## retention 0.4446 0.036902 0.033999 0.3808 0.51424 ## coercion 0.3111 0.024272 0.025168 0.2644 0.36354 ## area 0.9056 0.050730 0.067653 0.7784 1.04401 ## lag 1.8962 0.239631 0.198688 1.5244 2.30648 ## split.angle 51.7255 -1.133655 1.841465 47.8951 55.17196 ## hysteresis.x 0.4792 0.051312 0.042094 0.3982 0.56373 ## hysteresis.y 0.5318 0.093991 0.079967 0.3907 0.70569 ## ampx 0.6455 -0.013318 0.027425 0.5943 0.70208 ## ampy 0.9233 -0.015553 0.034047 0.8622 0.99635 ## rote.deg 56.1994 0.579308 1.686501 52.8728 59.52510 ## rote.rad 0.9809 0.010111 0.029435 0.9228 1.03891 ## semi.major 1.0966 -0.027721 0.037715 1.0281 1.17772 ## semi.minor 0.2614 0.023362 0.021739 0.2211 0.30632 ## focus.x 0.5929 -0.028488 0.033347 0.5288 0.66050 ## focus.y 0.8870 -0.025134 0.039144 0.8147 0.96948 ## eccentricity 0.9731 -0.009261 0.008375 0.9547 0.98732  plot(modeldirect,main="Direct Ellipse Example")  The ‘direct’ method uses different fundamental parameters than the ‘harmonic2’ method. However summary results for b.x, b.y, and retention are still available from the matrix of values produced by summary.ellipsefit. summodel$values

##              Orig.Estimate   B.q0.025    B.q0.25     B.q0.5    B.q0.75
## b.x             0.63222737  0.5943469  0.6265459  0.6447086  0.6634031
## b.y             0.76946317  0.7259104  0.7827870  0.8138749  0.8448268
## cx             -0.09181282 -0.1508035 -0.1184158 -0.1011185 -0.0830374
## cy             -0.14053096 -0.2216152 -0.1739340 -0.1517622 -0.1315260
## retention       0.48150576  0.3807594  0.4213035  0.4435704  0.4665853
## coercion        0.33537594  0.2643730  0.2937040  0.3100667  0.3273819
## area            0.95636716  0.7784350  0.8584757  0.9039469  0.9506348
## lag             2.13580220  1.5244426  1.7584421  1.8915312  2.0259208
## split.angle    50.59185554 47.8951299 50.5349709 51.8113952 52.9906854
## hysteresis.x    0.53046729  0.3982164  0.4502675  0.4789588  0.5072140
## hysteresis.y    0.62576844  0.3906621  0.4755504  0.5269423  0.5815030
## ampx            0.63222737  0.5943469  0.6265459  0.6447086  0.6634031
## ampy            0.90770114  0.8622446  0.8996684  0.9216696  0.9442699
## rote.deg       56.77867835 52.8727878 55.0799932 56.2060066 57.3031095
## rote.rad        0.99097488  0.9228042  0.9613272  0.9809799  1.0001279
## semi.major      1.06888762  1.0280588  1.0709429  1.0946092  1.1200236
## semi.minor      0.28480180  0.2210529  0.2465134  0.2607114  0.2755807
## focus.x         0.56444608  0.5287686  0.5705548  0.5921731  0.6145695
## focus.y         0.86186385  0.8147444  0.8601544  0.8858419  0.9118110
## eccentricity    0.96384960  0.9546789  0.9679234  0.9739154  0.9790191
##                 B.q0.975   Std.Error   Boot.Mean         Bias Boot.Estimate
## b.x           0.70207780 0.027425186  0.61890941 -0.013317963     0.6455453
## b.y           0.91088889 0.046435743  0.72458857 -0.044874596     0.8143378
## cx           -0.04655740 0.026720401 -0.08319326  0.008619562    -0.1004324
## cy           -0.09479259 0.032489284 -0.12748254  0.013048421    -0.1535794
## retention     0.51423640 0.033999301  0.51840781  0.036902045     0.4446037
## coercion      0.36353631 0.025167675  0.35964806  0.024272120     0.3111038
## area          1.04401490 0.067652794  1.00709691  0.050729753     0.9056374
## lag           2.30647951 0.198687776  2.37543342  0.239631219     1.8961710
## split.angle  55.17196272 1.841465396 49.45820056 -1.133654981    51.7255105
## hysteresis.x  0.56372856 0.042093547  0.58177899  0.051311707     0.4791556
## hysteresis.y  0.70568551 0.079967042  0.71975951  0.093991076     0.5317774
## ampx          0.70207780 0.027425186  0.61890941 -0.013317963     0.6455453
## ampy          0.99634787 0.034047346  0.89214800 -0.015553148     0.9232543
## rote.deg     59.52509606 1.686501393 57.35798677  0.579308420    56.1993699
## rote.rad      1.03890891 0.029435002  1.00108572  0.010110839     0.9808640
## semi.major    1.17771787 0.037714645  1.04116676 -0.027720857     1.0966085
## semi.minor    0.30632344 0.021738830  0.30816401  0.023362212     0.2614396
## focus.x       0.66050124 0.033346897  0.53595842 -0.028487664     0.5929337
## focus.y       0.96947844 0.039143872  0.83673023 -0.025133624     0.8869975
## eccentricity  0.98731742 0.008374862  0.95458873 -0.009260870     0.9731105


The four plots below illustrate a comparison of the four methods for fitting an ellipse to simulated data. The data points are based on the simulated red ellipse; the fitted ellipse is in black.

set.seed(13)
par(mfrow=c(2,2))
halfellipse <- mel(method=2,cx=20,cy=25,retention=1.2,b.x=14,b.y=0.8,sd.x=1,sd.y=0.2,period=24,n.points=16,phase.angle=pi/2)
halftrueellipse <- mel(method=2,cx=20,cy=25,retention=1.2,b.x=14,b.y=0.8,phase.angle=0,period=99,n.points=100)
harmodel<-fel(halfellipse$x,halfellipse$y,method="harmonic2",period=24,times="equal")
dirmodel<-fel(halfellipse$x,halfellipse$y,method="direct",period=24,times="equal")
lmmodel<-fel(halfellipse$x,halfellipse$y,method="lm",period=24,times="equal")
nlsmodel<-fel(halfellipse$x,halfellipse$y,method="nls",period=24,times="equal",control=c(n.iter=500))
plot(harmodel,main="Harmonic2 Model",xlim=c(5,34),ylim=c(23.4,26.9))
lines(halftrueellipse$x,halftrueellipse$y,col="red")
plot(dirmodel,main="Direct Model",xlim=c(5,34),ylim=c(23.4,26.9))
lines(halftrueellipse$x,halftrueellipse$y,col="red")
plot(lmmodel,main="Linear Model",xlim=c(5,34),ylim=c(23.4,26.9))
lines(halftrueellipse$x,halftrueellipse$y,col="red")
plot(nlsmodel,main="Non-Linear Model",xlim=c(5,34),ylim=c(23.4,26.9))
lines(halftrueellipse$x,halftrueellipse$y,col="red")


## Geometric Method

The geometric ellipse method used here is based on the work of Gander, Golub and Strebel, and the code used is an R translation of the Matlab code they provided. This method directly minimizes the Euclidean distances from the ellipse through Gauss-Newton minimization. Standard errors are obtainable through either the Delta Method and bootstrapping, however the use of bootstrapping through summary.ellipsefit is discouraged on these ellipses as the geometric method is extremely computationally expensive. The “geometric” method works best when sd.x=sd.y, and it is important to check for false convergence. One way to check for false convergence is to examine the plot after fitting the data.

set.seed(101)
ellip <- mel(rote.deg=45,semi.major=5,semi.minor=3,n.points=13,sd.x=0.4,sd.y=0.4)
true.ellip <- mel(rote.deg=45,semi.major=5,semi.minor=3,n.points=100,period=100)
ellip.geometric <- fel(ellip$x,ellip$y,method="geometric")
## write.table(summodels$Boot.Estimates,"file_name.txt")  ## Fitting Hysteresis Loops The function floop can be used to fit hysteresis loops with specific values of n and m as arguments to floop. Below is an example of a hysteresis loop with n=5, m=3. loop <- mloop(n=5, m=3,sd.x=0.02,sd.y=0.02) fitloop <- floop(loop$x,loop$y,n=5, m=3,period=24,times="equal")  ## Warning in floop(loop$x, loop$y, n = 5, m = 3, period = 24, times = "equal"): ## hysteresis.x and coercion only available if m=n  fitloop$Estimates

##                n                m              b.x              b.y
##      5.000000000      3.000000000      0.599422931      0.799635941
##      phase.angle               cx               cy        retention
##      0.012909411     -0.001077111     -0.007705282      0.205051341
##         coercion             area              lag beta.split.angle
##               NA      0.289605698      0.958833524      0.000000000
##     hysteresis.x     hysteresis.y
##               NA      0.256430871

plot(fitloop,main="Fitted Hysteresis Loop")


summary(fitloop)

## Summary Call:
## summary.fittedloop(object = fitloop)
## Call for Original Fit:
## floop(x = loop$x, y = loop$y, n = 5, m = 3, times = "equal",
##     period = 24)
##
## Bootstrapped Estimates:
##                  Boot.Estimate       Bias Std.Error  B.q0.025   B.q0.975
## n                     5.000000  0.000e+00  0.000000  5.000000  5.000e+00
## m                     3.000000  0.000e+00  0.000000  3.000000  3.000e+00
## b.x                   0.599290  1.329e-04  0.005540  0.588436  6.104e-01
## b.y                   0.799588  4.843e-05  0.007572  0.784037  8.146e-01
## phase.angle           0.025782 -1.287e-02  0.009265  0.008167  4.476e-02
## cx                   -0.001263  1.860e-04  0.004127 -0.009552  6.838e-03
## cy                   -0.007670 -3.496e-05  0.003743 -0.014682 -9.934e-05
## retention             0.205157 -1.053e-04  0.007242  0.191403  2.192e-01
## coercion                    NA         NA        NA        NA         NA
## area                  0.289689 -8.349e-05  0.010596  0.269101  3.096e-01
## lag                   0.959336 -5.025e-04  0.034190  0.892984  1.030e+00
## beta.split.angle      0.000000  0.000e+00  0.000000  0.000000  0.000e+00
## hysteresis.x                NA         NA        NA        NA         NA
## hysteresis.y          0.256549 -1.183e-04  0.009543  0.238095  2.764e-01


## Acknowledgments

NIFA MRF 25-008/W-2173 Impacts of Stress Factors on Performance, Health, and Well Being of Farm Animals