hyper.sigcor {hyper.fit}

Function to convert from biased sample sigma to unbiased population sigma

Description

Calculates the required corrections for transforming the biased estimate of sigma to an unbiased estimate, and for transforming the sample expectation to the population expectation.

Usage

  hyper.sigcor(N,parmDF)

Arguments

N

The number of data points in the sample where sigma is being estimated.

parmDF

The number of degrees of freedom for the parameters. In the case of fitting a 1d Gaussian to data via maximum likelihood (equivalent the measuring the standard deviation of the data) parmDF=1, when fitting a 1d line with intrinsic scatter parmDF=2 and when fitting a plane to 3d data with intrinsic scatter parmDF=3.

Value

A vector of length 3. The first element is the correction from the biased to unbiased estimate for sigma. The second element is the correction from the sample to population estimate for sigma. The third is the combination of the previous two (i.e. the total correction the user will typically want to apply to their data).

Author(s)

Aaron Robotham and Danail Obreschkow

References

Robotham, A.S.G., & Obreschkow, D., PASA, in press

See Also

hyper.basic, hyper.convert, hyper.fit-data, hyper.fit, hyper.plot, hyper.sigcor, hyper.summary

Examples

#The below will take *a long* time to run- of the order a few days for the LD tests.
## Not run: 
Ngen=1e3
sdsamp=3
testvec=c(5,10,20,50,100)

set.seed(650)

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))}
    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
  }
  convtest2dOpt=rbind(convtest2dOpt, c(N=Nsamp,Raw=mean(fittest[,3]),
  mean(fittest[,3])*hyper.sigcor(Nsamp, 2)))
}

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))}
    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',
    algo.func='LD', algo.method='GG', Specs=list(Grid=seq(-0.1,0.1, len=5), dparm=NULL,
    CPUs=1, Packages=NULL, Dyn.libs=NULL))$parm
    print(fittest[i,])
  }
  convtest2dLD=rbind(convtest2dLD, c(N=Nsamp, Raw=mean(fittest[,3]),
  mean(fittest[,3])*hyper.sigcor(Nsamp, 2)))
}

convtest1dNorm={}
for(Nsamp in testvec){
  print(paste('Nsamp=', Nsamp))
  normtest={}
  for(i in 1:Ngen){
    if(i %% 100==0){print(paste('Ngen=',i))}
    normtemp=rnorm(Nsamp, sd=sdsamp)
    normtest=c(normtest, sqrt(sum((normtemp-mean(normtemp))^2)/Nsamp))
  }
  convtest1dNorm=rbind(convtest1dNorm, c(N=Nsamp, Raw=mean(normtest),
  mean(normtest)*hyper.sigcor(Nsamp, 1)))
}

## End(Not run)
#The runs above have been pre-generated and can be loaded via

data(convtest2dOpt)
data(convtest2dLD)
data(convtest1dNorm)

magplot(convtest2dOpt[,c('N','Raw')],xlim=c(5,100),ylim=c(0,4),type='b',log='x')
lines(convtest2dOpt[,c('N','sampbias2popunbias')],type='b',lty=2,pch=4)
lines(convtest2dLD[,c('N','Raw')],type='b',col='blue')
lines(convtest2dLD[,c('N','bias2unbias')],type='b',lty=2,pch=4,col='blue')
lines(convtest1dNorm[,c('N','Raw')],type='b',col='red')
lines(convtest1dNorm[,c('N','sampbias2popunbias')],type='b',lty=2,pch=4,col='red')
legend('topleft', legend=c('2 DoF and optim fit','2 DoF and LD fit', '1 DoF and direct SD'),
col=c('black','blue','red'),pch=1)
legend('topright', legend=c('Raw intrinsic scatter', 'Corrected intrinsic scatter'),
lty=c(1,2))