From 64a31afe4c4ae8a7323bac8256587b5fff4a1ea1 Mon Sep 17 00:00:00 2001 From: Aaron Robotham Date: Thu, 6 Jun 2024 12:04:23 +0800 Subject: [PATCH] Various fixes, as suggested by Beni Altmann (CRAN). --- DESCRIPTION | 10 +++---- R/hyper.fit.R | 6 ++-- inst/CITATION | 3 +- man/hyper.basic.Rd | 6 ++-- man/hyper.convert.Rd | 2 +- man/{hyper.data.Rd => hyper.fit-data.Rd} | 6 ++-- man/hyper.fit-package.Rd | 6 ++-- man/hyper.fit.Rd | 37 +++++++++++++----------- man/hyper.like.Rd | 2 +- man/hyper.plot.Rd | 16 +++++----- man/hyper.sigcor.Rd | 14 ++++----- man/hyper.summary.Rd | 2 +- 12 files changed, 56 insertions(+), 54 deletions(-) rename man/{hyper.data.Rd => hyper.fit-data.Rd} (95%) diff --git a/DESCRIPTION b/DESCRIPTION index 647e5b8..3beec64 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,10 +1,10 @@ Package: hyper.fit Type: Package -Title: Generic N-Dimensional Hyperplane Fitting with Heteroscedastic Covariant Errors and Intrinsic Scatter -Version: 1.1.3 -Date: 2020-10-06 -Author: Aaron Robotham and Danail Obreschkow +Title: N-Dimensional Hyperplane Fitting with Heteroscedastic Errors +Version: 1.2.0 +Date: 2024-06-06 +Author: c(person("Aaron", "Robotham", role = c("cre", "aut","cph"), email = "aaron.robotham@uwa.edu.au"), person("Danail", "Obreschkow", role = c("aut", "ctb"))) Maintainer: Aaron Robotham -Description: Includes two main high level codes for hyperplane fitting (hyper.fit) and visualising (hyper.plot2d / hyper.plot3d). In simple terms this allows the user to produce robust 1D linear fits for 2D x vs y type data, and robust 2D plane fits to 3D x vs y vs z type data. This hyperplane fitting works generically for any N-1 hyperplane model being fit to a N dimension dataset. All fits include intrinsic scatter in the generative model orthogonal to the hyperplane. +Description: High level functions for hyperplane fitting (hyper.fit()) and visualising (hyper.plot2d() / hyper.plot3d()). In simple terms this allows the user to produce robust 1D linear fits for 2D x vs y type data, and robust 2D plane fits to 3D x vs y vs z type data. This hyperplane fitting works generically for any N-1 hyperplane model being fit to a N dimension dataset. All fits include intrinsic scatter in the generative model orthogonal to the hyperplane. License: GPL-3 Depends: R (>= 3.00), magicaxis, MASS, rgl, LaplacesDemon diff --git a/R/hyper.fit.R b/R/hyper.fit.R index cfa9761..cc3d302 100644 --- a/R/hyper.fit.R +++ b/R/hyper.fit.R @@ -47,10 +47,10 @@ hyper.fit=function(X,covarray,vars,parm,parm.coord,parm.beta,parm.scat,parm.erro } if(missing(parm) & missing(parm.coord) & missing(parm.beta) & missing(parm.scat)){ - makeformula=function(objname='X',dims,vert.axis,env=.GlobalEnv){ + makeformula=function(objname='X',dims,vert.axis){ usedims=(1:dims)[-vert.axis] text=paste(objname,'[,',vert.axis,']~',paste(paste(objname,'[,',usedims,']',sep=''),collapse ='+'),sep='') - return=list(text=text,form=as.formula(text,env=env)) + return=list(text=text,form=as.formula(text)) } if(any(covarray[vert.axis,vert.axis,]>0)){ startweights=1/sqrt(covarray[vert.axis,vert.axis,]) @@ -58,7 +58,7 @@ hyper.fit=function(X,covarray,vars,parm,parm.coord,parm.beta,parm.scat,parm.erro }else{ startweights=rep(1,N) } - startfit=lm(makeformula('X',dims=dims,vert.axis=vert.axis,env=environment())$form,weights=startweights*weights) + startfit=lm(makeformula('X',dims=dims,vert.axis=vert.axis)$form,weights=startweights*weights) startfit$coef[is.na(startfit$coef)]=0 start.alphas=startfit$coef[2:dims] start.beta.vert=startfit$coef[1] diff --git a/inst/CITATION b/inst/CITATION index 3542c7b..7954aba 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -10,4 +10,5 @@ Data with Multivariate Gaussian Uncertainties", journal = "PASA", volume = "32", pages = "33", - url = "http://hyperfit.icrar.org/") \ No newline at end of file + url = "https://hyperfit.icrar.org/") + \ No newline at end of file diff --git a/man/hyper.basic.Rd b/man/hyper.basic.Rd index 72c870a..214db5b 100644 --- a/man/hyper.basic.Rd +++ b/man/hyper.basic.Rd @@ -175,7 +175,7 @@ Robotham, A.S.G., & Obreschkow, D., PASA, in press Aaron Robotham and Danail Obreschkow } \seealso{ - \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} + \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \examples{ extheta=30 #Experiment with changing this @@ -247,9 +247,7 @@ print(arrayvecmult(array(1:24,c(2,3,4)),1:3)) #Note that the following is not allowed: -\dontrun{ -print(arrayvecmult(array(1:24,c(3,2,4)),1:3)) -} +try(arrayvecmult(array(1:24,c(3,2,4)),1:3)) } % Add one or more standard keywords, see file 'KEYWORDS' in the diff --git a/man/hyper.convert.Rd b/man/hyper.convert.Rd index 5b8353b..ec77976 100644 --- a/man/hyper.convert.Rd +++ b/man/hyper.convert.Rd @@ -110,7 +110,7 @@ Robotham, A.S.G., & Obreschkow, D., PASA, in press Aaron Robotham and Danail Obreschkow } \seealso{ - \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} + \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \examples{ diff --git a/man/hyper.data.Rd b/man/hyper.fit-data.Rd similarity index 95% rename from man/hyper.data.Rd rename to man/hyper.fit-data.Rd index e3e2517..f4c3c21 100644 --- a/man/hyper.data.Rd +++ b/man/hyper.fit-data.Rd @@ -1,5 +1,5 @@ -\name{hyper.data} -\alias{hyper.data} +\name{hyper.fit-data} +\alias{hyper.fit-data} \alias{hogg} \alias{intrin} \alias{trumpet} @@ -67,7 +67,7 @@ Obreschkow & Glazebrook, 2014, ApJ, 784, 26 Aaron Robotham and Danail Obreschkow } \seealso{ - \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} + \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \examples{ hogg=read.table(system.file('data/hogg.tab', package='hyper.fit'),header=TRUE) diff --git a/man/hyper.fit-package.Rd b/man/hyper.fit-package.Rd index 7c762eb..4127fd9 100644 --- a/man/hyper.fit-package.Rd +++ b/man/hyper.fit-package.Rd @@ -11,8 +11,8 @@ \tabular{ll}{ Package: \tab hyper.fit\cr Type: \tab Package\cr -Version: \tab 1.1.3\cr -Date: \tab 2020-10-06\cr +Version: \tab 1.2.0\cr +Date: \tab 2024-06-06\cr License: \tab GPL-3\cr } Most users will interact with this package through the high level hyper.fit and hyper.plot functions. The utility functions to construct covariance matrices and arrays will also be useful (see makecovmat2d and makecovarray2d). @@ -26,7 +26,7 @@ Maintainer: Aaron Robotham Robotham, A.S.G., & Obreschkow, D., PASA, 2015, 32, 33 } \seealso{ - \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} + \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \keyword{ package } \keyword{ fit } diff --git a/man/hyper.fit.Rd b/man/hyper.fit.Rd index 81bb4ca..654befb 100644 --- a/man/hyper.fit.Rd +++ b/man/hyper.fit.Rd @@ -145,7 +145,7 @@ Robotham, A.S.G., & Obreschkow, D., PASA, in press Aaron Robotham and Danail Obreschkow } \seealso{ - \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} + \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \examples{ #### A very simple 2D example #### @@ -215,7 +215,7 @@ plot(fitwitherrorandcor) #Setup the initial data: -\dontrun{ +\donttest{ set.seed(650) sampN=200 @@ -281,7 +281,7 @@ plot(fitline, parm.errorscale=fitlineerrscale$parm['errorscale'], trans=0.2, asp #Setup the initial data: -\dontrun{ +\donttest{ set.seed(650) @@ -359,7 +359,7 @@ plot(fitlinenoexp, trans=0.2, asp=1) #Setup the initial data: -\dontrun{ +\donttest{ set.seed(650) sampN=200 @@ -427,7 +427,7 @@ plot(fitplane) #Load data: -\dontrun{ +\donttest{ data(hogg) @@ -469,7 +469,7 @@ abline(v=quantile(fithoggtrimMCMC$fit$Posterior2[,3], c(0.95,0.99)), lty=2) #### Example using 'real' data with intrinsic scatter #### -\dontrun{ +\donttest{ data(intrin) @@ -487,7 +487,7 @@ plot(fitintrincor, trans=0.1, pch='.', asp=1) #### Example using flaring trumpet data #### -\dontrun{ +\donttest{ data(trumpet) fittrumpet=hyper.fit(X=cbind(trumpet$x, trumpet$y), covarray=makecovarray2d(trumpet$x_err, @@ -500,14 +500,17 @@ plot(fittrumpet, trans=0.1, pch='.', asp=1) } -\dontrun{ +\donttest{ #To find the likelihood of zero intrinsic scatter we will need to run LaplacesDemon. The -#following will take a couple of minutes to run: +#following will take a minute or so to run. Note we change the Specs for the +#default Griddy Gibbs sampler. set.seed(650) -fittrumpetMCMC=hyper.fit(X=cbind(trumpet$x, trumpet$y), covarray=makecovarray2d(trumpet$x_err, -trumpet$y_err, trumpet$corxy), coord.type='normvec', algo.func='LD', itermax=1e5) +fittrumpetMCMC = hyper.fit(X=cbind(trumpet$x, trumpet$y), covarray=makecovarray2d( + trumpet$x_err, trumpet$y_err, trumpet$corxy), coord.type='normvec', algo.func='LD', + itermax=1e4, Specs = list(Grid=seq(-0.001,0.001, len=5), dparm=NULL, CPUs=1, + Packages=NULL, Dyn.libs=NULL)) #Assuming the user has specified the same initial seed we should find that the data #has exactly zero intrinsic scatter with ~47\% likelihood: @@ -519,7 +522,7 @@ print(fittrumpetMCMC$zeroscatprob) set.seed(650) fittrumpetMCMCerrscale=hyper.fit(X=cbind(trumpet$x, trumpet$y), covarray=makecovarray2d( -trumpet$x_err, trumpet$y_err, trumpet$corxy), itermax=1e5, coord.type='normvec', algo.func='LD', +trumpet$x_err, trumpet$y_err, trumpet$corxy), itermax=1e4, coord.type='normvec', algo.func='LD', algo.method='CHARM', Specs=list(alpha.star = 0.44), doerrorscale=TRUE) #Assuming the user has specified the same initial seed we should find that the data @@ -531,7 +534,7 @@ print(fittrumpetMCMCerrscale$zeroscatprob) #### Example using 6dFGS Fundamental Plane data #### -\dontrun{ +\donttest{ data(FP6dFGS) @@ -547,7 +550,7 @@ plot(fitFP6dFGS, doellipse=FALSE, alpha=0.5) #Next we add the censoring weights provided by C. Magoulas: -\dontrun{ +\donttest{ fitFP6dFGSw=hyper.fit(FP6dFGS[,c('logIe_J', 'logsigma', 'logRe_J')], vars=FP6dFGS[,c('logIe_J_err', 'logsigma_err', 'logRe_J_err')]^2, weights=FP6dFGS[,'weights'], @@ -565,7 +568,7 @@ in.coord.type='alpha')) #### Example using GAMA mass-size relation data #### -\dontrun{ +\donttest{ data(GAMAsmVsize) fitGAMAsmVsize=hyper.fit(GAMAsmVsize[,c('logmstar', 'rekpc')], @@ -597,7 +600,7 @@ plot(fitGAMAsmVsizelogrenoerr, doellipse=FALSE, unlog='xy') ### Example using Tully-Fisher relation data ### -\dontrun{ +\donttest{ data(TFR) TFRfit=hyper.fit(X=TFR[,c('logv','M_K')],vars=TFR[,c('logv_err','M_K_err')]^2) @@ -607,7 +610,7 @@ plot(TFRfit, xlim=c(1.7,2.5), ylim=c(-19,-26)) ### Mase-Angular Momentum-Bulge/Total ### -\dontrun{ +\donttest{ data(MJB) MJBfit=hyper.fit(X=MJB[,c('logM','logj','B.T')], covarray=makecovarray3d(MJB$logM_err, diff --git a/man/hyper.like.Rd b/man/hyper.like.Rd index efaa9b6..b6817ec 100644 --- a/man/hyper.like.Rd +++ b/man/hyper.like.Rd @@ -47,7 +47,7 @@ Robotham, A.S.G., & Obreschkow, D., PASA, in press Aaron Robotham and Danail Obreschkow } \seealso{ -\code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} +\code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \examples{ #Setup the initial data: diff --git a/man/hyper.plot.Rd b/man/hyper.plot.Rd index 48acfa4..8f2db0d 100644 --- a/man/hyper.plot.Rd +++ b/man/hyper.plot.Rd @@ -93,7 +93,7 @@ Robotham, A.S.G., & Obreschkow, D., PASA, in press Aaron Robotham and Danail Obreschkow } \seealso{ - \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} + \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \examples{ #### A very simple 2D example #### @@ -189,7 +189,7 @@ plot(fitline, trans=0.2, asp=1) #### A 3D example with fitting a plane #### -\dontrun{ +\donttest{ #Setup the initial data: @@ -258,7 +258,7 @@ plot(fitplane) #Full data -\dontrun{ +\donttest{ data(hogg) fithogg=hyper.fit(X=cbind(hogg$x, hogg$y), covarray=makecovarray2d(hogg$x_err, hogg$y_err, @@ -288,7 +288,7 @@ hoggtrim$y_err, hoggtrim$corxy), fitobj=fithogg, trans=0.2, xlim=c(0, 300), ylim #### Example using 'real' data with intrinsic scatter #### -\dontrun{ +\donttest{ data(intrin) fitintrin=hyper.fit(X=cbind(intrin$x, intrin$y), vars=cbind(intrin$x_err, @@ -302,7 +302,7 @@ plot(fitintrin, trans=0.1, pch='.', asp=1) #### Example using flaring trumpet data #### -\dontrun{ +\donttest{ data(trumpet) fittrumpet=hyper.fit(X=cbind(trumpet$x, trumpet$y), covarray=makecovarray2d(trumpet$x_err, @@ -318,7 +318,7 @@ plot(fittrumpet, trans=0.1, pch='.', asp=1) #### Example using 6dFGS Fundamental Plane data #### -\dontrun{ +\donttest{ data(FP6dFGS) fitFP6dFGSw=hyper.fit(FP6dFGS[,c('logIe_J', 'logsigma', 'logRe_J')], @@ -331,7 +331,7 @@ plot(fitFP6dFGSw, doellipse=FALSE, alpha=0.5) #### Example using GAMA mass-size relation data #### -\dontrun{ +\donttest{ data(GAMAsmVsize) fitGAMAsmVsize=hyper.fit(GAMAsmVsize[,c('logmstar', 'rekpc')], @@ -360,7 +360,7 @@ plot(fitGAMAsmVsizelogrenoerr, doellipse=FALSE, unlog='xy') ### Example using Tully-Fisher relation data ### -\dontrun{ +\donttest{ data(TFR) TFRfit=hyper.fit(X=TFR[,c('logv','M_K')],vars=TFR[,c('logv_err','M_K_err')]^2) diff --git a/man/hyper.sigcor.Rd b/man/hyper.sigcor.Rd index 3fc6486..1be5641 100644 --- a/man/hyper.sigcor.Rd +++ b/man/hyper.sigcor.Rd @@ -30,7 +30,7 @@ Aaron Robotham and Danail Obreschkow } \seealso{ - \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} + \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \examples{ #The below will take *a long* time to run- of the order a few days for the LD tests. @@ -41,12 +41,12 @@ testvec=c(5,10,20,50,100) set.seed(650) -fittestmean={} +convtest2dOpt={} for(Nsamp in testvec){ print(paste('Nsamp=', Nsamp)) fittest=matrix(0, nrow=Ngen, ncol=3) for(i in 1:Ngen){ - if(i %% 100==0){print(paste('Ngen=',i))} + if(i \%\% 100==0){print(paste('Ngen=',i))} mockx=runif(Nsamp, -100, 100) mocky=mockx*tan(45*pi/180)+rnorm(Nsamp, sd=sdsamp) fittest[i,]=hyper.fit(X=cbind(mockx,mocky), coord.type='theta', scat.type='vert.axis')$parm @@ -55,12 +55,12 @@ for(Nsamp in testvec){ mean(fittest[,3])*hyper.sigcor(Nsamp, 2))) } -fittestmeanLD={} +convtest2dLD={} for(Nsamp in testvec){ print(paste('Nsamp=', Nsamp)) fittest=matrix(0, nrow=Ngen, ncol=3) for(i in 1:Ngen){ - if(i %% 100==0){print(paste('Ngen=',i))} + if(i \%\% 100==0){print(paste('Ngen=',i))} mockx=runif(Nsamp, -100, 100) mocky=mockx*tan(45*pi/180)+rnorm(Nsamp, sd=sdsamp) fittest[i,]=hyper.fit(X=cbind(mockx,mocky), coord.type='theta', scat.type='vert.axis', @@ -72,12 +72,12 @@ for(Nsamp in testvec){ mean(fittest[,3])*hyper.sigcor(Nsamp, 2))) } -normtestmean={} +convtest1dNorm={} for(Nsamp in testvec){ print(paste('Nsamp=', Nsamp)) normtest={} for(i in 1:Ngen){ - if(i %% 100==0){print(paste('Ngen=',i))} + if(i \%\% 100==0){print(paste('Ngen=',i))} normtemp=rnorm(Nsamp, sd=sdsamp) normtest=c(normtest, sqrt(sum((normtemp-mean(normtemp))^2)/Nsamp)) } diff --git a/man/hyper.summary.Rd b/man/hyper.summary.Rd index 08136eb..5199f1a 100644 --- a/man/hyper.summary.Rd +++ b/man/hyper.summary.Rd @@ -33,7 +33,7 @@ Prints out basic summary information for hyper.fit objects output by the hyper.f Aaron Robotham and Danail Obreschkow } \seealso{ - \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} + \code{\link{hyper.basic}}, \code{\link{hyper.convert}}, \code{\link{hyper.fit-data}}, \code{\link{hyper.fit}}, \code{\link{hyper.plot}}, \code{\link{hyper.sigcor}}, \code{\link{hyper.summary}} } \examples{