Code / Scripts

This page will feature some of the open source scripts that I created over the years for various purposes. Some of these scripts are written in python, some are visual basic subroutines for computing a custom statistic, some are macros for batch image processing in Imagej/Fiji, some were written for text mining specific to interlaced video formatting analysis in win-CLI (not the best way of going about this business but gets the job done), and others are written in Matlab.

1. Dodonov’s Distance-Weighted Estimator (VBA)

2. Listwise deletion: dynamic repositioning after deleting flagged data (Python)

3. TextMiner: Batch text mining in CLI for video interlacing determination

4. Batch Thresholding Confocal images for Fractal Pre-processing

5. Batch Save all open images in Imagej/Fiji

6. Time Lapse Image capture at specified intervals (Matlab)

7. Linear Regression (Lunge and Wing Extension data; Matlab)

____________________________________

 
 
 

1. Dodonov’s Distance-Weighted Estimator

This visual basic subroutine will implement a simple distance-weighted estimation of mean and weighted SD. It’ll interactively allow you to select a range of cells with numerical data, compute the DWE statistic and allow you to interactively store distance weighted mean and SD at specified address / cells.

Please refer to the 2011 article by Dodonov for comparison of various measures of central tendency for stability across conditions, and the rationale behind using a distance-weighted estimator:
Dodonov, Y. S., & Dodonova, Y. A. (2011). Robust measures of central tendency: weighting as a possible alternative to trimming in responsetime
data analysis. Psikhologicheskie Issledovaniya, 5(19).


Sub DWE()

