MACRO
NMODE X

#Minitab Local Macro. Call in session window of Minitab.
#AMC Software. ANALYTICAL METHODS COMMITTEE, ROYAL SOCIETY OF CHEMISTRY.
#______________________________________________________________________________
#ACTION: Plots a 128-point kernel density of univariate data in X with selected values of H.
#Also optionally produces a plot of number of modes versus H.
#PURPOSE: to assist visual judgement in producing an optimal representation of skewed, heavy tailed or multimodal data.
#METHOD: convolutes a Fourier transform of the discretised data with FT of normal distribution with a variance of H-squared.
#LAST MODIFIED: 22/02/06
#AUTHOR: M Thompson, Birkbeck College, London.
#SAVE as file called NMODE14.MAC as file-type ‘All files’.
#CALL with command %NMODE ‘DATA’ if saved in subdirectory ‘Macros’ of Minitab.
#CALL with command %<pathname><filename> if saved otherwise.

#______________________________________________________________________________
#i/o declarations
MCOLUMNS X 	#input data, unchanged by procedure

#scratch declarations
MCOLUMNS HLIST 	#list of h-values for plot
MCOLUMNS NMODE	#corresponding list of number of modes
MCOLUMNS TK 	#discretised x-axis, 128 points
MCOLUMNS KDEX	#correspoding densities
MCONSTANTS H	# smoothing parameter
MCOLUMNS XX	#discretised data

MCOLUMNS Y L SL SLGO
MCOLUMNS ZERO REAL IMAG REALA IMAGA REALB IMAGB
MCOLUMNS HVALUES RHVAL HALPHA
MCOLUMNS CTEMP


MCONSTANTS M S
MCONSTANTS BOTTOM TOP
MCONSTANTS MODES	
MCONSTANTS XMIN XMAX DELTA
MCONSTANTS NY NYA POINTER DENOM TEMP
MCONSTANTS I K KK KOUNT
MCONSTANTS NH ANSWER
MCONSTANTS ITITLE CNAME INOTE INOTE2
MCONSTANTS HAL IMAX IMIN GAP

#____________________________________________________________________________

NOTITLE
NOTE
NOTE MACRO RUNNING--PLEASE WAIT
NOTE
NOTE OUTLIERS IN THE DATA MAY BE OMITTED FROM KERNEL DENSITY PLOTS

#set up plot title
LET ITITLE="Kernel densities for column "
KKNAME CNAME X
KKCAT ITITLE CNAME ITITLE

#locate and trim data
CALL H15 X M S
LET BOTTOM=M-3.5*S
LET TOP=M+3.5*S
COPY X Y;
	USE X BOTTOM:TOP.

#warnings for small datasets
LET NY=N(Y)
NAME NY 'NUMBER ='
IF NY<10
	NOTE
	NOTE PROCEDURE NOT RECOMMENDED FOR SMALL DATASETS.
	PRINT NY
	NOTE EXITING FROM MACRO NMODE.
	TITLE
	EXIT
ELSEIF NY<20
	NOTE
	NOTE WARNING - PROCEDURE MAY BE SUSPECT FOR SMALL DATASETS.
	PRINT NY
ENDIF
	
#discretise x-axis points
LET XMIN=M-5*S
LET XMAX=M+5*S
LET DELTA=(XMAX-XMIN)/127

SET TK
XMIN:XMAX/DELTA
END

#initialising
SORT Y Y
LET POINTER=2
DO I=1:128
	LET XX(I)=0
ENDDO
LET DENOM=NY*DELTA*DELTA

# discretise the data to adjacent points 

DO I=1:NY
	LET TEMP=Y(I)

	DO K=POINTER:128
		IF TEMP >= TK(K)
			NEXT
		ELSE
			LET KK=K-1
			LET POINTER=K
			LET XX(KK)=XX(KK)+(TK(K)-TEMP)/DENOM
			LET XX(K)=XX(K)+(TEMP-TK(KK))/DENOM
			BREAK
		ENDIF
	ENDDO
ENDDO

# fourier transform of data 
SET ZERO
128(0)
END

CALL FFT XX ZERO REAL IMAG

# basis for ft of standard normal deviate

SET SLGO
1:64 64:1
END

