Nothing
R/delta.error.R
In hysteresis: Tools for Modeling Rate-Dependent Hysteretic Processes and Ellipses
Defines functions
delta.error
delta.error <- function(g){
if (g$method=="nls") {
lag.se.delta <- as.vector(deltaMethod(coef(g$fit), "(-atan (((sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))/ ((sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))) - atan((-(sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))/( (sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))))/(pi*2)",vcov(g$fit))[[2]]*g$fit.statistics["period"])
area.se.delta <- deltaMethod(coef(g$fit),"sqrt((((ampx)^2+(ampy)^2+((ampx)^2-(ampy)^2)/((cos(theta))^2-(sin(theta))^2))/2) )*sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2)*pi",vcov(g$fit))[[2]]
coercion.se.delta <- deltaMethod(coef(g$fit), "ampx*sin(-atan (((sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))/ ((sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))) - atan((-(sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))/( (sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))))",vcov(g$fit) )[[2]]
retention.se.delta <- deltaMethod(coef(g$fit), "ampy*sin(-atan (((sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))/ ((sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))) - atan((-(sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))/( (sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))))",vcov(g$fit) )[[2]]
rotated.se.angle <- coef(summary(g$fit))[5,2]*180/pi
split.angle <- deltaMethod(coef(g$fit), "(atan(sqrt(ampy^2-(ampy*sin(atan (((sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))/ ((sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))) - atan((-(sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))/( (-sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta)))))^2)/ampx))*180/pi",vcov(g$fit))[[2]]
SE <- c(coef(summary(g$fit))[c(3,4),2],"retention"=retention.se.delta,"coercion"=coercion.se.delta,"area"=area.se.delta,"lag"=lag.se.delta,"split.angle"=split.angle,coef(summary(g$fit))[c(1,2),2],"rote.deg"=rotated.se.angle)}
if (g$method=="geometric") {
vmat <- (as.vector(crossprod(g$residuals)/(g$fit.statistics["n"]-5))*solve(g$fit$hessian$J))
vmat <- vmat[(g$fit.statistics["n"]+1):(g$fit.statistics["n"]+5),(g$fit.statistics["n"]+1):(g$fit.statistics["n"]+5)]
pars <- g$fit$values
names(pars) <- c("theta","semi.major","semi.minor","cx","cy")
lag.se.delta<-as.vector(deltaMethod(pars, "atan((-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/((semi.major*sin(theta)-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))*sin(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)))/(pi*2)",vmat)[[2]]*g$fit.statistics["period"])
area.se.delta<-deltaMethod(pars,"semi.major*semi.minor*pi",vmat)[[2]]
ampx.delta<-deltaMethod(pars,"sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2)",vmat)[[2]]
ampy.delta<-deltaMethod(pars,"sqrt(((semi.major*sin(theta)-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))*sin(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2+(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2)",vmat)[[2]]
coercion.se.delta <- deltaMethod(pars,"sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2)/(sqrt(1+(((semi.major*sin(theta)-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))*sin(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)))^2))",vmat )[[2]]
retention.se.delta <- deltaMethod(pars, "-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)" ,vmat)[[2]]
rotated.se.angle <- sqrt(abs(vmat[1,1]))*180/pi
split.angle = deltaMethod(pars, "atan((sqrt((sqrt(((semi.major*sin(theta)-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))*sin(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2+(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2))^2-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2))/(sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2)))",vmat)[[2]]*180/pi
SE<-c("cx"=sqrt(abs(vmat[4,4])),"cy"=sqrt(abs(vmat[5,5])),"phase.angle"=sqrt(vmat[1,1]), "retention"=retention.se.delta,"coercion"=coercion.se.delta,"area"=area.se.delta,"lag"=lag.se.delta,"split.angle"=split.angle,"ampx"=ampx.delta,"ampy"=ampy.delta,"rote.deg"=rotated.se.angle)}
else if (g$method=="lm") {
angle.SE.delta<-deltaMethod(g$fit,"atan(xy/(1-y2))*90/pi")[[2]]
cx.SE.delta<-deltaMethod(g$fit,"-(2*y2*x-xy*y)/(4*y2-xy*xy)")[[2]]
cy.SE.delta<-deltaMethod(g$fit,"-(2*y-xy*x)/(4*y2-xy*xy)")[[2]]
area.SE.delta<-deltaMethod(g$fit,"(1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))) * (1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))) *pi")[[2]]
if (g$values["rote.deg"] < 180) {
ampx.SE.delta<-deltaMethod(g$fit,"sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2)")[[2]]
ampy.SE.delta<-deltaMethod(g$fit,"sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2)")[[2]]
lag.se.delta<-as.vector(deltaMethod(g$fit,"( -atan((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))/(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))) - atan((-((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))/ (((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))))/(pi*2)")[[2]]*g$fit.statistics["period"])
coercion.se.delta<- deltaMethod(g$fit, "(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2)")[[2]]
retention.se.delta<- deltaMethod(g$fit,"(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2)")[[2]]
split.angle<-deltaMethod(g$fit, "atan(sqrt((sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))^2-((sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2))^2)/(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2)))*180/pi")[[2]]
}
else {
ampx.SE.delta<-deltaMethod(g$fit,"sqrt((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2)")[[2]]
ampy.SE.delta<-deltaMethod(g$fit,"sqrt((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2)")[[2]]
as.vector(lag.se.delta<-deltaMethod(g$fit, "( -atan((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))/( ((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))) - atan((-((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))/ (((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))))/(pi*2)")[[2]]*g$fit.statistics["period"])
coercion.se.delta<- deltaMethod(g$fit,"(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2)")[[2]]
retention.se.delta<- deltaMethod(g$fit,"(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2)")[[2]]
split.angle<-deltaMethod(g$fit, "atan(sqrt((sqrt((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))^2-((sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2))^2)/(sqrt((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2)))*180/pi")[[2]]
}
SE <- c("cx"=cx.SE.delta,"cy"=cy.SE.delta,"retention"=retention.se.delta,"coercion"=coercion.se.delta,
"area"=area.SE.delta,"lag"=lag.se.delta,"split.angle"=split.angle,"ampx"=ampx.SE.delta,
"ampy"=ampy.SE.delta,"rote.deg"=angle.SE.delta)
}
else if (g$method=="direct") {
stop("Delta errors not currently available for direct method ellipses. Use summary() to
obtain bootstrap standard errors instead.")
}
else if(g$method=="harmonic2"){
z <- g$fit[[1]]
rss <- sum(z$residuals^2)
p <- z$rank
p1 <- 1L:p
resvar1 <- rss/z$df.residual
R1 <- chol2inv(z$qr$qr[p1, p1, drop = FALSE])
se1 <- sqrt(diag(R1) * resvar1)
z2 <- g$fit[[2]]
rss <- sum(z2$residuals^2)
p <- z2$rank
p1 <- 1L:p
resvar2 <- rss/z2$df.residual
R2 <- chol2inv(z2$qr$qr[p1, p1, drop = FALSE])
se2 <- sqrt(diag(R2) * resvar2)
cov.Ta<- R1*resvar1
cov.Tb<- R2*resvar2
cov.matrix=matrix(c(cov.Ta[2,2],cov.Ta[2,3],0,0,cov.Ta[2,3],cov.Ta[3,3],0,0,0,0,cov.Tb[2,2],cov.Tb[2,3],0,0,cov.Tb[2,3],cov.Tb[3,3]),nrow=4)
b.x <- deltamethod(~sqrt(x1^2+x2^2),c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
phase.angle <- deltamethod(~ atan(x2/x1),c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
area<-deltamethod(~pi*sqrt(x1^2+x2^2)*x3,c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
coercion<-deltamethod(~sqrt(x1^2+x2^2)/sqrt(1+(x4/x3)^2),c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
lag<-deltamethod(~atan(x3/x4)/pi,c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)*g$fit.statistics["period"]/2
ampx<-deltamethod(~sqrt(x1^2+x2^2),z$coefficients[-1], cov.matrix[1:2,1:2])
ampy<-deltamethod(~sqrt(x1^2+x2^2),z2$coefficients[-1], cov.matrix[3:4,3:4])
split.angle<- deltamethod(~atan(x4/sqrt(x1^2+x2^2)), c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
theta <- deltamethod(~ (-asin(2/(tan(atan(-x1/x2)-atan(-x3/x4))*((sqrt((x1^2+x2^2+x3^2+x4^2-sqrt((x1^2+x2^2+x3^2+x4^2)^2-4*(x1^2+x2^2)*(x3^2+x4^2)*sin(atan(-x1/x2)-atan(-x3/x4))^2))/2))/(sqrt((x1^2+x2^2+x3^2+x4^2+sqrt((x1^2+x2^2+x3^2+x4^2)^2-4*(x1^2+x2^2)*(x3^2+x4^2)*sin(atan(-x1/x2)-atan(-x3/x4))^2))/2))-(sqrt((x1^2+x2^2+x3^2+x4^2+sqrt((x1^2+x2^2+x3^2+x4^2)^2-4*(x1^2+x2^2)*(x3^2+x4^2)*sin(atan(-x1/x2)-atan(-x3/x4))^2))/2))/(sqrt((x1^2+x2^2+x3^2+x4^2-sqrt((x1^2+x2^2+x3^2+x4^2)^2-4*(x1^2+x2^2)*(x3^2+x4^2)*sin(atan(-x1/x2)-atan(-x3/x4))^2))/2)))))/2)*180/pi, c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
SE <- c("b.x"=b.x,"b.y"=se2[3],"phase.angle"=phase.angle ,"cx"=se1[1],"cy"=se2[1],"retention"=se2[2],"coercion"=coercion,"area"=area,"lag"=as.vector(lag),
"split.angle"=split.angle,"ampx"=ampx,"ampy"=ampy,"rote.deg"=theta)
}
SE
}
Try the hysteresis package in your browser
Any scripts or data that you put into this service are public.
hysteresis documentation built on May 15, 2021, 1:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions delta.error
delta.error <- function(g){
if (g$method=="nls") {
lag.se.delta <- as.vector(deltaMethod(coef(g$fit), "(-atan (((sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))/ ((sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))) - atan((-(sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))/( (sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))))/(pi*2)",vcov(g$fit))[[2]]*g$fit.statistics["period"])
area.se.delta <- deltaMethod(coef(g$fit),"sqrt((((ampx)^2+(ampy)^2+((ampx)^2-(ampy)^2)/((cos(theta))^2-(sin(theta))^2))/2) )*sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2)*pi",vcov(g$fit))[[2]]
coercion.se.delta <- deltaMethod(coef(g$fit), "ampx*sin(-atan (((sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))/ ((sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))) - atan((-(sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))/( (sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))))",vcov(g$fit) )[[2]]
retention.se.delta <- deltaMethod(coef(g$fit), "ampy*sin(-atan (((sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))/ ((sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))) - atan((-(sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))/( (sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))))",vcov(g$fit) )[[2]]
rotated.se.angle <- coef(summary(g$fit))[5,2]*180/pi
split.angle <- deltaMethod(coef(g$fit), "(atan(sqrt(ampy^2-(ampy*sin(atan (((sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta))/ ((sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))) - atan((-(sqrt((ampx^2+ampy^2+(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*sin(theta))/( (-sqrt((ampx^2+ampy^2-(ampx^2-ampy^2)/((cos(theta))^2-(sin(theta))^2))/2))*cos(theta)))))^2)/ampx))*180/pi",vcov(g$fit))[[2]]
SE <- c(coef(summary(g$fit))[c(3,4),2],"retention"=retention.se.delta,"coercion"=coercion.se.delta,"area"=area.se.delta,"lag"=lag.se.delta,"split.angle"=split.angle,coef(summary(g$fit))[c(1,2),2],"rote.deg"=rotated.se.angle)}
if (g$method=="geometric") {
vmat <- (as.vector(crossprod(g$residuals)/(g$fit.statistics["n"]-5))*solve(g$fit$hessian$J))
vmat <- vmat[(g$fit.statistics["n"]+1):(g$fit.statistics["n"]+5),(g$fit.statistics["n"]+1):(g$fit.statistics["n"]+5)]
pars <- g$fit$values
names(pars) <- c("theta","semi.major","semi.minor","cx","cy")
lag.se.delta<-as.vector(deltaMethod(pars, "atan((-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/((semi.major*sin(theta)-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))*sin(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)))/(pi*2)",vmat)[[2]]*g$fit.statistics["period"])
area.se.delta<-deltaMethod(pars,"semi.major*semi.minor*pi",vmat)[[2]]
ampx.delta<-deltaMethod(pars,"sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2)",vmat)[[2]]
ampy.delta<-deltaMethod(pars,"sqrt(((semi.major*sin(theta)-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))*sin(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2+(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2)",vmat)[[2]]
coercion.se.delta <- deltaMethod(pars,"sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2)/(sqrt(1+(((semi.major*sin(theta)-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))*sin(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)))^2))",vmat )[[2]]
retention.se.delta <- deltaMethod(pars, "-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)" ,vmat)[[2]]
rotated.se.angle <- sqrt(abs(vmat[1,1]))*180/pi
split.angle = deltaMethod(pars, "atan((sqrt((sqrt(((semi.major*sin(theta)-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))*sin(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))/cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2+(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2))^2-(-semi.major*sin(theta)*cos(pi/2-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2)+semi.minor*cos(theta)*cos(-asin(semi.major*cos(theta)/sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2))+pi/2))^2))/(sqrt((semi.major*cos(theta))^2+(semi.minor*sin(theta))^2)))",vmat)[[2]]*180/pi
SE<-c("cx"=sqrt(abs(vmat[4,4])),"cy"=sqrt(abs(vmat[5,5])),"phase.angle"=sqrt(vmat[1,1]), "retention"=retention.se.delta,"coercion"=coercion.se.delta,"area"=area.se.delta,"lag"=lag.se.delta,"split.angle"=split.angle,"ampx"=ampx.delta,"ampy"=ampy.delta,"rote.deg"=rotated.se.angle)}
else if (g$method=="lm") {
angle.SE.delta<-deltaMethod(g$fit,"atan(xy/(1-y2))*90/pi")[[2]]
cx.SE.delta<-deltaMethod(g$fit,"-(2*y2*x-xy*y)/(4*y2-xy*xy)")[[2]]
cy.SE.delta<-deltaMethod(g$fit,"-(2*y-xy*x)/(4*y2-xy*xy)")[[2]]
area.SE.delta<-deltaMethod(g$fit,"(1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))) * (1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))) *pi")[[2]]
if (g$values["rote.deg"] < 180) {
ampx.SE.delta<-deltaMethod(g$fit,"sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2)")[[2]]
ampy.SE.delta<-deltaMethod(g$fit,"sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2)")[[2]]
lag.se.delta<-as.vector(deltaMethod(g$fit,"( -atan((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))/(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))) - atan((-((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))/ (((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))))/(pi*2)")[[2]]*g$fit.statistics["period"])
coercion.se.delta<- deltaMethod(g$fit, "(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2)")[[2]]
retention.se.delta<- deltaMethod(g$fit,"(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2)")[[2]]
split.angle<-deltaMethod(g$fit, "atan(sqrt((sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))^2-((sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2))^2)/(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2)))*180/pi")[[2]]
}
else {
ampx.SE.delta<-deltaMethod(g$fit,"sqrt((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2)")[[2]]
ampy.SE.delta<-deltaMethod(g$fit,"sqrt((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2)")[[2]]
as.vector(lag.se.delta<-deltaMethod(g$fit, "( -atan((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))/( ((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))) - atan((-((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))/ (((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))))/(pi*2)")[[2]]*g$fit.statistics["period"])
coercion.se.delta<- deltaMethod(g$fit,"(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2)")[[2]]
retention.se.delta<- deltaMethod(g$fit,"(sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2)")[[2]]
split.angle<-deltaMethod(g$fit, "atan(sqrt((sqrt((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))^2-((sqrt((((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2))*sin(atan(xy/(1-y2))/2-pi/2))^2)/(sqrt((((1/sqrt((1*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) - xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*cos(atan(xy/(1-y2))/2-pi/2))^2+(((1/sqrt((1*cos(atan(xy/(1-y2))/2-pi/2)*cos(atan(xy/(1-y2))/2-pi/2) + xy*cos(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2) + y2*sin(atan(xy/(1-y2))/2-pi/2)*sin(atan(xy/(1-y2))/2-pi/2)) / (1*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y2*x-xy*y)/(4*y2-xy*xy) + xy*(2*y2*x-xy*y)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) + y2*(2*y-xy*x)/(4*y2-xy*xy)*(2*y-xy*x)/(4*y2-xy*xy) - (int) ))))*sin(atan(xy/(1-y2))/2-pi/2))^2)))*180/pi")[[2]]
}
SE <- c("cx"=cx.SE.delta,"cy"=cy.SE.delta,"retention"=retention.se.delta,"coercion"=coercion.se.delta,
"area"=area.SE.delta,"lag"=lag.se.delta,"split.angle"=split.angle,"ampx"=ampx.SE.delta,
"ampy"=ampy.SE.delta,"rote.deg"=angle.SE.delta)
}
else if (g$method=="direct") {
stop("Delta errors not currently available for direct method ellipses. Use summary() to
obtain bootstrap standard errors instead.")
}
else if(g$method=="harmonic2"){
z <- g$fit[[1]]
rss <- sum(z$residuals^2)
p <- z$rank
p1 <- 1L:p
resvar1 <- rss/z$df.residual
R1 <- chol2inv(z$qr$qr[p1, p1, drop = FALSE])
se1 <- sqrt(diag(R1) * resvar1)
z2 <- g$fit[[2]]
rss <- sum(z2$residuals^2)
p <- z2$rank
p1 <- 1L:p
resvar2 <- rss/z2$df.residual
R2 <- chol2inv(z2$qr$qr[p1, p1, drop = FALSE])
se2 <- sqrt(diag(R2) * resvar2)
cov.Ta<- R1*resvar1
cov.Tb<- R2*resvar2
cov.matrix=matrix(c(cov.Ta[2,2],cov.Ta[2,3],0,0,cov.Ta[2,3],cov.Ta[3,3],0,0,0,0,cov.Tb[2,2],cov.Tb[2,3],0,0,cov.Tb[2,3],cov.Tb[3,3]),nrow=4)
b.x <- deltamethod(~sqrt(x1^2+x2^2),c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
phase.angle <- deltamethod(~ atan(x2/x1),c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
area<-deltamethod(~pi*sqrt(x1^2+x2^2)*x3,c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
coercion<-deltamethod(~sqrt(x1^2+x2^2)/sqrt(1+(x4/x3)^2),c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
lag<-deltamethod(~atan(x3/x4)/pi,c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)*g$fit.statistics["period"]/2
ampx<-deltamethod(~sqrt(x1^2+x2^2),z$coefficients[-1], cov.matrix[1:2,1:2])
ampy<-deltamethod(~sqrt(x1^2+x2^2),z2$coefficients[-1], cov.matrix[3:4,3:4])
split.angle<- deltamethod(~atan(x4/sqrt(x1^2+x2^2)), c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
theta <- deltamethod(~ (-asin(2/(tan(atan(-x1/x2)-atan(-x3/x4))*((sqrt((x1^2+x2^2+x3^2+x4^2-sqrt((x1^2+x2^2+x3^2+x4^2)^2-4*(x1^2+x2^2)*(x3^2+x4^2)*sin(atan(-x1/x2)-atan(-x3/x4))^2))/2))/(sqrt((x1^2+x2^2+x3^2+x4^2+sqrt((x1^2+x2^2+x3^2+x4^2)^2-4*(x1^2+x2^2)*(x3^2+x4^2)*sin(atan(-x1/x2)-atan(-x3/x4))^2))/2))-(sqrt((x1^2+x2^2+x3^2+x4^2+sqrt((x1^2+x2^2+x3^2+x4^2)^2-4*(x1^2+x2^2)*(x3^2+x4^2)*sin(atan(-x1/x2)-atan(-x3/x4))^2))/2))/(sqrt((x1^2+x2^2+x3^2+x4^2-sqrt((x1^2+x2^2+x3^2+x4^2)^2-4*(x1^2+x2^2)*(x3^2+x4^2)*sin(atan(-x1/x2)-atan(-x3/x4))^2))/2)))))/2)*180/pi, c(z$coefficients[-1],z2$coefficients[-1]), cov.matrix)
SE <- c("b.x"=b.x,"b.y"=se2[3],"phase.angle"=phase.angle ,"cx"=se1[1],"cy"=se2[1],"retention"=se2[2],"coercion"=coercion,"area"=area,"lag"=as.vector(lag),
"split.angle"=split.angle,"ampx"=ampx,"ampy"=ampy,"rote.deg"=theta)
}
SE
}
Try the hysteresis package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.