' ----------------------------------
' This subroutine estimates Dodonov's distance-weighted mean and SD
' Reference: Robust measures of central tendency: weighting as a possible alternative to trimming in response-time data analysis
'            Yury S. Dodonov, Yulia A. Dodonova (Psikhologicheskie Issledovaniya. ISSN 2075-7999)
' Tarun Gupta, UVM (2016)
' -----------------------------------


    Dim rRange As Variant
    Dim sub_down As Variant
    Dim sub_up As Variant
    Dim i As Variant
    Dim Cdown_val As Variant
    Dim Cup_val As Variant
    Dim sub_Total As Variant
    Dim Weight As Variant
    Dim NumberOfCells As Variant
    Dim Out_file As Variant
    Dim Sum_weighted_vals As Variant
    Dim Weighted_val() As Variant
    Dim Weights_All() As Variant
    Dim Weighted_Mean As Variant
    Dim Weighted_Stdev As Variant
    Dim val_Sub_wtMean() As Double
    Dim Cl_count As Variant



    On Error Resume Next


        Application.DisplayAlerts = False

            'Specificy Range of cells on Excel sheet
            Set rRange = Application.InputBox(Prompt:= _
                "Please select a range with your Mouse for calulating DWE.", _
                    Title:="SPECIFY RANGE for distance-weighted estimator", Type:=8)


    On Error GoTo 0

        Application.DisplayAlerts = True


        If rRange Is Nothing Then
           Exit Sub
        Else
            NumberOfCells = rRange.Cells.Count                'Counts # of cells in the specified Range
            MsgBox "Number of Cells is " & NumberOfCells
            Cell_count = 0
            Sum_weighted_vals = 0


            For Each Cl In rRange.Cells                        ' Loop through each element in the specified range
                C_val = Cl.Value
                sub_down = 0
                sub_up = 0
                Weight = 0
                Cdown_val = 0
                Cup_val = 0
                Cell_count = Cell_count + 1
                Weighted_Mean = 0


                For i = 1 To NumberOfCells - 1              'Subloop within the main loop to estimate distance of given value in relation to all other values


                    If rRange.Cells(i + 1) <> "" Then                                 '-- check for non-empty cells below current cell
                        If IsNumeric(rRange.Cells(i + 1)) = True Then                 '-- If the cell contains a numeric value
                            If i >= Cell_count Then                                   '-- If it's positioned below the current cell
                                Cdown_val = rRange.Cells(i + 1)
                                sub_down = sub_down + Abs(C_val - Cdown_val)
                            Else
                                sub_down = sub_down + 0
                            End If
                        End If
                    End If
             
                    If i < Cell_count Then
                        If rRange.Cells(i) <> "" Then
                            If IsNumeric(rRange.Cells(i)) = True Then
                                Cup_val = rRange.Cells(i)
                                sub_up = sub_up + Abs(C_val - Cup_val)
                            End If
                        End If
                    Else
                        sub_up = sub_up + 0


                    End If
                    sub_Total = sub_down + sub_up
                                                 


                Next i


                If sub_Total <> 0 Then
                    Weight_Cl = Weight_Cl + (((NumberOfCells - 1) / (sub_Total)))
                Else
                    Weight_Cl = Weight_Cl
                End If


                ReDim Preserve Weights_All(0 To NumberOfCells - 1)
                    Weights_All(Cell_count - 1) = Weight_Cl        'Array of all Weights


                ReDim Preserve Weighted_val(0 To NumberOfCells - 1)
                    Weighted_val(Cell_count - 1) = Cl.Value * Weight_Cl        'Array of all weighted values (value * weight)

            Next Cl
            
        Weighted_Mean = Application.WorksheetFunction.Sum(Weighted_val) / Application.WorksheetFunction.Sum(Weights_All)

        Cl_count = 0

            For Each Cl In rRange.Cells
               
                ReDim Preserve val_Sub_wtMean(0 To NumberOfCells - 1)
                    val_Sub_wtMean(Cl_count) = Weights_All(Cl_count) * ((Cl.Value - Weighted_Mean) * (Cl.Value - Weighted_Mean))

                Cl_count = Cl_count + 1

            Next Cl

        Weighted_Stdev = Sqr(Application.WorksheetFunction.Sum(val_Sub_wtMean) / Application.WorksheetFunction.Sum(Weights_All))

        MsgBox "Dodonov's Distance Weighted Mean for specified range is: " & Weighted_Mean & " +/- S.D. " & Weighted_Stdev

        ' Specify cells to output Distance weighted mean and SD
        Set OutRange = Application.InputBox(Prompt:= _
                "Please select two cells for writing DWE Mean and Sdev.", _
                    Title:="SPECIFY 2 Cells for writing DWE weighted Mean and Stdev", Type:=8)


        OutRange.Cells(1) = Weighted_Mean
        OutRange.Cells(2) = Weighted_Stdev
                            
        
        End If
End Sub


 

_____________ *** _____________

 
 
 

2. Listwise deletion Module with dynamic repositioning (Python)

This python function takes two dictionaries, one containing original data and another containing missing data flagged with Python None (NoneType) as arguments, unzips the dictionary to a list of keys and a tuple of lists. Then it extracts index of data values flagged with Python None and deletes values in each list that specific index while dynamically adjusting the index for deleted rows. It returns sanitized data as a dictionary. The Github link for this file is: https://github.com/tguptaUVM/Listwise_deletion