# optional skip plot of number of modes vs h
NOTE
NOTE DO YOU WANT A PRELIMINARY PLOT OF NUMBER OF MODES VS. H? (Y/N)
NOTE MAY TAKE SEVERAL MINUTES.


YESNO ANSWER
IF ANSWER<>1
	GOTO 10
ENDIF

# initialise things for loop through required values of H 
NOTE
NOTE LONG PROCEDURE--PLEASE WAIT

ERASE HLIST
ERASE NMODE
LET KOUNT=0

# loop through required values of H

DO I=0.1:1/0.02

	IF I=0.2
		NOTE
		NOTE PROCEDURE ABOUT 20% COMPLETE
	ENDIF

	LET H=I*S
	#calc ft of normal distribution of variance H squared
	LET SL=6.28318*SLGO/(XMAX-XMIN)
	LET SL= EXPO(-0.5*(H*SL)**2)

	# convolution of data with normal distribution
	LET REALA=REAL*SL
	LET IMAGA=IMAG*SL

	# inverse fft to get kernel density estimate
	CALL FFT REALA IMAGA REALB IMAGB;
		INVERSE 1.
	
	#find corresponding number of modes
	CALL KOUNTMO REALB MODES
	LET KOUNT=KOUNT+1
	LET HLIST(KOUNT)=H
	LET NMODE(KOUNT)=MODES

	#shortcut to end when only one mode detected
	IF MODES=1
		LET KOUNT=KOUNT+1
		LET HLIST(KOUNT)=S
		LET NMODE(KOUNT)=1
		BREAK
	ENDIF

ENDDO


PLOT NMODE*HLIST;
	CONNECT;
		STEP 0;
	TITLE ITITLE;
	SCALE 1;
		MIN 0;
	SCALE 2;
		MIN 0;
		MAX 10;
	FOOTNOTE "High end of line is at the";
	FOOTNOTE "robust standard deviation of data.";
	AXLABEL 1 "Value of h";
	AXLABEL 2 "Number of Modes".

#optional plots of kernels
NOTE
NOTE DO YOU WANT PLOTS OF KERNEL DENSITIES? (Y/N)

YESNO ANSWER
IF ANSWER<>1
	TITLE
	EXIT
ENDIF

#input selected h-values
MLABEL 10
NOTE
NOTE ENTER H-VALUES FOR YOUR SELECTED KERNEL DENSITIES.
NOTE TYPE SPACE DELIMITED NUMBERS, THEN PRESS <ENTER>.
NOTE THEN TYPE <END> AND PRESS <ENTER>.

SET HVALUES;
FILE 'TERMINAL'.

NOTE
NOTE MACRO RUNNING--PLEASE WAIT

# compose footnoteS with h-values and n-data information
LET NH=N(HVALUES)
SORT HVALUES HVALUES
LET RHVAL=ROUND(HVALUES*1000)/1000
NTOA RHVAL HALPHA
LET INOTE="H-values are: "
LET GAP="  "
DO I=1:NH
	LET HAL=HALPHA(I)
	KKCAT INOTE HAL INOTE
	KKCAT INOTE GAP INOTE
ENDDO

LET CTEMP(1)=NY
NTOA CTEMP CTEMP
LET INOTE2="Number of data = "
LET NYA=CTEMP(1)
KKCAT INOTE2 NYA INOTE2

BRIEF 0
LAYOUT;
	TITLE ITITLE.

#loop through selected H-values and plot kernel density
DO I=1:NH
	LET H=HVALUES(I)

	#calc ft of normal distribution of variance H squared
	LET SL=6.28318*SLGO/(XMAX-XMIN)
	LET SL= EXPO(-0.5*(H*SL)**2)

	# convolution of data with normal distribution
	LET REALA=REAL*SL
	LET IMAGA=IMAG*SL

	# inverse fft to get kernel density estimate
	CALL FFT REALA IMAGA REALB IMAGB;
		INVERSE 1.

	IF I=1
		LET IMAX=1.2*MAXIMUM(REALB)
		LET IMIN=-IMAX/20
	
		PLOT REALB*TK;
		CONNECT;
		NODTITLE;
		DATA 0.2 0.9 0.2 0.9;
		FOOTNOTE INOTE;
		FOOTNOTE INOTE2;
		SCALE 2;	
			MIN IMIN;
			MAX IMAX;
		AXLABEL 1 "Variable";
		AXLABEL 2 "Density".
	ELSE
		PLOT REALB*TK;
		CONNECT;
		DATA 0.2 0.9 0.2 0.9;
			TYPE 0;
		NODTITLE;
		FOOTNOTE INOTE;
		FOOTNOTE INOTE2;
		SCALE 2;	
			MIN IMIN;
			MAX IMAX;
		AXLABEL 1 " ";
		AXLABEL 2 " ".
	ENDIF


