R/mcskew.R

In musto101/wilcox_R: What the package does (short line)

mcskew <-
function(z)
{
	n=length(z)
	y1=0
	y2=0
	left=0
	right=0
	q=0
	p=0
	eps=0.0000000000001
	z=-z
	xmed=pull(z,n,floor(n/2)+1)
	if (n%%2 == 0)
	{
		xmed=(xmed+pull(z,n,floor(n/2)))/2
	}
	z=z-xmed
	y=-sort(z)
	y1=y[y>-eps]
	y2=y[y<=eps]
	h1=length(y1)
	h2=length(y2)	
	left[1:h2]=1
	right[1:h2]=h1
	nl=0
	nr=h1*h2
	knew=floor(nr/2)+1
	IsFound=0
	while ((nr-nl>n) & (IsFound==0))
	{
		weight=0
		work=0
		j=1
		for (i in 1:h2)
		{
			if (left[i]<=right[i])
			{
				weight[j]=right[i]-left[i]+1
				k=left[i]+floor(weight[j]/2)
				work[j]=calwork(y1[k],y2[i],k,i,h1+1,eps)
				j=j+1
			}
		}
		trial=whimed(work,weight,j-1)
		j=1
		for (i in h2:1)
		{
			while ((j<=h1)&(calwork(y1[min(j,h1)],y2[i],j,i,h1+1,eps)>trial))
			{
				j=j+1
			}
			p[i]=j-1
		}
		j=h1
		for (i in 1:h2)
		{
			while ((j>=1)&(calwork(y1[max(j,1)],y2[i],j,i,h1+1,eps)<trial))
			{
				j=j-1
			}
			q[i]=j+1
		}
		sump=sum(p[1:h2])
		sumq=sum(q[1:h2])-h2
		if (knew<=sump)
		{
			right[1:h2]=p[1:h2]
			nr=sump
		}
		else
		{
			if (knew>sumq)
			{
				left[1:h2]=q[1:h2]
				nl=sumq
			}
			else
			{
				medc=trial
				IsFound=1
			}
		}
	}
	if (IsFound==0)
	{work=0
		j=1
		for (i in 1:h2)
		{
			if (left[i]<=right[i])
			{
				for (jj in left[i]:right[i])
				{
					work[j]=0-calwork(y1[jj],y2[i],jj,i,h1+1,eps)
					j=j+1
				}
			}
		}
		medc=0-pull(work,j-1,knew-nl)
	}
	medc
}