def listwise_deletion(flagged_data, original_data):
    """
   This python function takes two dictionaries, one containing original data and another containing missing data flagged with Python None (NoneType) as arguments, unzips the dictionary to a list of keys and a tuple of lists. Then it extracts index of data values flagged with Python None and deletes values in each list that specific index while dynamically adjusting the index for deleted rows. Tarun Gupta (UVM, 2016)
    """
    
    temp_flagged_data = flagged_data.copy() # shallow copy flagged data to temp container  
    s_keys, s_values = zip(*temp_flagged_data.items())    # generate list of keys and tuple of values
    y_index = set()    
    
    # Get the index of NoneType values in a set
    for x in list(s_values):
        y_index_temp = set([idx for idx, val_ in enumerate(list(x)) if type(val_) is type(None)])
        y_index.update(y_index_temp)
    
    # Loop through index of None and dynamically update 
    ## based on updated index post deletion 
    count_indices = 0
    for each_ in y_index:
        if count_indices >= 1:
            if (each_ - count_indices) >= 0:
                each_ =  each_ - count_indices
    
        for eachlist_ in range(len(list(s_values))):
            if each_ <= len(list(s_values[eachlist_:eachlist_+1])[0]):
                del (list(s_values[eachlist_:eachlist_+1])[0][each_])
        count_indices += 1
                 

    sanitized_data = dict(zip(s_keys, s_values))
          
    # Sanity Check: is sanitized_data actually sanitized of flagged values?
    for key in sanitized_data.keys():
        for val in sanitized_data[key]:
            if val == 'None':
                print("\nError: sanitized_data still contains python None values in\
                ", key)
                break
                
    print("\n\nSanitized dataset has following entries now: \n\nkey\tlength\t\
    Original # of 0's\t\tAny Leftover Nones\n","_" * 55, "\n")
    for key, val in flagged_data.items():
        print(key, ":\t", len(sanitized_data[key]), "\t\t", original_data[key].count(0), \
        "\t\t\t", sanitized_data[key].count(None))
    
    input("review counter values... Press Enter to Continue.")
        
    return sanitized_data



 

_____________ *** _____________

 
 
 

3. TextMiner: Batch text mining in windows command line interface

This file will ask you for location of files, file type, and search query.
echo It’ll mine those files for the specified query and save results in a file in specified folder. This is most definitely not the best way of going about this business and python and R provide much better support for something like this but I wanted to try it out in CLI. Anyways, the code wrapper on this page thinks that the following code is written in bash, not cli (for which there’s no support; so ignore the color formatting for bash) The Github link for this file is: https://github.com/tguptaUVM/Win_CLI-batch_processing