ENDDO
ENDLAYOUT
BRIEF 2

NOTE
NOTE DO YOU WANT MORE KERNEL DENSITY PLOTS? (Y/N)
YESNO ANSWER
IF ANSWER=1
	GOTO 10
ENDIF

ENDMACRO
#__________________________________________________________________________






MACRO
FFT XX YY XTRAN YTRAN;
INVERSE FLAG.

#Minitab Local Macro. Call in session window of Mintab.
#AMC Software. ANALYTICAL METHODS COMMITTEE, ROYAL SOCIETY OF CHEMISTRY.
#__________________________________________________________________________

#ACTION: Performs fast Fourier transform (forward or inverse) of lists of complex numbers XX+i*YY.
#Transformed numbers appear in columns XTRAN YTRAN.
#CALL without subcommand for forward transform.
#CALL with subcommand INVERSE 1 for inverse FFT:
#LAST MODIFIED: 22/02/06
#AUTHOR: M Thompson, Birkbeck College, London.
#__________________________________________________________________________

#i/o declarations
MCOLUMNS XX YY 	# input data, unchanged by procedure. XX is real part, YY imaginary part.
		# these columns MUST be of length of an integer power of 2.

MCONSTANTS FLAG	# Default value of -1 for Forward FFT. Set FLAG to 1 with subcommand for inverse FFT.

MCOLUMNS XTRAN YTRAN	# output data (Fourier transformed). XTRAN is real part, YTRAN is imaginary part.

#scratch declarations
MCOLUMNS X Y Z
MCONSTANTS KOUNT M N1 N2 N3
MCONSTANTS I J K JJ KK
MCONSTANTS ANGLE ARG C S
MCONSTANTS XT YT TEMP 

DEFAULT FLAG=-1
#___________________________________________________________________________

COPY XX X
COPY YY Y

LET KOUNT=N(X)
LET M=LOGE(KOUNT)/LOGE(2)
LET N2=KOUNT

DO K=1:M
	LET N1=N2
	LET N2=N2/2
	LET N3=N2-1
	LET ANGLE=0
	LET ARG=6.28318/N1

	DO J=0:N3
		LET C=COS(ANGLE)
		LET S=FLAG*SIN(ANGLE)
		LET JJ=J+1

			DO I=JJ:KOUNT/N1
				LET KK=I+N2
				LET XT=X(I)-X(KK)
				LET YT=Y(I)-Y(KK)
				LET X(I)=X(I)+X(KK)
				LET Y(I)=Y(I)+Y(KK)
				LET X(KK)=XT*C-YT*S
				LET Y(KK)=YT*C+XT*S
			ENDDO

			LET ANGLE=(J+1)*ARG
		ENDDO
ENDDO

CALL BITREV KOUNT Z

DO I=1:KOUNT
	LET TEMP=Z(I)+1
	LET XTRAN(I)=X(TEMP)
	LET YTRAN(I)=Y(TEMP)
ENDDO

IF FLAG=1
	LET XTRAN=XTRAN/KOUNT
	LET YTRAN=YTRAN/KOUNT
ENDIF

ENDMACRO
#_____________________________________________________________________________





MACRO
BITREV N X

#Minitab Local Macro. Call in session window of Mintab.
#AMC Software. ANALYTICAL METHODS COMMITTEE, ROYAL SOCIETY OF CHEMISTRY.
#_____________________________________________________________________________

#PURPOSE: produces decimal equivalent of bit-reversed binary numbers 0:n-1.
#n MUST be an integer power of 2.
#LAST MODIFIED: 22/02/06
#AUTHOR: M Thompson, Birkbeck College, London.

MCOLUMNS X TEMP
MCONSTANTS I M MM N

LET M=LOGE(N)/LOGE(2)
LET MM=M-1

ERASE X
LET X(1)=0
LET X(2)=1

DO I=1:MM
	LET X=X*2
	LET TEMP=X+1
	STACK X TEMP X
