Description Usage Arguments Value Author(s) References See Also Examples
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.
1 | hyper.sigcor(N,parmDF)
|
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. |
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).
Aaron Robotham and Danail Obreschkow
Robotham, A.S.G., & Obreschkow, D., PASA, in press
hyper.basic
, hyper.convert
, hyper.data
, hyper.fit
, hyper.plot
, hyper.sigcor
, hyper.summary
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | #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)
fittestmean={}
for(Nsamp in testvec){
print(paste('Nsamp=', Nsamp))
fittest=matrix(0, nrow=Ngen, ncol=3)
for(i in 1:Ngen){
if(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)))
}
fittestmeanLD={}
for(Nsamp in testvec){
print(paste('Nsamp=', Nsamp))
fittest=matrix(0, nrow=Ngen, ncol=3)
for(i in 1:Ngen){
if(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)))
}
normtestmean={}
for(Nsamp in testvec){
print(paste('Nsamp=', Nsamp))
normtest={}
for(i in 1:Ngen){
if(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))
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.