@if (@a==@b) @end /*
cls
@echo off
setlocal EnableDelayedExpansion
TITLE &amp;quot;TextMiner (Tarun Gupta, UVM 2016)&amp;quot;
Color 60

set scriptlocation=%cd%

echo          -------------------------------------------------------------
echo                     TEXT MINER (Tarun GUPTA, Univ of Vermont 2016)   
echo           ----------------------------------------------------------
echo -------------------------------------------------------------

echo Current session was started on %DATE% at %TIME%
echo .            
echo . 

echo Description: This file will ask you for location of log files previously generated by ffmpeg's Idet filter for detecting interlaced video files. Input file format: .log 
echo .
echo It'll mine those files for the specified query and save results in a *.txt file in the same folder
echo .            
echo .          
echo -------------------------------------------------------------

echo Please enter location of folder containing log files you'd like to parse?:
pause

CALL :Open_Folder

CHDIR /D %folderlocation%

echo -------------------------------------------------------------
set /p searchkeyword=  Enter the keyword you'd like to mine these files for? :  
echo -------------------------------------------------------------
set /p outputfilename= Please enter the output summary file name:  
echo -------------------------------------------------------------
findstr &amp;quot;%searchkeyword%&amp;quot; *.log&amp;gt; %folderlocation%\%outputfilename%.txt
echo -------------------------------------------------------------
echo .            
echo .
echo .            
echo .
echo .            
echo .
echo .            
echo .
echo Done! Check %outputfilename% in location: %folderlocation%

CHDIR /D  %scriptlocation%
echo -------------------------------------------------------------

CALL :AllChoices
Echo All done!
GOTO :eof

:AllChoices
echo .
echo .            
echo .
echo What would you like to do now? Select options for further processing:
echo 1. Open the processed file
echo 2. Identify non-interlaced videos
echo 3. Exit

set /p Choice= Enter value: .. 
echo.%Choice:~,1%| findstr /r /i &amp;quot;[1-3]&amp;quot;
if %errorlevel% equ 0 (
if %Choice%==1 ( 
CALL :openfile1 %folderlocation% %outputfilename%
) ELSE IF %Choice%==2 ( 
CALL :InterlaceParse %folderlocation% %outputfilename%
) ELSE IF %Choice%==3 ( 
echo Tata! 
pause
Exit
)
) else (
  echo Please provide a valid choice
)

Goto :eof


:openfile1
echo Alright, here it is!
%folderlocation%\%outputfilename%.txt
GOTO :AllChoices


:InterlaceParse
set Prog_list=
echo ------------------------------------------------------------- 
echo This function will analyze first 500 frames of the video in order to identify files that meet following cretieria:
echo .
echo 1) TFF=BFF= Zero. That is, no TFF or BFF frames.
echo 2) Number of progressive frames is greater than number of undetermined frames
echo.
echo The result file will only contain file names that meet both of these criteria. 
echo .
echo .
pause

set inputlocation=%folderlocation%\%outputfilename%


for /f &amp;quot;tokens=1,13,15 delims=: &amp;quot; %%A in (
  'findstr /r /c:&amp;quot;[TFF:][ ]*0[ ]*BFF:[ ]*0&amp;quot; %inputlocation%.txt'
) do (
echo .
echo .
echo .
echo %%A %%B %%C
set progressive=%%B
set undetermined=%%C
IF %%B GTR %%C (
for %%a in (%%A) do ( 
set Prog_list=!Prog_list! %%a
) 
) ELSE (
for %%g in (%%A) do ( 
set undeter_list=!undeter_list! %%g
)
)
)
(for %%i in (%Prog_list%) do echo '%%i') &amp;gt; %inputlocation%_Progressive.txt
(for %%k in (%undeter_list%) do echo '%%k') &amp;gt; %inputlocation%_Undetermined.txt
echo .
echo .
echo ---------------------------------------------------
echo Please note that the number of progressive frames and undetermined frames could be different
echo for Single frame detection and multi frame detection for the same video file. Consequently, only 
echo the files that appear in the progressive file but not in the undetermined files can be considered
echo to contain only propgressively-acquired frames.
echo .
echo .
echo ---------------------------------------------------
echo Following files were saved:
echo .
echo 1) %inputlocation%_Progressive.txt
echo 2) %inputlocation%_Undetermined.txt
echo ---------------------------------------------------
echo .
echo .
pause
Exit

:Open_Folder
@if (@a==@b) @end /*

@echo off
setlocal

for /f &amp;quot;delims=&amp;quot; %%I in ('cscript /nologo /e:jscript &amp;quot;%~f0&amp;quot;') do (
    	echo You chose %%I
	ENDLOCAL &amp;amp; SET folderlocation= %%I
)

goto :EOF

:: JScript portion */

var shl = new ActiveXObject(&amp;quot;Shell.Application&amp;quot;);
var folder = shl.BrowseForFolder(0, &amp;quot;Please choose a folder.&amp;quot;, 0, 0x00);
WSH.Echo(folder ? folder.self.path : '');


 

_____________ *** _____________

 
 
 

4. Batch Thresholding: Fractal Pre-processing

This macro with batch convert all image files to binary outlines after Mean-based and OTSU based thresholding in Fiji/ImageJ. It requires the user to specify the folder that contains all raw image files. After that, the script creates two subfolders within the specified folder, duplicates the original raw image files for OTSU (often tracks cell morphology, dendritic branches etc better than other thresholding methods) and Mean-based (relatively more sensitive to noise – sometimes preferred if “noise” is biologically meaningful) thresholding. These files are then converted to binary outlines and adds appropriate processing reference to the file names for tracking to the raw source. This Macro is very useful for preprocessing large number of images for fractal analysis.


// This Macro will batch process all image (.tiff) files in a specified folder for OTSU and 
// Mean-based thresholding. thresholding. Copies of processed files are saved in a subdirectory 
// within the folder. Processed files will retain the original file names followed by processing-specific suffix.
// Tarun Gupta, Ph.D. University of Vermont, Burlington, VT (July 19, 2016) 


