diff --git a/applications/lazstats/source/forms/analysis/nonparametric/chisqrunit.pas b/applications/lazstats/source/forms/analysis/nonparametric/chisqrunit.pas index 5a051c754..dae87f80b 100644 --- a/applications/lazstats/source/forms/analysis/nonparametric/chisqrunit.pas +++ b/applications/lazstats/source/forms/analysis/nonparametric/chisqrunit.pas @@ -59,6 +59,35 @@ type FRowColPropsReportFrame: TReportFrame; FCellChiSqrReportFrame: TReportFrame; + function CalcGStatistic(const AFrequencies: IntDyneMat; + const AExpected: DblDyneMat): Double; + + function CalcLikelihoodRatio(const AFrequencies: IntDyneMat; + const AExpected: DblDyneMat): Double; + + function CalcPearsonR(const AFrequencies: IntDyneMat): Double; + + procedure GetExpectedAndCellChiSqr(const AFrequencies: IntDyneMat; + out AExpected, ACellChiSqr: DblDyneMat; + out AChiSqr: Double); + + procedure GetFrequencies(const AColNoSelected: IntDyneVec; + ARowIndex, AColIndex, ADepIndex: Integer; + out ANumRows, ANumCols, ANumCases: Integer; + out AFrequencies: IntDyneMat); + + procedure GetYatesCorrection(const AFrequencies: IntDyneMat; + out AdjChiSqr: Double); + + procedure GetColProportions(const AFrequencies: IntDyneMat; + out AProportions: DblDyneMat); + + procedure GetRowProportions(const AFrequencies: IntDyneMat; + out AProportions: DblDyneMat); + + procedure GetTotalProportions(const AFrequencies: IntDyneMat; + out AProportions: DblDyneMat); + protected procedure AdjustConstraints; override; procedure Compute; override; @@ -79,7 +108,7 @@ implementation uses Math, - Utils, GridProcs; + Utils, GridProcs, MatrixUnit; { TChiSqrFrm } @@ -145,6 +174,75 @@ begin end; +function TChiSqrFrm.CalcGStatistic(const AFrequencies: IntDyneMat; + const AExpected: DblDyneMat): Double; +var + numRows, numCols: Integer; + i, j: Integer; +begin + MatSize(AExpected, numRows, numCols); + Result := 0.0; + for i := 0 to numRows-1 do + for j := 0 to numCols-1 do + if (AExpected[i, j] > 0) then + Result := Result + AFrequencies[i, j] * ln(AFrequencies[i, j] / AExpected[i, j]); + Result := 2.0 * Result; +end; + + +function TChiSqrFrm.CalcLikelihoodRatio(const AFrequencies: IntDyneMat; + const AExpected: DblDyneMat): Double; +var + numRows, numCols: Integer; + i, j: Integer; +begin + MatSize(AExpected, numRows, numCols); + Result := 0.0; + for i := 0 to numRows-1 do + for j := 0 to numCols-1 do + if (AFrequencies[i, j] > 0.0) then + Result := Result + AFrequencies[i, j] * ln(AExpected[i, j] / AFrequencies[i, j]); + Result := -2.0 * Result; +end; + + +function TChiSqrFrm.CalcPearsonR(const AFrequencies: IntDyneMat): Double; +var + numRows, numCols, numCases: Integer; + sumX, sumY: Double; + varX, varY: Double; + i, j: Integer; +begin + MatSize(AFrequencies, numRows, numCols); + dec(numRows); // Do not iterate into the totals row and column + dec(numCols); + numCases := AFrequencies[numRows, numCols]; + + SumX := 0; + SumY := 0; + VarX := 0; + VarY := 0; + + for i := 0 to numRows-1 do + sumX := sumX + ( (i+1) * AFrequencies[i, numCols] ); + for j := 0 to numCols-1 do + sumY := sumY + ( (j+1) * AFrequencies[numRows, j] ); + for i := 0 to numRows-1 do + varX := varX + ( sqr(i+1) * AFrequencies[i, numCols] ); + for j := 0 to numCols-1 do + varY := varY + ( sqr(j+1) * AFrequencies[numRows, j] ); + varX := varX - sqr(sumX) / numCases; + varY := varY - sqr(sumY) / numCases; + + Result := 0; + for i := 0 to numRows-1 do + for j := 0 to numCols-1 do + Result := Result + ((i+1)*(j+1) * AFrequencies[i, j]); + + Result := Result - (sumX * sumY / numCases); + Result := Result / sqrt(varX * varY); +end; + procedure TChiSqrFrm.ColInClick(Sender: TObject); var index: integer; @@ -179,20 +277,19 @@ var CellChi: DblDyneMat = nil; RowLabels: StrDyneVec = nil; ColLabels: StrDyneVec = nil; - yates : boolean = false; - + yates : boolean; NoSelected, NCases, NRows, NCols: Integer; - i, j, RowNo, ColNo, DepNo, MinRow, MaxRow, MinCol, MaxCol: integer; - Row, Col, FObs, df: integer; - PObs, ChiSquare, ProbChi, phi, SumX, SumY, VarX, VarY, likelihood: double; + i, j, rowNo, colNo, depNo: integer; + Row, Col, df: integer; + ChiSquare, probChi, phi: double; title : string; - AdjChiSqr, AdjProbChi, probLikelihood, G, pearsonr, MantelHaenszel, MHprob: double; + AdjChiSqr, AdjProbChi, pearsonr, G, likelihood, MantelHaenszel, prob: double; CoefCont, CramerV: double; lReport: TStrings; begin - RowNo := GetVariableIndex(OS3MainFrm.DataGrid, RowEdit.Text); - ColNo := GetVariableIndex(OS3MainFrm.DataGrid, ColEdit.Text); - DepNo := GetVariableIndex(OS3MainFrm.DataGrid, DepEdit.Text); + rowNo := GetVariableIndex(OS3MainFrm.DataGrid, RowEdit.Text); + colNo := GetVariableIndex(OS3MainFrm.DataGrid, ColEdit.Text); + depNo := GetVariableIndex(OS3MainFrm.DataGrid, DepEdit.Text); SetLength(ColNoSelected, NoVariables); ColNoSelected[0] := RowNo; @@ -203,126 +300,29 @@ begin NoSelected := 3; ColNoSelected[2] := DepNo; end; - SetLength(ColNoSelected, NoSelected); + SetLength(ColNoSelected, NoSelected); // trim length - // Get min and max of row and col numbers - MinRow := MaxInt; - MaxRow := -MinRow; - MinCol := MaxInt; - MaxCol := -MinCol; - for i := 1 to NoCases do - begin - if not GoodRecord(OS3MainFrm.DataGrid, i, ColNoSelected) then continue; - Row := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[RowNo, i]))); - Col := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ColNo, i]))); - if Row > MaxRow then MaxRow := Row; - if Row < MinRow then MinRow := Row; - if Col > MaxCol then MaxCol := Col; - if Col < MinCol then MinCol := Col; - end; - NRows := MaxRow - MinRow + 1; - NCols := MaxCol - MinCol + 1; + // Get frequencies from the grid + GetFrequencies(ColNoselected, rowNo, colNo, depNo, nRows, nCols, nCases, Freq); + df := (nRows - 1) * (nCols - 1); - // allocate and initialize - SetLength(Freq, NRows+1, NCols+1); - SetLength(Prop, NRows+1, NCols+1); - SetLength(Expected, NRows, NCols); - SetLength(CellChi, NRows, NCols); - SetLength(RowLabels, NRows+1); - SetLength(ColLabels, NCols+1); - - for i := 1 to NRows + 1 do - for j := 1 to NCols + 1 do Freq[i-1, j-1] := 0; - - // get cell data - NCases := 0; - case InputGrp.ItemIndex of - 0 : begin // count number of cases in each row and column combination - for i := 1 to NoCases do - begin - if not GoodRecord(OS3MainFrm.DataGrid, i, ColNoSelected) then continue; - NCases := NCases + 1; - Row := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[RowNo, i]))); - Col := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ColNo, i]))); - Row := Row - MinRow + 1; - Col := Col - MinCol + 1; - Freq[Row-1, Col-1] := Freq[Row-1, Col-1] + 1; - end; - end; - 1 : begin // read frequencies data from grid - for i := 1 to NoCases do - begin - if not GoodRecord(OS3MainFrm.DataGrid, i, ColNoSelected) then continue; - Row := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[RowNo, i]))); - Col := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ColNo, i]))); - Row := Row - MinRow + 1; - Col := Col - MinCol + 1; - FObs := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[DepNo, i]))); - Freq[Row-1, Col-1] := Freq[Row-1, Col-1] + FObs; - NCases := NCases + FObs; - end; - end; - 2 : begin // get no. of cases and proportions for each cell - NCases := StrToInt(NCasesEdit.Text); - for i := 1 to NoCases do - begin - if not GoodRecord(OS3MainFrm.Datagrid, i, ColNoSelected) then continue; - Row := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[RowNo, i]))); - Col := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ColNo, i]))); - Row := Row - MinRow + 1; - Col := Col - MinCol + 1; - PObs := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[DepNo, i])); - Freq[Row-1, Col-1] := Freq[Row-1, Col-1] + round(PObs * NCases); - end; - end; - end; // end case - Freq[Nrows, Ncols] := NCases; - - // Calculate expected values - // Get row totals first - for i := 1 to NRows do - for j := 1 to Ncols do - Freq[i-1, NCols] := Freq[i-1, NCols] + Freq[i-1, j-1]; - - // Get col totals next - for j := 1 to NCols do - for i := 1 to NRows do - Freq[NRows, j-1] := Freq[NRows, j-1] + Freq[i-1, j-1]; - - // Then get expected values and cell chi-squares - ChiSquare := 0.0; - AdjChisqr := 0.0; - - if (YatesChk.Checked) and (NRows = 2) and (NCols = 2) then - yates := true; - - for i := 1 to NRows do - begin - for j := 1 to Ncols do - begin - Expected[i-1, j-1] := Freq[NRows, j-1] * Freq[i-1, NCols] / NCases; - if Expected[i-1, j-1] > 0.0 then - CellChi[i-1, j-1] := sqr(Freq[i-1, j-1] - Expected[i-1, j-1]) / Expected[i-1, j-1] - else begin - ErrorMsg('Zero expected value found.'); - CellChi[i-1,j-1] := 0.0; - end; - ChiSquare := ChiSquare + CellChi[i-1, j-1]; - end; - end; - df := (Nrows - 1) * (Ncols - 1); - if yates then // 2 x 2 corrected chi-square - begin - AdjChiSqr := abs((Freq[0,0] * Freq[1,1]) - (Freq[0,1] * Freq[1,0])); - AdjChiSqr := sqr(AdjChiSqr - NCases / 2.0) * NCases; // numerator - AdjChiSqr := AdjChiSqr / (Freq[0,2] * Freq[1,2] * Freq[2,0] * Freq[2,1]); - AdjProbChi := 1.0 - ChiSquaredProb(AdjChiSqr, df); - end; + // Calculate expected values and cell chi-squares + GetExpectedAndCellChiSqr(Freq, Expected, CellChi, ChiSquare); ProbChi := 1.0 - ChiSquaredProb(ChiSquare, df); // prob. larger chi + yates := YatesChk.Checked and (nRows = 2) and (nCols = 2); + if yates then begin + GetYatesCorrection(Freq, AdjChiSqr); + AdjProbChi := 1.0 - ChiSquaredProb(AdjChiSqr, df); + end else + AdjChiSqr := 0; + + // Get row and column labels for reports + SetLength(RowLabels, NRows+1); for i := 1 to NRows do RowLabels[i-1] := Format('Row %d', [i]); RowLabels[NRows] := 'Total'; + SetLength(ColLabels, NCols+1); for j := 1 to NCols do ColLabels[j-1] := Format('Col.%d', [j]); ColLabels[NCols] := 'Total'; @@ -356,53 +356,24 @@ begin if PropsChk.Checked then begin RowColPage.TabVisible := true; - title := 'ROW PROPORTIONS'; - for i := 1 to NRows + 1 do - begin - for j := 1 to NCols do - begin - if Freq[i-1, NCols] > 0.0 then - Prop[i-1, j-1] := Freq[i-1, j-1] / Freq[i-1, NCols] - else - Prop[i-1, j-1] := 0.0; - end; - if Freq[i-1, NCols] > 0.0 then - Prop[i-1, NCols] := 1.0 - else - Prop[i-1, NCols] := 0.0; - end; + GetRowProportions(Freq, Prop); lReport.Add('CHI-SQUARE ANALYSIS RESULTS'); lReport.Add(''); + title := 'ROW PROPORTIONS'; MatPrint(Prop, NRows+1, NCols+1, title, RowLabels, ColLabels, NCases, lReport); lReport.Add(DIVIDER_SMALL_AUTO); lReport.Add(''); + GetColProportions(Freq, Prop); title := 'COLUMN PROPORTIONS'; - for j := 1 to Ncols + 1 do - begin - for i := 1 to Nrows do - begin - if Freq[Nrows,j-1] > 0.0 then - Prop[i-1,j-1] := Freq[i-1,j-1] / Freq[Nrows,j-1] - else - Prop[i-1,j-1] := 0.0; - end; - if Freq[Nrows,j-1] > 0.0 then - Prop[NRows,j-1] := 1.0 - else - Prop[NRows,j-1] := 0.0; - end; MatPrint(Prop, NRows+1, NCols+1, title, RowLabels, ColLabels, NCases, lReport); lReport.Add(DIVIDER_SMALL_AUTO); lReport.Add(''); + GetTotalProportions(Freq, Prop); Title := 'PROPORTIONS OF TOTAL N'; - for i := 1 to NRows + 1 do - for j := 1 to NCols + 1 do - Prop[i-1,j-1] := Freq[i-1,j-1] / NCases; - Prop[Nrows, Ncols] := 1.0; MatPrint(Prop, NRows+1, NCols+1, title, RowLabels, ColLabels, NCases, lReport); lReport.Add(DIVIDER_SMALL_AUTO); @@ -438,26 +409,16 @@ begin lReport.Add(''); end; - likelihood := 0.0; - for i := 0 to NRows-1 do - for j := 0 to NCols-1 do - if (Freq[i, j] > 0.0) then - likelihood := Likelihood + (Freq[i, j] * (ln(Expected[i, j] / Freq[i,j ]))); - likelihood := -2.0 * likelihood; - probLikelihood := 1.0 - ChiSquaredProb(likelihood, df); + likelihood := CalcLikelihoodRatio(Freq, Expected); + prob := 1.0 - ChiSquaredProb(likelihood, df); lReport.Add('Likelihood Ratio: %.3f', [likelihood]); - lReport.Add(' Probability > value is %.4f', [probLikelihood]); + lReport.Add(' Probability > value is %.4f', [prob]); lReport.Add(''); - G := 0.0; - for i := 0 to NRows-1 do - for j := 0 to NCols-1 do - if (Expected[i, j] > 0) then - G := G + Freq[i, j] * (ln(Freq[i, j] / Expected[i, j])); - G := 2.0 * G; - probLikelihood := 1.0 - ChiSquaredProb(G, df); + G := CalcGStatistic(Freq, Expected); + prob := 1.0 - ChiSquaredProb(G, df); lReport.Add('G statistic: %.3f ', [G]); - lReport.Add(' Probability > value is %.4f', [G, probLikelihood]); + lReport.Add(' Probability > value is %.4f', [prob]); lReport.Add(''); if ((NRows > 1) and (NCols > 1)) then @@ -466,29 +427,14 @@ begin lReport.Add('phi correlation: %.4f', [phi]); lReport.Add(''); - pearsonr := 0.0; - SumX := 0.0; - SumY := 0.0; - VarX := 0.0; - VarY := 0.0; - for i := 0 to NRows-1 do SumX := SumX + ( (i+1) * Freq[i, NCols] ); - for j := 0 to NCols-1 do SumY := SumY + ( (j+1) * Freq[NRows, j] ); - for i := 0 to NRows-1 do VarX := VarX + ( sqr(i+1) * Freq[i, NCols] ); - for j := 0 to NCols-1 do VarY := VarY + ( sqr(j+1) * Freq[NRows, j] ); - VarX := VarX - sqr(SumX) / NCases; - VarY := VarY - sqr(SumY) / NCases; - for i := 0 to NRows-1 do - for j := 0 to NCols-1 do - pearsonR := pearsonR + ((i+1)*(j+1) * Freq[i, j]); - pearsonR := pearsonR - (SumX * SumY / Ncases); - pearsonR := pearsonR / sqrt(VarX * VarY); + pearsonR := CalcPearsonR(Freq); lReport.Add('Pearson Correlation r: %.4f', [pearsonR]); lReport.Add(''); MantelHaenszel := (NCases-1) * sqr(pearsonR); - MHprob := 1.0 - ChiSquaredProb(MantelHaenszel, 1); + prob := 1.0 - ChiSquaredProb(MantelHaenszel, 1); lReport.Add('Mantel-Haenszel Test of Linear Association: %.3f', [MantelHaenszel]); - lReport.Add(' Probability > value is %.4f', [MantelHaenszel, MHprob]); + lReport.Add(' Probability > value is %.4f', [prob]); lReport.Add(''); CoefCont := sqrt(ChiSquare / (ChiSquare + NCases)); @@ -496,9 +442,9 @@ begin lReport.Add(''); if (Nrows < Ncols) then - CramerV := sqrt(ChiSquare / (NCases * ((NRows-1)))) + CramerV := sqrt(ChiSquare / (NCases * ((NRows-1)))) else - CramerV := sqrt(ChiSquare / (NCases * ((NCols-1)))); + CramerV := sqrt(ChiSquare / (NCases * ((NCols-1)))); lReport.Add('Cramers V is %.3f', [CramerV]); lReport.Add(''); end; @@ -511,57 +457,48 @@ begin // save frequency data file if elected if SaveFChk.Checked then begin - OS3MainFrm.mnuFileCloseClick(self); - OS3MainFrm.FileNameEdit.Text := ''; - for i := 1 to DictionaryFrm.DictGrid.RowCount - 1 do - for j := 0 to 7 do DictionaryFrm.DictGrid.Cells[j,i] := ''; - DictionaryFrm.DictGrid.RowCount := 1; - // DictionaryFrm.FileNameEdit.Text := ''; + OS3MainFrm.mnuFileCloseClick(self); + OS3MainFrm.FileNameEdit.Text := ''; + for i := 1 to DictionaryFrm.DictGrid.RowCount - 1 do + for j := 0 to 7 do DictionaryFrm.DictGrid.Cells[j,i] := ''; + DictionaryFrm.DictGrid.RowCount := 1; - // get labels for new file - ColLabels[0] := 'ROW'; - ColLabels[1] := 'COL'; - ColLabels[2] := 'FREQ'; - // create new variables - Row := 0; - OS3MainFrm.DataGrid.ColCount := 4; - DictionaryFrm.DictGrid.ColCount := 8; - NoVariables := 0; - for i := 1 to 3 do - begin - col := NoVariables + 1; - DictionaryFrm.NewVar(col); - DictionaryFrm.DictGrid.Cells[1,col] := ColLabels[i-1]; - OS3MainFrm.DataGrid.Cells[col,0] := ColLabels[i-1]; - NoVariables := NoVariables + 1; - end; - OS3MainFrm.DataGrid.RowCount := (Nrows * NCols) + 1; - for i := 1 to Nrows do - begin - for j := 1 to Ncols do - begin - Row := Row + 1; - OS3MainFrm.DataGrid.Cells[0,Row] := Format('Case:%d',[Row]); - OS3MainFrm.DataGrid.Cells[1,Row] := IntToStr(i); - OS3MainFrm.DataGrid.Cells[2,Row] := IntToStr(j); - OS3MainFrm.DataGrid.Cells[3,Row] := IntToStr(Freq[i-1,j-1]); - end; - end; - NoCases := Row; - OS3MainFrm.FileNameEdit.Text := 'ChiSqrFreq.LAZ'; - OS3MainFrm.NoCasesEdit.Text := IntToStr(NoCases); - OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables); - // OS3MainFrm.SaveFileBtnClick(self); + // get labels for new file + ColLabels[0] := 'ROW'; + ColLabels[1] := 'COL'; + ColLabels[2] := 'FREQ'; + + // create new variables + Row := 0; + OS3MainFrm.DataGrid.ColCount := 4; + DictionaryFrm.DictGrid.ColCount := 8; + NoVariables := 0; + for i := 1 to 3 do + begin + col := NoVariables + 1; + DictionaryFrm.NewVar(col); + DictionaryFrm.DictGrid.Cells[1,col] := ColLabels[i-1]; + OS3MainFrm.DataGrid.Cells[col,0] := ColLabels[i-1]; + NoVariables := NoVariables + 1; + end; + OS3MainFrm.DataGrid.RowCount := (Nrows * NCols) + 1; + + for i := 1 to Nrows do + begin + for j := 1 to Ncols do + begin + Row := Row + 1; + OS3MainFrm.DataGrid.Cells[0,Row] := Format('Case:%d',[Row]); + OS3MainFrm.DataGrid.Cells[1,Row] := IntToStr(i); + OS3MainFrm.DataGrid.Cells[2,Row] := IntToStr(j); + OS3MainFrm.DataGrid.Cells[3,Row] := IntToStr(Freq[i-1,j-1]); + end; + end; + NoCases := Row; + OS3MainFrm.FileNameEdit.Text := 'ChiSqrFreq.LAZ'; + OS3MainFrm.NoCasesEdit.Text := IntToStr(NoCases); + OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables); end; - - //clean up - ColLabels := nil; - RowLabels := nil; - CellChi := nil; - Expected := nil; - Prop := nil; - Freq := nil; - ColNoSelected := nil; end; @@ -590,6 +527,207 @@ begin end; +procedure TChiSqrFrm.GetExpectedAndCellChiSqr(const AFrequencies: IntDyneMat; + out AExpected, ACellChiSqr: DblDyneMat; + out AChiSqr: Double); +var + numRows, numCols, numCases: Integer; + i, j: Integer; +begin + MatSize(AFrequencies, numRows, numCols); // contains the totals row/col + numCases := AFrequencies[numRows-1, numCols-1]; + + SetLength(AExpected, numRows-1, numCols-1); // -1: we don't need the totals here + SetLength(ACellChiSqr, numRows-1, numCols-1); + + AChiSqr := 0; + for i := 0 to numRows-2 do // -2 instead of -1 to skip the totals row + begin + for j := 0 to numCols-2 do // -2 instead of -1 to skip the totals column + begin + AExpected[i, j] := AFrequencies[numRows-1, j] * AFrequencies[i, numCols-1] / numCases; + if AExpected[i, j] > 0 then + ACellChiSqr[i, j] := sqr(AFrequencies[i, j] - AExpected[i, j]) / AExpected[i, j] + else begin + ErrorMsg('Zero expected value found.'); + ACellChiSqr[i, j] := 0; + end; + AChiSqr := AChiSqr + ACellChiSqr[i, j]; + end; + end; +end; + + +procedure TChiSqrFrm.GetFrequencies(const AColNoSelected: IntDyneVec; + ARowIndex, AColIndex, ADepIndex: Integer; out ANumRows, ANumCols, ANumCases: Integer; + out AFrequencies: IntDyneMat); +var + i, j: Integer; + row, col: Integer; + minRow, maxRow, minCol, maxCol: Integer; + FObs: Integer; + PObs: Double; +begin + // Get min and max of row and col numbers + minRow := MaxInt; + maxRow := -MinRow; + minCol := MaxInt; + maxCol := -MinCol; + for i := 1 to NoCases do + begin + if not GoodRecord(OS3MainFrm.DataGrid, i, AColNoSelected) then continue; + row := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ARowIndex, i]))); + col := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[AColIndex, i]))); + if row > maxRow then maxRow := row; + if row < minRow then minRow := row; + if col > maxCol then maxCol := col; + if col < minCol then minCol := col; + end; + ANumRows := maxRow - minRow + 1; + ANumCols := maxCol - minCol + 1; + + AFrequencies := nil; + SetLength(AFrequencies, ANumRows+1, ANumCols+1); // +1 for row and column totals + + ANumCases := 0; + case InputGrp.ItemIndex of + 0 : begin // count number of cases in each row and column combination + for i := 1 to NoCases do + begin + if not GoodRecord(OS3MainFrm.DataGrid, i, AColNoSelected) then continue; + row := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ARowIndex, i]))); + col := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[AColIndex, i]))); + row := row - minRow; + col := col - minCol; + AFrequencies[row, col] := AFrequencies[row, col] + 1; + ANumCases := ANumCases + 1; + end; + end; + 1 : begin // read frequencies data from grid + for i := 1 to NoCases do + begin + if not GoodRecord(OS3MainFrm.DataGrid, i, AColNoSelected) then continue; + row := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ARowIndex, i]))); + col := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ARowIndex, i]))); + row := row - minRow; + col := col - minCol; + FObs := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ADepIndex, i]))); + AFrequencies[row, col] := AFrequencies[row, col] + FObs; + ANumCases := ANumCases + FObs; + end; + end; + 2 : begin // get no. of cases and proportions for each cell + ANumCases := StrToInt(NCasesEdit.Text); + for i := 1 to NoCases do + begin + if not GoodRecord(OS3MainFrm.Datagrid, i, AColNoSelected) then continue; + row := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ARowIndex, i]))); + col := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[AColIndex, i]))); + row := row - minRow + 1; + col := col - minCol + 1; + PObs := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ADepIndex, i])); + AFrequencies[row, col] := AFrequencies[row, col] + round(PObs * ANumCases); + end; + end; + end; + + // Get row totals + for i := 0 to ANumRows-1 do + for j := 0 to ANumCols-1 do + AFrequencies[i, ANumCols] := AFrequencies[i, ANumCols] + AFrequencies[i, j]; + + // Get col totals + for j := 0 to ANumCols-1 do + for i := 0 to ANumRows-1 do + AFrequencies[ANumRows, j] := AFrequencies[ANumRows, j] + AFrequencies[i, j]; + + // Grand total + AFrequencies[ANumRows, ANumCols] := ANumCases; +end; + + +// 2 x 2 corrected chi-square +procedure TChiSqrFrm.GetYatesCorrection(const AFrequencies: IntDyneMat; + out AdjChiSqr: Double); +var + numRows, numCols, numCases: Integer; +begin + MatSize(AFrequencies, numRows, numCols); + numCases := AFrequencies[numRows-1, numCols-1]; + + AdjChiSqr := abs((AFrequencies[0,0] * AFrequencies[1,1]) - (AFrequencies[0,1] * AFrequencies[1,0])); + AdjChiSqr := sqr(AdjChiSqr - numCases / 2.0) * numCases; // numerator + AdjChiSqr := AdjChiSqr / (AFrequencies[0,2] * AFrequencies[1,2] * AFrequencies[2,0] * AFrequencies[2,1]); +end; + + +procedure TChiSqrFrm.GetColProportions(const AFrequencies: IntDyneMat; + out AProportions: DblDyneMat); +var + numRows, numCols: Integer; + i, j: Integer; +begin + MatSize(AFrequencies, numRows, numCols); // totals in last row/col + SetLength(AProportions, numRows, numCols); + for j := 0 to numCols-1 do + begin + for i := 0 to numRows-1 do // Do not process the totals row + begin + if AFrequencies[numRows-1, j] > 0.0 then + AProportions[i, j] := AFrequencies[i, j] / AFrequencies[numRows-1, j] + else + AProportions[i, j] := 0.0; + end; + if AFrequencies[numRows-1, j] > 0.0 then + AProportions[numRows-1,j] := 1.0 + else + AProportions[numRows-1,j] := 0.0; + end; +end; + + +procedure TChiSqrFrm.GetRowProportions(const AFrequencies: IntDyneMat; + out AProportions: DblDyneMat); +var + numRows, numCols: Integer; + i, j: Integer; +begin + MatSize(AFrequencies, numRows, numCols); // totals in last row/col + SetLength(AProportions, numRows, numCols); + for i := 0 to numRows-1 do + begin + for j := 0 to numCols-1 do // do not touch the totals column here + begin + if AFrequencies[i, numCols-1] > 0.0 then + AProportions[i, j] := AFrequencies[i, j] / AFrequencies[i, numCols-1] + else + AProportions[i-1, j-1] := 0.0; + end; + if AFrequencies[i, numCols-1] > 0.0 then + AProportions[i, numCols-1] := 1.0 + else + AProportions[i, numCols-1] := 0.0; + end; +end; + + +procedure TChiSqrFrm.GetTotalProportions(const AFrequencies: IntDyneMat; + out AProportions: DblDyneMat); +var + numRows, numCols, numCases: Integer; + i, j: Integer; +begin + MatSize(AFrequencies, numRows, numCols); // totals in last row/col + numCases := AFrequencies[numRows-1, numCols-1]; + + SetLength(AProportions, numRows, numCols); + for i := 0 to numRows-1 do + for j := 0 to numCols-1 do + AProportions[i, j] := AFrequencies[i, j] / numCases; + AProportions[numRows-1, numCols-1] := 1.0; +end; + + procedure TChiSqrFrm.InputGrpClick(Sender: TObject); begin case InputGrp.ItemIndex of diff --git a/applications/lazstats/source/units/matrixunit.pas b/applications/lazstats/source/units/matrixunit.pas index 78ac80745..e9bbe4992 100644 --- a/applications/lazstats/source/units/matrixunit.pas +++ b/applications/lazstats/source/units/matrixunit.pas @@ -14,6 +14,8 @@ type TDblMatrix = DblDyneMat; TDblVector = DblDyneVec; + TIntMatrix = IntDyneMat; + EMatrix = class(Exception); // Vectors @@ -79,6 +81,7 @@ function MatNumRows(A: TDblMatrix): Integer; function MatRowMeans(A: TDblMatrix): TDblVector; function MatRowVector(A: TDblMatrix; ARowIndex: Integer): TDblVector; procedure MatSize(A: TDblMatrix; out n, m: Integer); +procedure MatSize(A: TIntMatrix; out n, m: Integer); function MatTransposed(A: TDblMatrix): TDblMatrix; function SubMatrix(A: TDblMatrix; i1,j1, i2,j2: Integer): TDblMatrix; @@ -699,6 +702,13 @@ begin end; +procedure MatSize(A: TIntMatrix; out n, m: Integer); +begin + n := Length(A); + m := Length(A[0]); +end; + + function MatTransposed(A: TDblMatrix): TDblMatrix; var n, m, i, j: Integer;