ENDDO

ENDMACRO
#___________________________________________________________________________





MACRO
KOUNTMO Y NMODE

#Minitab Local Macro. Call in session window of Mintab.
#AMC Software. ANALYTICAL METHODS COMMITTEE, ROYAL SOCIETY OF CHEMISTRY.
#____________________________________________________________________________

#PURPOSE: counts modes in a pseudo-smooth arbitrary function
#LAST MODIFIED: 22/02/06
#AUTHOR: M Thompson, Birkbeck College, London.
#____________________________________________________________________________

#I/O
MCOLUMNS Y # list of Y-values in order of increasing X 
MCONSTANTS NMODE	# number of modes in Y

#Other declarations
MCONSTANTS I  
MCONSTANTS NOBS 
MCONSTANTS A  B  C 
MCOLUMNS Z 
#____________________________________________________________________________

DIFFERENCE Y Z
LET A = Z(2)
LET B = Z(3)
LET NMODE = 0
LET NOBS=N(y)

DO I = 4:NOBS
	LET C = Z(I)
	IF A>0
		IF B<0
			LET NMODE=NMODE+1			
		ELSEIF B=0 AND C<0
			LET NMODE=NMODE+1			
		ENDIF
	ENDIF
	LET A=B
	LET B=C
ENDDO
ENDMACRO
#____________________________________________________________________________






MACRO  
H15 X M S  
#Minitab Local Macro. Call in session window of Mintab.
#AMC Software. ANALYTICAL METHODS COMMITTEE, ROYAL SOCIETY OF CHEMISTRY.
#____________________________________________________________________________
 
#Takes input column X and returns Huber's H15 robust mean and standard deviation.   
#The input column is not changed. Follows the AMC algorithm. 
#Author: M Thompson, Birkbeck College 
#Last modified:  04/02/00 
#____________________________________________________________________________ 
 
#i/o arguments 
MCOLUMN X 	#column of input data 
MCONSTANTS M S 	#robust mean and sd 

#other declarations  
MCOLUMNS DIFF Y	#scratch columns 
MCONSTANT SA 	#scratch constant 
MCONSTANT MAD 	#median absolute difference 
MCONSTANT TOL 	#terminates iteration if < 0.00001 
MCONSTANT HIGH LOW  #limits for winsorisation 
MCONSTANT K 	#counter 
MCONSTANT MAX MIN 
#____________________________________________________________________________

# set up starting values  
IF N(X)< 6			 
	NOTE
	NOTE TOO FEW DATA. EXITING FROM MACRO H15.
	EXIT  
ELSE					 
	LET M = MEDIAN(X) 
	LET DIFF = X - M 
	LET MAD= MEDIAN(ABSOLUTE(DIFF)) 
	LET S = MAD/0.6745 
 
	LET SA = ABSOLUTE(S/M)
	IF SA < 0.000001 
		NOTE STANDARD DEVIATION VERY LOW. 
		NOTE THERE MAY BE PROBLEMS WITH DATA - CHECK! 
		NOTE EXITING FROM MACRO H15. 
		EXIT
	ENDIF 
ENDIF 
 
# Avoid H15 for n<11 or flag 
IF N(X) < 11			 
	RETURN 
ENDIF 
 
# initialise for iteration	 
LET TOL = 0.1 
LET K=0 
LET MAX=MAXIMUM(X)+1 
LET MIN=MINIMUM(X)-1 
 
# iteration 
WHILE TOL>0.0001 AND K<100		 
 
	COPY X Y  
	LET SA=S 
	LET LOW=M-1.5*S 
	LET HIGH=M+1.5*S 
 
	#winsorise raw data 
	CODE (MIN:LOW)LOW  (HIGH:MAX)HIGH Y Y 
 
	LET M=MEAN(Y) 
	LET S=STDEV(Y)/0.882 
 
	IF SA>0 
		LET TOL=ABSOLUTE(SA-S)/SA 
	ELSE
		NOTE 
		NOTE PROBLEM WITH NEAR ZERO SD IN MACRO H15
	ENDIF 
 
	LET K=K+1 
ENDWHILE 
 
IF  K>99 
	NOTE WARNING - ITERATION NOT CONVERGED IN MACRO H15
ENDIF 
 
RETURN 
ENDMACRO 
#________________________________________________________________________