Dialog.create(&amp;amp;amp;quot;Auto Thresholding (OTSU / Mean)&amp;amp;amp;quot;);
Dialog.addMessage(&amp;amp;amp;quot;(Tarun Gupta, UVM 2016) \n \n \n This Macro will batch process all image (.tiff) files in a specified folder for OTSU \n and Mean-based thresholding. Copies of processed files are saved in a subdirectory within the folder \n Processed files will retain the original file names followed by processing-specific suffix. \n Please specifiy the folder with raw .tiff files for processing.&amp;amp;amp;quot;);
Dialog.show()


input = getDirectory(&amp;amp;amp;quot;Input directory&amp;amp;amp;quot;);
  if (input==&amp;amp;amp;quot;&amp;amp;amp;quot;)
      exit(&amp;amp;amp;quot;No input directory specified&amp;amp;amp;quot;);

Dialog.create(&amp;amp;amp;quot;File type&amp;amp;amp;quot;);
Dialog.addString(&amp;amp;amp;quot;File suffix: &amp;amp;amp;quot;, &amp;amp;amp;quot;.tif&amp;amp;amp;quot;, 10);
Dialog.show();
suffix = Dialog.getString();

processFolder(input);

function processFolder(input) 

{
	list = getFileList(input);
	for (i = 0; i &amp;amp;amp;lt; list.length; i++) {
		if(File.isDirectory(input + list[i]))
			processFolder(&amp;amp;amp;quot;&amp;amp;amp;quot; + input + list[i]);
		if(endsWith(list[i], suffix))
			processFile(input, list[i]);
	}
}


function processFile(input, file) 

{
open(file);
originalfilename=getTitle;

//Uses input image information to create subdirectory for storing processed files
outputdir=getDirectory(&amp;amp;amp;quot;image&amp;amp;amp;quot;);
outputnameOTSU=&amp;amp;amp;quot;Binary_OTSU&amp;amp;amp;quot;;  //Name of the output folder for OTSU thresholded files
outputnameMEAN=&amp;amp;amp;quot;Binary_MEAN&amp;amp;amp;quot;;  //Name of the output folder for Mean thresholded files
outputpathOTSU=outputdir+outputnameOTSU+File.separator;
outputpathMEAN=outputdir+outputnameMEAN+File.separator;
  File.makeDirectory(outputpathOTSU);
  if (!File.exists(outputpathOTSU))
      exit(&amp;amp;amp;quot;Unable to create directory : already exists&amp;amp;amp;quot;);

    File.makeDirectory(outputpathMEAN);
  if (!File.exists(outputpathMEAN))
      exit(&amp;amp;amp;quot;Unable to create directory : already exists&amp;amp;amp;quot;);

// Duplicates original file and renames the copy  
originalfilename=getTitle;
newbinaryOTSU = File.nameWithoutExtension+&amp;amp;amp;quot;binaryOTSU&amp;amp;amp;quot;;
newbinaryMEAN = File.nameWithoutExtension+&amp;amp;amp;quot;binaryMEAN&amp;amp;amp;quot;;
run(&amp;amp;amp;quot;Duplicate...&amp;amp;amp;quot;, &amp;amp;amp;quot; &amp;amp;amp;quot;);
rename(newbinaryOTSU);
selectWindow(originalfilename);
run(&amp;amp;amp;quot;Duplicate...&amp;amp;amp;quot;, &amp;amp;amp;quot; &amp;amp;amp;quot;);
rename(newbinaryMEAN);

// Processing Type: OTSU thresholding
selectWindow(newbinaryOTSU);
run(&amp;amp;amp;quot;Auto Threshold&amp;amp;amp;quot;, &amp;amp;amp;quot;method=Otsu ignore_black white&amp;amp;amp;quot;);
setOption(&amp;amp;amp;quot;BlackBackground&amp;amp;amp;quot;, true);
run(&amp;amp;amp;quot;Make Binary&amp;amp;amp;quot;);
run(&amp;amp;amp;quot;Outline&amp;amp;amp;quot;);
rename(newbinaryOTSU);

//Saves OTSU binary outline file at subdirectory created previously
saveAs(&amp;amp;amp;quot;Tiff&amp;amp;amp;quot;, outputpathOTSU+newbinaryOTSU);
print(&amp;amp;amp;quot;&amp;amp;amp;quot;);
print(&amp;amp;amp;quot;OTSU thresholded binary outline&amp;amp;amp;quot;, newbinaryOTSU, &amp;amp;amp;quot;saved to:&amp;amp;amp;quot;, outputpathOTSU);
close();

// Processing Type: MEAN thresholding
selectWindow(newbinaryMEAN);
run(&amp;amp;amp;quot;Auto Threshold&amp;amp;amp;quot;, &amp;amp;amp;quot;method=Mean ignore_black white&amp;amp;amp;quot;);
setOption(&amp;amp;amp;quot;BlackBackground&amp;amp;amp;quot;, true);
run(&amp;amp;amp;quot;Make Binary&amp;amp;amp;quot;);
run(&amp;amp;amp;quot;Outline&amp;amp;amp;quot;);
rename(newbinaryMEAN);

//Saves MEAN binary outline file at subdirectory created previously
saveAs(&amp;amp;amp;quot;Tiff&amp;amp;amp;quot;, outputpathMEAN+newbinaryMEAN);
print(&amp;amp;amp;quot;&amp;amp;amp;quot;);
print(&amp;amp;amp;quot;MEAN thresholded binary outline&amp;amp;amp;quot;, newbinaryMEAN, &amp;amp;amp;quot;saved to:&amp;amp;amp;quot;, outputpathMEAN);
close();

selectWindow(originalfilename);
close();

}

