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
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.