Dialog.create(&amp;amp;amp;quot;Auto Thresholding (OTSU / Mean)&amp;amp;amp;quot;);
Dialog.addMessage(&amp;amp;amp;quot;All Images thresholded and saved succesfully. \n See Log for details&amp;amp;amp;quot;);
Dialog.show()

 

_____________ *** _____________

 
 
 

5. Save all open images in Imagej/Fiji

This stupid little macro will save all open images as tiff files to a specified format. Helpful if you have a lot of images that need to be saved independently. Also helpful when used in conjunction with Macro for batch converting image files to binary outlines after Mean-based and OTSU based thresholding (see above).

// This Macro saves all open images as tiff files

Dialog.create(&amp;amp;amp;quot;Save ALL ROIs&amp;amp;amp;quot;);
Dialog.addMessage(&amp;amp;amp;quot;(Tarun Gupta, UVM 2016) \n \n \n This Macro saves all open images as tiff files. \n Please specifiy the folder to save.&amp;amp;amp;quot;);
Dialog.show();

outputdir = getDirectory(&amp;amp;amp;quot;Output directory&amp;amp;amp;quot;);
  if (outputdir==&amp;amp;amp;quot;&amp;amp;amp;quot;)
      exit(&amp;amp;amp;quot;No output directory specified&amp;amp;amp;quot;);

while (nImages&amp;amp;amp;gt;0) 
{ 
outputtitle = getTitle(); 
saveAs(&amp;amp;amp;quot;Tiff&amp;amp;amp;quot;, outputdir+outputtitle);
print(&amp;amp;amp;quot;File saved: &amp;amp;amp;quot;, outputtitle, &amp;amp;amp;quot;in folder:&amp;amp;amp;quot;, outputdir);  
close();
}

Dialog.create(&amp;amp;amp;quot;Result!&amp;amp;amp;quot;);
Dialog.addMessage(&amp;amp;amp;quot;All ROIs saved succesfully. See Log for details.&amp;amp;amp;quot;);
Dialog.show();

 

_____________ *** _____________

 
 
 

6. Time Lapse capture (Matlab)

3 frames every hour; each frame 1 minute apart; Used successfully for estimating egg laying rhythm in a 5.6cm (D) x 7.6cm (H) chamber containing 20 females and 10 male CS flies by frame subtraction method. .m files are available on request.

% Tarun Gupta, UVM; 2016 March 29
% Acquire images at regular intervals.
% (this script requires AcquireAndsave function script)

global vid
imaqhwinfo;
info=imaqhwinfo('macvideo'); 
vid = videoinput('macvideo', 1, 'YUY2_1600x1200');
vid.ReturnedColorspace = 'rgb';
triggerconfig(vid, 'immediate');

ImageTimer = timer;
set(ImageTimer,'TimerFcn', 'AcquireAndsave(Mainctr, framectr)',...
    'Startdelay', 10,'ExecutionMode', 'singleshot');

start(vid);

for Mainctr = 1:25   % 24 for 24 hr acqusition.
    MainTime = clock;
    preview(vid);
     
    for framectr = 1:3   % will capture 3 frames, each one minute apart. 
        frameTime = clock;     
        start(ImageTimer);
        elapsed = etime(clock, frameTime);
        pause(60 - elapsed); % 1 frame every minute, value = 60
    end
    
    stoppreview(vid);
    closepreview(vid);
    stop(vid);
    Mainelapsed = etime(clock, MainTime);
    pause(3600 - Mainelapsed); % 1 hour main counter, value = 3600
  
end
stop(vid);
clear all;

 

_____________ *** _____________

AcquireAndSave Function: This function is required for the time lapse imaging code to function

function AcquireAndsave (Mainctr, framectr)
vid = videoinput('macvideo', 1, 'YUY2_1600x1200');
vid.ReturnedColorspace = 'rgb';

img = getsnapshot(vid);
frameTimename = datestr(now,'HH_MM_SS_AM a');
fname = ['Image' num2str(Mainctr), '-' num2str(frameTimename),...
    '_Frame' num2str(framectr)]; 
imwrite(img, fname, 'tiff');
disp(['Just saved image ' fname]);
datetime('now','TimeZone','local','Format','eeee, MMMM d, yyyy hh:mm:ss:ms a')
end

 

_____________ *** _____________

 
 
 

7. Linear regression (Matlab)

Example data – Number of Lunges and Single Wing Extensions; Data files and *.m files are available on request.

clear;
% Defines the Lunge and Wing extension (WExt) vectors and plots
% Lunge vs. WExt.  The &amp;amp;amp;quot;xlim&amp;amp;amp;quot; option controls
% the limits used in plotting the x-axis. Change file name to import data 
%importfile('Filename.xls');
%Lunges = data(1:31,4); % you'll have to change this acc. to excel file
%WExt= log((data(1:31,6)))';  % change this acc. to excel file data
%% Import the data
[~, ~, raw] = xlsread('/Users/mpk/Desktop/121813 preak ss model/Filename.xlsx','TbhM6 and TbhM18 data.txt');
%raw = raw(2:28,3:4);   % reads 1st and 2nd column
raw = raw(2:end,3:end); % reads 3rd and 4th column)
%% Create output variable
data = reshape([raw{:}],size(raw));

%% Allocate imported array to column variable names
Lunges = data(:,1);
WExt = data(:,2)';

%% Clear temporary variables
clearvars data raw;
plot(WExt,Lunges,'bo')
xlim([10 15])
xlabel('WExt')
ylabel('Lunges')

%% Defines the design matrix X, and regresses y on X using the &amp;amp;amp;quot;\&amp;amp;amp;quot;
% function to compute the parameter estimates in &amp;amp;amp;quot;coef&amp;amp;amp;quot;.  The
% predicted y's are then calculated by multiplying X by the coef
% vector, and the regression line is plotted with &amp;amp;amp;quot;WExt,yhat&amp;amp;amp;quot; in
% the second half of the plot command.

%X = [ones(size(WExt)) WExt];

X = [ones(size(WExt))' WExt'];
fit.coef = X\Lunges;
fit.yhat = X*fit.coef;
plot(WExt,Lunges,'bo',WExt,fit.yhat)
xlim([10 15])
xlabel('Average Number of Single Wing Extensions (M8)')
ylabel('Average Total Lunges (M8)')

n = length(Lunges);                       % Calculates the sample size
fit.resid = Lunges-fit.yhat;              % Calculate the residuals
fit.mse = sum(fit.resid.^2)/(n-2);      % Computes the mean squared error
xbar = mean(WExt);                       % Computes the mean of the x's
fit.seb1 = sqrt(fit.mse/sum((WExt-xbar).^2)); % Computes the SE(beta1).
fit.seb0 = sqrt(fit.mse*(1/n + xbar^2/...    % Computes the SE(beta0).
  sum((WExt-xbar).^2)));

b0 = fit.coef(1);                      % Parameter estimate for b0
b1 = fit.coef(2);                      % Parameter estimate for b1
df = n-2;                              % Computes degrees of freedom
tstar = tinv(0.975,df);                % Computes t* critical value
b0ci = [b0-tstar*fit.seb0 b0+tstar*fit.seb0];      % CI for b0
b1ci = [b1-tstar*fit.seb1 b1+tstar*fit.seb1];      % CI for b1

out = regstats(Lunges,WExt);          % Regresses Lunges on WExt
out.tstat                              % Prints t-test statistics

xh = 10:0.1:15;                        % xh vector from 10 to 15 by 0.1
yhath = b0+b1*xh;                      % Predicted values at xh
sxx = sum((WExt-mean(WExt)).^2);     % Total WExt variance (SS WExt)
moec = tstar*sqrt(fit.mse*(1/n + ...   % Computes margins of error for 95%
       (xh-mean(WExt)).^2/sxx));      %   confidence intervals for E(yh)
moep = tstar*sqrt(fit.mse*(1 + 1/n +... % Computes margins of error for 95%
       (xh-mean(WExt)).^2/sxx));      %   prediction intervals for yh
ci_mean = [yhath-moec;yhath+moec]';    % Confidence intervals for E(yh)
pi_yh = [yhath-moep;yhath+moep]';      % Prediction intervals for yh
h = plot(WExt,Lunges,'bo',xh,yhath,... % Plots Lunges vs. WExt with fit,
         xh,ci_mean(:,1),'r--',...     %   lower confidence band for E(yh),
         xh,ci_mean(:,2),'r--',...     %   upper confidence band for E(yh),
         xh,pi_yh(:,1),'m-.',...       %   lower prediction band for yh, &amp;amp;amp;amp;
         xh,pi_yh(:,2),'m-.')          %   upper prediction band for yh
legend(h([2,3,5]),'Fitted Line', ...   % Puts legend in lower LH corner (3),
    'Confidence Bands', ...            %   with 3 text strings corresponding
    'Prediction Bands',3);             %   to 2nd, 3rd, &amp;amp;amp;amp; 5th plot types
xlim([10 15])                          % X-axis plotting limits
xlabel('Average Number of Single Wing Extensions (M8)')
ylabel('Average Total Lunges (M8)')

% Plots the confidence and prediction bands around the fitted line for
% Lunges vs. WExt at a confidence level of 95%, and plots the simultaneous
% confidence region for the parameters as well as the joint Bonferroni
% confidence region.
figure
xlab = 'Average Number of Single Wing Extensions (M8)';
ylab = 'Average Total Lunges (M8)';
Confregplot(WExt,Lunges,xlab,ylab,95);

figure
subplot(1,2,1)
bmn = normplot(WExt)
title('Normality plot for Wing Extensions (M8)')
xlabel ('notice the outliers')
subplot(1,2,2)
tsn = normplot(Lunges);
title('Normality plot for Total Lunges (M8)')

disp('T-stat:')
disp(out.tstat.t)

disp('T-stat - p-value:')
disp(out.tstat.pval)

disp('R-square (M8):')
disp(out.rsquare)

disp('F-stat: ')
disp(out.fstat)

disp('degree of freedom: ')
disp(out.tstat.dfe)

disp('   Lunges    Lunges-yhat  Lunges-residuals')
disp([Lunges fit.yhat fit.resid])


 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

_____________ *** _____________

Advertisements