diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/NOTICE.txt b/NOTICE.txt new file mode 100644 --- /dev/null +++ b/NOTICE.txt @@ -0,0 +1,48 @@ +Copyright 2020 Ideas project team + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. + +PACKAGE + + Interactive domain reasoner for statistics, based on the feedback services + of the Ideas framework. Provides feedback and hints on step-by-step + construction of hypothesis tests by students. Test types supported include + z-tests, t-tests for one group, dependent groups and independent groups, + correlation, ANOVA and Chi-square tests. The domain reasoner is used by the + learning environment Numworx (www.numworx.nl) and is available for + implementation in other learning environments as well. + +PROJECT + + Ideas (Interactive Domain-specific Exercise Assistants) is a joint research + project between the Open University of the Netherlands and Utrecht + University. The project's goal is to use software and compiler technology to + build state-of-the-art components of intelligent tutoring systems (ITS) and + learning environments. The 'ideas' software package provides a generic + framework for constructing the expert knowledge module (also known as a + domain reasoner) for an ITS or learning environment. More information about + the project can be found on the project's homepage: + + http://ideas.cs.uu.nl/ + +AUTHORS + + The authors and copyright holders of the 'ideas-statistics' software package are: + - Bastiaan Heeren (bastiaan.heeren@ou.nl, maintainer) + - Sietske Tacoma + - Noeri Huisman + +CREDITS + + Paul Drijvers, Johan Jeuring, Martijn Fleuren, Wim van Velthoven, Peter Boon, + Jeltje Wassenberg-Severijnen, Corine Geurts, Alex Gerdes diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,4 @@ +#! /usr/bin/env runhaskell + +> import Distribution.Simple +> main = defaultMain diff --git a/ideas-statistics.cabal b/ideas-statistics.cabal new file mode 100644 --- /dev/null +++ b/ideas-statistics.cabal @@ -0,0 +1,62 @@ +name: ideas-statistics +version: 1.0 +synopsis: Interactive domain reasoner for statistics +homepage: http://ideas.cs.uu.nl/ +description: + + Interactive domain reasoner for statistics, based on the feedback services + of the Ideas framework. Provides feedback and hints on step-by-step construction + of hypothesis tests by students. Test types supported include z-tests, t-tests + for one group, dependent groups and independent groups, correlation, ANOVA and + Chi-square tests. The domain reasoner is used by the learning environment Numworx + (www.numworx.nl) and is available for implementation in other learning + environments as well. + +category: Education +copyright: (c) 2020 +license: Apache-2.0 +license-file: LICENSE.txt +author: Bastiaan Heeren, Sietske Tacoma, Noeri Huisman +maintainer: bastiaan.heeren@ou.nl +stability: provisional +extra-source-files: NOTICE.txt +build-type: Simple +cabal-version: >= 1.8.0.2 +tested-with: GHC == 8.6.5 + +source-repository head + type: svn + location: https://ideastest.science.uu.nl/svn/ideas/Tutors/statistics + +-------------------------------------------------------------------------------- + +executable statistics.cgi + main-is: Main.hs + hs-source-dirs: src + ghc-options: -Wall + + build-depends: + base >= 4.8 && < 5, + containers, + ideas >= 1.8, + ideas-math-types >= 1.1 + + other-modules: + Domain.Hypothesis.BuggyRules, + Domain.Hypothesis.Common, + Domain.Hypothesis.Constraints, + Domain.Hypothesis.Examples, + Domain.Hypothesis.Exercises, + Domain.Hypothesis.Rules, + Domain.Hypothesis.Strategies, + Domain.Hypothesis.Tables, + Domain.Statistics.Component, + Domain.Statistics.ComponentSet, + Domain.Statistics.Data, + Domain.Statistics.Parser, + Domain.Statistics.Rules, + Domain.Statistics.Symbols, + Domain.Statistics.Views, + Service.ConstraintServices + +-------------------------------------------------------------------------------- diff --git a/src/Domain/Hypothesis/BuggyRules.hs b/src/Domain/Hypothesis/BuggyRules.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Hypothesis/BuggyRules.hs @@ -0,0 +1,406 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Hypothesis.BuggyRules + ( buggyRules + , buggyTValueSided, buggyRValueSided, buggyZValueSided + , buggyTValueTestValue, buggyRValueTestValue, buggyZValueTestValue + , buggyChiValueSided, buggyChiValueTestValue + ) where + +import Control.Monad +import Data.List +import Domain.Hypothesis.Common +import Domain.Math.Data.Relation +import Domain.Math.Expr hiding ((.*.), (./.)) +import Domain.Math.Numeric.Views +import Domain.Hypothesis.Rules +import Domain.Statistics.ComponentSet +import Domain.Statistics.Views +import Ideas.Common.Library + +---------------------------------------------------------- +-- Buggy rules + +-- Transform the steps described below into a rule, the pattern is factored out. +-- The list can be freely appended with new buggy rules +buggyRules:: [Rule ComponentSet] +buggyRules = + let _pat (f, name, descr) = (describe descr . buggyRule name) f + _name_common x = "buggy.component." ++ x + in map _pat + [ ( wrongHARule + , _name_common "ha" + , "Choosing wrong alternative hypothesis" + ) + , ( wrongAlphaRule + , _name_common "alpha" + , "Buggy rule for adding the wrong alpha" + ) + , ( wrongDf + , _name_common "df" + , "Buggy rule for adding the wrong degrees of freedom" + ) + ] ++ [buggyHypothesesSampleMean, buggyHASampleMean, buggyH0SampleMean] + ++ buggyRulesTest + ++ [ buggyTestZValueRule, buggyTestTValueRule + , buggyTValuePositive + , buggyZValueAlpha, buggyZValueSided, buggyZValueTestValue + , buggyTValueAlpha, buggyTValueSided, buggyTValueDf, buggyTValueTestValue + , buggyRValueAlpha, buggyRValueSided, buggyRValueDf, buggyRValueTestValue + , buggyFValueAlpha, buggyFValueDf, buggyFValueTestValue + , buggyChiValueAlpha, buggyChiValueSided, buggyChiValueDf, buggyChiValueTestValue + , buggyTestValueCritical + ] + +alphaValues :: [Double] +alphaValues = [0.01, 0.05, 0.10] + +alphaChiValues :: [Double] +alphaChiValues = [0.10, 0.05, 0.025, 0.01, 0.005] + +sidedValues :: [Sided] +sidedValues = [TwoSided, LeftSided, RightSided] + +incorrectDf :: Double -> [Double] +incorrectDf correctDf = [correctDf - 1, correctDf + 1, correctDf + 2] + +confusingTestValue :: ComponentSet -> [Double] +confusingTestValue cs = do + r <- inferTestValue cs + matchM doubleView (rightHandSide r) + +buggyTestValueCritical :: Rule ComponentSet +buggyTestValueCritical = + describe "Buggy rule for confusing the test value with critical value" $ + buggyRule "buggy.component.test-value.critical" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` TestValue) + test <- inferTestChoice cs + var <- inferVar cs + alpha <- getExpr SignificanceLevel cs >>= match doubleView + crit <- + case test of + Anova -> do + (dfBetween, dfWithin) <- inferDfBetweenWithin cs + computeCriticalF dfBetween dfWithin alpha + ZTest -> do + sided <- inferSidedness cs + return $ computeCriticalZ sided alpha + RPearson -> do + sided <- inferSidedness cs + df <- inferDf cs >>= match doubleView + return $ computeCriticalR sided alpha df + _ | isTTest test -> do + sided <- inferSidedness cs + df <- inferDf cs >>= match doubleView + return $ computeCriticalT sided alpha df + _ -> Nothing + return $ append TestValue (CRelation $ var .==. toExpr crit) cs + +------------------------------------------------------ + +buggyZValueAlpha :: Rule ComponentSet +buggyZValueAlpha = + describe "Buggy rule for looking up a z-value with wrong alpha" . + buggyRule "buggy.component.critical.z-value.alpha" $ + inferCriticalZWith $ \sided alpha -> + [ computeCriticalZ sided a | a <- alphaValues, a /= alpha ] + +buggyZValueSided :: Rule ComponentSet +buggyZValueSided = + describe "Buggy rule for looking up a z-value with wrong sidedness" . + buggyRule "buggy.component.critical.z-value.sided" $ + inferCriticalZWith $ \sided alpha -> + [ computeCriticalZ s alpha | s <- sidedValues, s /= sided ] + +buggyZValueTestValue :: Rule ComponentSet +buggyZValueTestValue = + describe "Buggy rule for confusing z-value with test-value" . + buggyRule "buggy.component.critical.z-value.test-value" $ \cs -> + inferCriticalZWith (\_ _ -> confusingTestValue cs) cs + +------------------------------------------------------ + +buggyTValueAlpha :: Rule ComponentSet +buggyTValueAlpha = + describe "Buggy rule for looking up a t-value with wrong alpha" . + buggyRule "buggy.component.critical.t-value.alpha" $ + inferCriticalTWith $ \sided alpha df -> + [ computeCriticalT sided a df | a <- alphaValues, a /= alpha ] + +buggyTValueSided :: Rule ComponentSet +buggyTValueSided = + describe "Buggy rule for looking up a t-value with wrong sidedness" . + buggyRule "buggy.component.critical.t-value.sided" $ + inferCriticalTWith $ \sided alpha df -> + [ computeCriticalT s alpha df | s <- sidedValues, s /= sided ] + +buggyTValueDf :: Rule ComponentSet +buggyTValueDf = + describe "Buggy rule for looking up a t-value with wrong df" . + buggyRule "buggy.component.critical.t-value.df" $ + inferCriticalTWith $ \sided alpha df -> + [ computeCriticalT sided alpha dfWrong | dfWrong <- incorrectDf df ] + +buggyTValueTestValue :: Rule ComponentSet +buggyTValueTestValue = + describe "Buggy rule for confusing t-value with test-value" . + buggyRule "buggy.component.critical.t-value.test-value" $ \cs -> + inferCriticalTWith (\_ _ _ -> confusingTestValue cs) cs + +------------------------------------------------------ + +buggyRValueAlpha :: Rule ComponentSet +buggyRValueAlpha = + describe "Buggy rule for looking up a r-value with wrong alpha" . + buggyRule "buggy.component.critical.r-value.alpha" $ + inferCriticalRWith $ \sided alpha df -> + [ computeCriticalR sided a df | a <- alphaValues, a /= alpha ] + +buggyRValueSided :: Rule ComponentSet +buggyRValueSided = + describe "Buggy rule for looking up a r-value with wrong sidedness" . + buggyRule "buggy.component.critical.r-value.sided" $ + inferCriticalRWith $ \sided alpha df -> + [ computeCriticalR s alpha df | s <- sidedValues, s /= sided ] + +buggyRValueDf :: Rule ComponentSet +buggyRValueDf = + describe "Buggy rule for looking up a r-value with wrong df" . + buggyRule "buggy.component.critical.r-value.df" $ + inferCriticalRWith $ \sided alpha df -> + [ computeCriticalR sided alpha dfWrong | dfWrong <- incorrectDf df ] + +buggyRValueTestValue :: Rule ComponentSet +buggyRValueTestValue = + describe "Buggy rule for confusing r-value with test-value" . + buggyRule "buggy.component.critical.r-value.test-value" $ \cs -> + inferCriticalRWith (\_ _ _ -> confusingTestValue cs) cs + +------------------------------------------------------ + +buggyFValueAlpha :: Rule ComponentSet +buggyFValueAlpha = + describe "Buggy rule for looking up f-value with wrong alpha" . + buggyRule "buggy.component.critical.f-value.alpha" $ + inferCriticalFWith $ \dfBetween dfWithin alpha -> + [ result + | a <- [0.01, 0.05], a /= alpha + , result <- computeCriticalF dfBetween dfWithin a + ] + +buggyFValueDf :: Rule ComponentSet +buggyFValueDf = + describe "Buggy rule for looking up f-value with wrong df" . + buggyRule "buggy.component.critical.f-value.df" $ + inferCriticalFWith $ \dfBetween dfWithin alpha -> + [ result + | dfB <- [dfBetween, 2] + , dfW <- [dfWithin, 40] + , dfBetween /= dfB || dfWithin /= dfW + , result <- computeCriticalF dfB dfW alpha + ] + +buggyFValueTestValue :: Rule ComponentSet +buggyFValueTestValue = + describe "Buggy rule for confusing f-value with test-value" . + buggyRule "buggy.component.critical.f-value.test-value" $ \cs -> + inferCriticalFWith (\_ _ _ -> confusingTestValue cs) cs + +------------------------------------------------------ + +buggyChiValueSided :: Rule ComponentSet +buggyChiValueSided = + describe "Buggy rule for looking up a chi-value with wrong sidedness" . + buggyRule "buggy.component.critical.chi-value.sided" $ + inferCriticalChiWith $ \sided alpha df -> + [ result + | s <- sidedValues + , s /= sided + , result <- computeCriticalChi s alpha df + ] + +buggyChiValueAlpha :: Rule ComponentSet +buggyChiValueAlpha = + describe "Buggy rule for looking up chi-value with wrong alpha" . + buggyRule "buggy.component.critical.chi-value.alpha" $ + inferCriticalChiWith $ \sided alpha df -> + [ result + | a <- alphaChiValues, a /= alpha + , result <- computeCriticalChi sided a df + ] + +buggyChiValueDf :: Rule ComponentSet +buggyChiValueDf = + describe "Buggy rule for looking up chi-value with wrong df" . + buggyRule "buggy.component.critical.chi-value.df" $ + inferCriticalChiWith $ \sided alpha df -> + [ result + | dfWrong <- incorrectDf df + , result <- computeCriticalChi sided alpha dfWrong + ] + +buggyChiValueTestValue :: Rule ComponentSet +buggyChiValueTestValue = + describe "Buggy rule for confusing chi-value with test-value" . + buggyRule "buggy.component.critical.chi-value.test-value" $ \cs -> + inferCriticalChiWith (\_ _ _ -> confusingTestValue cs) cs + +------------------------------------------------------ + +buggyHypothesesSampleMean :: Rule ComponentSet +buggyHypothesesSampleMean = + describe "Buggy rule for adding one or two hypotheses using sample mean instead of population mean" $ + buggyRule "buggy.component.hypotheses-samplemean" f + where + f :: ComponentSet -> [ComponentSet] + f cs = do + guard (derived cs `doesNotContain` NullHypothesis) + guard (derived cs `doesNotContain` AlternativeHypothesis) + ha <- getRelation AlternativeHypothesis cs + let h0 = h0FromHA ha + sm <- getExpr SampleMean cs + rtp <- nub [EqualTo, relationType h0] + let h0' = makeType rtp (leftHandSide h0) (rightHandSide h0) + buggyH0 = makeType rtp (leftHandSide h0) sm + buggyHA = makeType (relationType ha) (leftHandSide ha) sm + results = + [ append AlternativeHypothesis (CRelation buggyHA) $ + append NullHypothesis (CRelation buggyH0) cs + , append AlternativeHypothesis (CRelation ha) $ + append NullHypothesis (CRelation buggyH0) cs + , append AlternativeHypothesis (CRelation buggyHA) $ + append NullHypothesis (CRelation h0') cs + ] + results + +buggyH0SampleMean :: Rule ComponentSet +buggyH0SampleMean = + describe "Buggy rule for adding h0 using sample mean instead of population mean" $ + buggyRule "buggy.component.h0-samplemean" f + where + f :: ComponentSet -> [ComponentSet] + f cs = do + guard (derived cs `doesNotContain` NullHypothesis) + ha <- getRelation AlternativeHypothesis cs + let h0 = h0FromHA ha + sm <- getExpr SampleMean cs + rtp <- nub [EqualTo, relationType h0] + let buggyH0 = makeType rtp (leftHandSide h0) sm + return $ append NullHypothesis (CRelation buggyH0) cs + +buggyHASampleMean :: Rule ComponentSet +buggyHASampleMean = + describe "Buggy rule for adding ha using sample mean instead of population mean" $ + buggyRule "buggy.component.ha-samplemean" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` AlternativeHypothesis) + ha <- getRelation AlternativeHypothesis cs + sm <- getExpr SampleMean cs + let buggyHA = makeType (relationType ha) (leftHandSide ha) sm + return $ append AlternativeHypothesis (CRelation buggyHA) cs + +buggyTestZValueRule :: Rule ComponentSet +buggyTestZValueRule = + describe "Standard error instead of standard deviation" $ + buggyRule "buggy.component.test-z-value" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` TestValue) + let expr = (Var "M" - Var "mu") / Var "sigma" + val <- match doubleView (getSubstitution cs |-> expr) + var <- inferVar cs + return $ append TestValue (CRelation (var .==. toExpr val)) cs + +-- see buggyTestZValueRule: only difference is 's' instead of 'sigma' +buggyTestTValueRule :: Rule ComponentSet +buggyTestTValueRule = + describe "Standard error instead of standard deviation" $ + buggyRule "buggy.component.test-t-value" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` TestValue) + let expr = (Var "M" - Var "mu") / Var "s" + val <- match doubleView (getSubstitution cs |-> expr) + var <- inferVar cs + return $ append TestValue (CRelation (var .==. toExpr val)) cs + +buggyTValuePositive :: Rule ComponentSet +buggyTValuePositive = + buggyRule "buggy.component.critical.t-value-positive" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + cs' <- apply lookupTValueRule cs + rel <- getRelation Critical cs' + val <- match doubleView (rightHandSide rel) + guard (val < 0) + return $ append Critical (CRelation $ Var "tcrit" .==. toExpr (abs val)) cs + +buggyRulesTest :: [Rule ComponentSet] +buggyRulesTest = + [ buggyRule (makeId wrong correct) (wrongTest wrong correct) + | wrong <- types + , correct <- types + , wrong /= correct + ] + where + types = [TTestOne, TTestTwo, TTestPaired, ZTest] + makeId t1 t2 = "buggy.component." ++ show t1 ++ ".should-be-" ++ show t2 + +wrongTest :: TestType -> TestType -> ComponentSet -> Maybe ComponentSet +wrongTest wrong correct cs = do + guard (derived cs `doesNotContain` TestChoice) + let tests = validTests cs + guard (correct `elem` tests) + return $ append TestChoice (CChoice $ TestType wrong) cs + +-- Wrong sidedness in alternative hypothesis (due to sample mean) +wrongHARule :: ComponentSet -> Maybe ComponentSet +wrongHARule cs = + let relFromMean :: Double -> Double -> RelationType + relFromMean mean mu + | mean < mu = LessThanOrEqualTo + | otherwise = GreaterThanOrEqualTo + in do + guard (cs `doesNotContain` AlternativeHypothesis) + ha <- getRelation AlternativeHypothesis cs -- used to be Claim + let h0 = h0FromHA ha + h0Rel = relationType h0 + guard (h0Rel == EqualTo) + let lhs = leftHandSide ha + rhs = rightHandSide ha + rhs' <- match doubleView rhs + mean <- match doubleView =<< getExpr SampleMean cs + let rel = relFromMean mean rhs' + return . append AlternativeHypothesis + (CRelation (makeType rel lhs rhs)) $ cs + +wrongAlphaRule :: ComponentSet -> Maybe ComponentSet +wrongAlphaRule cs = do + guard (cs `doesNotContain` SignificanceLevel) + return $ append SignificanceLevel (CExpr $ Number 0.1) cs + +-- | Use the wrong df: t-tests are swapped so n instead of n-1 or n_1 + n_2 - 1 +-- instead of n_1 + n_2 - 2 +wrongDf :: ComponentSet -> Maybe ComponentSet +wrongDf cs = + let dfWrong :: TestType -> Maybe Expr + dfWrong TTestOne = Just (Var "n") + dfWrong TTestPaired = Just (Var "n") + dfWrong TTestTwo = Just (Var "n1" + Var "n2" - 1) + dfWrong _ = Nothing + in do + guard (cs `doesNotContain` Df) + test <- inferTestChoice cs + df <- dfWrong test + return $ append Df (CExpr df) cs diff --git a/src/Domain/Hypothesis/Common.hs b/src/Domain/Hypothesis/Common.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Hypothesis/Common.hs @@ -0,0 +1,267 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Hypothesis.Common + ( TypeOfTest(..), + dependent, independent, chiSquared, + pickAlpha, sidedFromHA, + h0FromHA, h0FromHAEqualSign, + validTests, testFormulaFromTest, + degreesOfFreedomFromTest, + computeCritical, computeCriticalT, computeCriticalZ + , computeCriticalR, computeCriticalF, computeCriticalChi, + computePValue, computePValueT, computePValueZ, + chooseTypeOfTest, isTTest + ) where + +import Data.Maybe +import Ideas.Common.Library (Term(TCon)) +import Domain.Statistics.Symbols +import Domain.Hypothesis.Tables +import Domain.Math.Expr hiding ((^)) +import Domain.Math.Data.Relation +import Domain.Statistics.Component +import Domain.Statistics.ComponentSet +import Ideas.Common.View +import Domain.Math.Numeric.Views + +data TypeOfTest = TestMean + | CompareMeans + | CompareMeansPaired + | TestProportion + | CompareProportions + | TestCorrelation + deriving (Eq, Show) + +dependent, independent, chiSquared :: Expr +dependent = toExpr (TCon dependentSymbol []) +independent = toExpr (TCon independentSymbol []) +chiSquared = Var "chisq" + +-- Pick a value for alpha +pickAlpha :: Double +pickAlpha = 0.05 + +-- Determine the sidedness based on the alternative hypothesis. +sidedFromHA :: Relation Expr -> Sided +sidedFromHA ha = fromRelation (relationType ha) + where + fromRelation :: RelationType -> Sided + fromRelation LessThan = LeftSided + fromRelation LessThanOrEqualTo = LeftSided + fromRelation GreaterThan = RightSided + fromRelation GreaterThanOrEqualTo = RightSided + fromRelation EqualTo = TwoSided + fromRelation NotEqualTo = TwoSided + fromRelation rt = error $ "sidedFromHA: " ++ show rt + +h0FromHA :: Relation Expr -> Relation Expr +h0FromHA ha = let h0Rel = inverseRelType (relationType ha) + in makeType h0Rel (leftHandSide ha) (rightHandSide ha) + +h0FromHAEqualSign :: Relation Expr -> Relation Expr +h0FromHAEqualSign ha = makeType EqualTo (leftHandSide ha) (rightHandSide ha) + +inverseRelType :: RelationType -> RelationType +inverseRelType relType = fromMaybe relType (lookup relType table) + where + table = pairs ++ map (\(a,b) -> (b,a)) pairs + pairs = [(LessThan, GreaterThanOrEqualTo), (LessThanOrEqualTo, GreaterThan), (EqualTo, NotEqualTo)] + +-- Returns the valid tests for this state +-- Corner case for two samples (n1,n2 instead of n) +validTests :: ComponentSet -> [TestType] +validTests cs = + case getTestType TestChoice (initials cs) of + Just x -> [x] + _ -> select (chooseTypeOfTest cs) + where + n = fromMaybe 0 $ do + expr <- getExpr SampleSize cs + match naturalView expr + + psdKnown = cs `contains` PopulationSdev + + select TestMean + | psdKnown = [ZTest] + | n > 100 = [ZTest] -- Sietske: ignore TTestOne for now + | otherwise = [TTestOne] + select CompareMeans = [TTestTwo] + select CompareMeansPaired = [TTestPaired] + select TestProportion = [ZTest] + select CompareProportions = [ZTest] + select TestCorrelation = [RPearson, TTestOne] + +-- Returns the test statistic formula from the chose test. +-- Note: this only works for the case of testing the mean, other cases should +-- give the formula as part of the exercise +testFormulaFromTest :: Monad m => TestType -> TypeOfTest -> m Expr +testFormulaFromTest testType typeOfTest = + case (testType, typeOfTest) of + (TTestOne, TestMean) -> + return $ (Var "M" - Var "mu") / (Var "s" / sqrt (Var "n")) + (ZTest, TestMean) -> + return $ (Var "M" - Var "mu") / (Var "sigma" / sqrt (Var "n")) + (TTestTwo, CompareMeans) -> + return $ (mean1 - mean2) / sqrt (toExpr PooledVariance * (1 / Var "n1" + 1 / Var "n2")) + (TTestPaired, TestMean) -> + return $ (Var "M" - Var "mu") / (Var "s" / sqrt (Var "n")) + (TTestPaired, CompareMeansPaired) -> + return $ (mean1 - mean2) / (toExpr SampleSdev * sqrt (1 / Var "n")) -- to do: ask Sietske + (_, TestProportion) -> + return $ (Var "p" - Var "p0") / sqrt (Var "p0" * (1.0 - Var "p0") / Var "n") + (_, CompareProportions) -> + return $ (Var "p1" - Var "p2" - Var "d0") / sqrt (Var "p0" * (1.0 - Var "p0") / (Var "n" / 2)) + (TTestOne, TestCorrelation) -> + return $ (Var "r" * sqrt (Var "n" - 2)) / sqrt (1 - Var "r"**2) + (RPearson, TestCorrelation) -> + return $ Var "r" + _ -> + fail $ "teststatisticFromTest " ++ show (testType, typeOfTest) + where + mean1 = toExpr (One SampleMean) + mean2 = toExpr (Two SampleMean) + +-- Returns the formula for the degrees of degreesOfFreedom +degreesOfFreedomFromTest :: Monad m => TestType -> TypeOfTest -> m Expr +degreesOfFreedomFromTest TTestOne TestCorrelation = return $ Var "n" - 2 +degreesOfFreedomFromTest TTestOne _ = return $ Var "n" - 1 +degreesOfFreedomFromTest TTestPaired _ = return $ Var "n" - 1 +degreesOfFreedomFromTest TTestTwo _ = return $ Var "n1" + Var "n2" - 2 +degreesOfFreedomFromTest RPearson _ = return $ Var "n" - 2 +degreesOfFreedomFromTest _ _ = fail "degrees of freedom test failed" + +-- Returns the critical value from the given test and alpha +computeCritical :: Monad m => TestType -> Sided -> Double -> Maybe Double -> m Double +computeCritical test sided alpha mdf + | isTTest test = case mdf of + Just df -> return $ computeCriticalT sided alpha df + Nothing -> fail "df missing" + | test == RPearson = case mdf of + Just df -> return $ computeCriticalR sided alpha df + Nothing -> fail "df missing" + | test == ZTest = return $ computeCriticalZ sided alpha + | otherwise = fail "unknown test" + +isTTest :: TestType -> Bool +isTTest TTestOne = True +isTTest TTestTwo = True +isTTest TTestPaired = True +isTTest _ = False + +computeCriticalR :: Sided -> Double -> Double -> Double +computeCriticalR sided alpha df = sqrt (t ^ (2 :: Int) / (t ^ (2 :: Int) + df)) + where + t = computeCriticalT sided alpha df + +computeCriticalF :: Monad m => Double -> Double -> Double -> m Double +computeCriticalF dfBetween dfWithin alpha = + maybe (fail "unknown critical-f value") return (fTable dfBetween dfWithin alpha) -- TODO Sietske + + +computeCriticalChi :: Monad m => Sided -> Double -> Double -> m Double +computeCriticalChi TwoSided alpha df = chivalue' df (alpha / 2) +computeCriticalChi LeftSided alpha df = negate <$> chivalue' df alpha +computeCriticalChi RightSided alpha df = chivalue' df alpha + +chivalue' :: Monad m => Double -> Double -> m Double +chivalue' df alpha = + maybe (fail "unknown critical-chi value") return $ chiTable alpha (round df) + +computeCriticalT :: Sided -> Double -> Double -> Double +computeCriticalT TwoSided alpha df = tvalue' df (alpha / 2) +computeCriticalT LeftSided alpha df = - tvalue' df alpha +computeCriticalT RightSided alpha df = tvalue' df alpha + +tvalue' :: Double -> Double -> Double +tvalue' df alpha | isJust tableLookup = fromJust tableLookup + | otherwise = fromInteger (round $ findValue (tvalue df) 0.00005 (0.5 - alpha) * 1000) / 1000.0 + where + tableLookup = tTable alpha (round df) + +computeCriticalZ :: Sided -> Double -> Double +computeCriticalZ TwoSided alpha = zvalue' (alpha / 2) +computeCriticalZ LeftSided alpha = - zvalue' alpha +computeCriticalZ RightSided alpha = zvalue' alpha + +zvalue' :: Double -> Double +zvalue' alpha | isJust tableLookup = fromJust tableLookup + | otherwise = fromInteger (round $ findValue zvalue 0.00005 (0.5 - alpha) * 1000) / 1000.0 + where + tableLookup = zTable alpha + +-- | Utils for computing t-/z-/p-values +zvalue :: Double -> Double +zvalue x = (1.0 / sqrt (2.0 * pi)) * exp 1 ** negate (x**2/2) + +tvalue :: Double -> (Double -> Double) +tvalue df x = 1 / (sqrt df * beta 0.5 (df / 2.0)) * (1.0 + (x*x) / df) ** negate ((df + 1.0) / 2.0) + +-- Utility functions for finding the t-value +-- Source: https://wiki.haskell.org/index.php?title=Gamma_and_Beta_function +cof :: [Double] +cof = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953] + +ser :: Double +ser = 1.000000000190015 + +gammaln :: Double -> Double +gammaln xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5) + ser' = foldl (+) ser $ map (\(y,c) -> c/(xx+y)) $ zip [1..] cof + in -tmp' + log(2.5066282746310005 * ser' / xx) + +beta :: Double -> Double -> Double +beta z w = exp (gammaln z + gammaln w - gammaln (z+w)) + +findValue :: (Double -> Double) -> Double -> Double -> Double +findValue f stepSize target = fst $ until (\(_, x) -> x >= target) (\(a, x) -> (a + stepSize, x + stepSize * f a)) (0, 0) + +findValue' :: (Double -> Double) -> Double -> Double -> Double +findValue' f stepSize target = snd $ until (\(a, _) -> a >= target) (\(a, x) -> (a + stepSize, x + stepSize * f a)) (0, 0) + + +-- Compute the P-value +computePValue :: Monad m => TestType -> Sided -> Double -> Maybe Double -> m Double +computePValue test sided ts (Just df) + | isTTest test = return $ computePValueT sided ts df + | otherwise = fail "cannot compute p-value" +computePValue ZTest sided ts _ = + return $ computePValueZ sided ts +computePValue _ _ _ _ = + fail "cannot compute p-value" + +computePValueT :: Sided -> Double -> Double -> Double +computePValueT TwoSided ts df = (0.5 - pvalueT df (abs ts)) * 2 +computePValueT LeftSided ts df | ts < 0 = 0.5 - pvalueT df (abs ts) + | otherwise = 0.5 + pvalueT df ts +computePValueT RightSided ts df | ts < 0 = 1.0 - (0.5 - pvalueT df (abs ts)) + | otherwise = 0.5 - pvalueT df ts + +pvalueT :: Double -> Double -> Double +pvalueT df testStatistic | testStatistic > 4.0 = 0.5 -- Performance gain, assuming precision is not needed for exceptionally high values + | otherwise = findValue' (tvalue df) 0.00005 testStatistic + +computePValueZ :: Sided -> Double -> Double +computePValueZ TwoSided ts = (0.5 - pvalueZ (abs ts)) * 2 +computePValueZ LeftSided ts | ts < 0 = 0.5 - pvalueZ (abs ts) + | otherwise = 0.5 + pvalueZ ts +computePValueZ RightSided ts | ts < 0 = 1.0 - (0.5 - pvalueZ (abs ts)) + | otherwise = 0.5 - pvalueZ ts + +pvalueZ :: Double -> Double +pvalueZ testStatistic | testStatistic > 4.0 = 0.5 -- Performance gain, assuming precision is not needed for exceptionally high values + | otherwise = findValue' zvalue 0.00005 testStatistic + +-- | Utility function +chooseTypeOfTest :: ComponentSet -> TypeOfTest +chooseTypeOfTest cs + | contains cs Correlation = TestCorrelation + | contains cs (Two SampleMean) = if contains cs (Two SampleSize) + then CompareMeans + else CompareMeansPaired + | contains cs (Two Proportion) = CompareProportions + | contains cs Proportion = TestProportion + | contains cs SampleMean = TestMean + | otherwise = TestMean diff --git a/src/Domain/Hypothesis/Constraints.hs b/src/Domain/Hypothesis/Constraints.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Hypothesis/Constraints.hs @@ -0,0 +1,434 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +-- to do: split constraints in check against initial value, and inferred value +module Domain.Hypothesis.Constraints + ( hypothesisConstraints, equalDouble + ) where + +import Control.Monad +import Control.Applicative +import Data.Maybe +import Domain.Hypothesis.Common +import Domain.Hypothesis.Rules +import Domain.Math.Data.Relation +import Domain.Math.Numeric.Views +import Domain.Statistics.ComponentSet +import Ideas.Text.OpenMath.Dictionary.Arith1 +import Ideas.Common.Constraint +import Domain.Statistics.Views +import Ideas.Common.Library hiding (Predicate) +import Domain.Math.Expr.Data +import Prelude hiding (until) + +hypothesisConstraints :: [Constraint ComponentSet] +hypothesisConstraints = + [ checkHA + , checkH0 + , checkAlpha + , checkSided -- check this before checkCritical + , checkTest -- check this before checkCritical + , checkDf -- check this before checkCritical + , checkDfAnova -- check this before checkCritical + , checkGroups + , checkCritical, checkRejectionCritical + , checkSampleMean, checkSampleSdev, checkPopulationSdev + , checkPopulationMean, checkSampleSize, checkCorrelation + , checkStandardError + , checkTestFormula, checkTestValue + , checkPValue + -- checking conclusions should be last + , checkConclusionPValue, checkConclusionCritical, checkConclusionHypotheses + , checkObservedTotals, checkExpectedFrequencies + ] + +checkSampleMean, checkSampleSdev, checkPopulationSdev, + checkPopulationMean, checkSampleSize, checkCorrelation :: Constraint ComponentSet +checkSampleMean = checkForInitial "samplemean" SampleMean +checkSampleSdev = checkForInitial "samplesdev" SampleSdev +checkPopulationSdev = checkForInitial "populationsdev" PopulationSdev +checkPopulationMean = checkForInitial "populationmean" PopulationMean +checkSampleSize = checkForInitial "samplesize" SampleSize +checkCorrelation = checkForInitial "correlation" Correlation + +checkForInitial :: String -> ComponentId -> Constraint ComponentSet +checkForInitial s cid = makeConstraint ("check" # s) $ \cs -> do + relevance $ guard $ isInitial cid cs && isDerived cid cs + x <- get cid (initials cs) + y <- get cid (derived cs) + guard $ x == y + +checkAlpha :: Constraint ComponentSet +checkAlpha = makeConstraint "check.alpha" $ \cs -> do + relevance $ guard $ isDerived SignificanceLevel cs + if isInitial SignificanceLevel cs then do + -- initial alpha overrules the default value 0.05 + x <- get SignificanceLevel (initials cs) + y <- get SignificanceLevel (derived cs) + guard $ x == y + else do + -- use the default value 0.05 + alpha <- getExpr SignificanceLevel cs + val <- matchM doubleView alpha + guard $ val == 0.05 + +checkH0 :: Constraint ComponentSet +checkH0 = makeConstraint "check.h0" $ \cs -> do + relevance $ guard $ isDerived NullHypothesis cs + test <- inferTestChoice cs + case test of + ChiSquared -> do + h0 <- getExpr NullHypothesis cs + guard (h0 == independent) + _ -> do + h0 <- getRelation NullHypothesis cs + h0' <- h0FromHA <$> getRelation AlternativeHypothesis (initials cs) + unless (leftHandSide h0 == leftHandSide h0') $ + fail "parameter mismatch" + unless (rightHandSide h0 == rightHandSide h0') $ + fail "value mismatch" + unless (relationType h0 == relationType h0' || relationType h0 == EqualTo) $ + fail "sign mismatch" + +checkHA :: Constraint ComponentSet +checkHA = makeConstraint "check.ha" $ \cs -> do + relevance $ guard $ isDerived AlternativeHypothesis cs + test <- inferTestChoice cs + case test of + ChiSquared -> do + ha <- getExpr AlternativeHypothesis cs + guard (ha == dependent) + _ -> do + ha <- getRelation AlternativeHypothesis cs + ha' <- case get AlternativeHypothesis (initials cs) of + Just (CRelation x) -> return x + _ -> empty + unless (leftHandSide ha == leftHandSide ha') $ + fail "parameter mismatch" + unless (rightHandSide ha == rightHandSide ha') $ + fail "value mismatch" + unless (relationType ha == relationType ha') $ + fail "sign mismatch" + +checkConclusionPValue :: Constraint ComponentSet +checkConclusionPValue = makeConstraint "check.conclusion.p-value" $ \cs -> do + relevance $ guard $ isDerived ConclusionPValue cs + conc <- getRelation ConclusionPValue cs + let cs' = substitute cs + p <- matchM doubleView =<< getExpr PValue cs' + alpha <- matchM doubleView =<< getExpr SignificanceLevel cs' + guard $ if p <= alpha + then relationType conc == LessThanOrEqualTo + else relationType conc == GreaterThan + +checkSided :: Constraint ComponentSet +checkSided = makeConstraint "check.sided" $ \cs -> do + relevance $ guard $ isDerived Sidedness cs + -- TODO: make order of checking constraints explicit + -- && isNothing (isViolated checkHA cs) + sided <- getSided Sidedness cs + case getTestType TestChoice (initials cs) of + Just Anova -> + unless (sided == RightSided) (fail "anova") + Just ChiSquared -> + unless (sided == RightSided) (fail "chi-squared") + _ -> do + ha <- getRelation AlternativeHypothesis cs + guard $ sided == sidedFromHA ha + +checkTest :: Constraint ComponentSet +checkTest = makeConstraint "check.test" $ \cs -> do + relevance $ guard $ isDerived TestChoice cs + if isInitial TestChoice cs then do + -- initial testType overrules the validTests function + x <- get TestChoice (initials cs) + y <- get TestChoice (derived cs) + guard $ x == y + else do + test <- getTestType TestChoice cs + guard $ test `elem` validTests cs + +checkTestValue, checkTestFormula :: Constraint ComponentSet +checkTestFormula = checkTestValueFor "test-formula" TestFormula +checkTestValue = checkTestValueFor "test-value" TestValue + +checkTestValueFor :: String -> ComponentId -> Constraint ComponentSet +checkTestValueFor s cid = makeConstraint ("check" # s) $ \cs -> do + relevance $ guard $ isDerived cid cs + rel <- getRelation cid cs + let testTypes = inferTestChoices cs + testType <- + case leftHandSide rel of + Var "z" | ZTest `elem` testTypes -> return ZTest + Var "r" | RPearson `elem` testTypes -> return RPearson + Var "F" | Anova `elem` testTypes -> return Anova + Var "chisq" | ChiSquared `elem` testTypes -> return ChiSquared + Var "t" -> + case filter isTTest testTypes of + [] -> fail "test mismatch" + hd:_ -> return hd + _ -> fail "test mismatch" + let tv = rightHandSide rel + sub = getSubstitution cs + tv' <- matchM doubleView (sub |-> tv) + -- special case: when checking test value, first check wether an initial + -- test value is present + case getRhsExpr TestValue (initials cs) of + Just initialTv | cid == TestValue -> do + val <- matchM doubleView initialTv + testDoubles tv' val + _ -> do + target <- case getRhsExpr TestFormula (initials cs) of + Just rhs -> return rhs + _ | testType == ChiSquared -> rightHandSide <$> chiSquaredTestValue cs + _ -> testFormulaFromTest testType (chooseTypeOfTest cs) + target' <- matchM doubleView (sub |-> target) + testDoubles tv' target' + +checkDf :: Constraint ComponentSet +checkDf = makeConstraint "check.df" $ \cs -> do + relevance $ guard $ isDerived Df cs + let cs' = substitute cs + df <- matchM doubleView =<< getExpr Df cs' + case getExpr Df (initials cs) of + Just initialDf -> do + val <- matchM doubleView initialDf + guard $ equalDouble df val + _ -> do + test <- inferTestChoice cs + target <- case test of + ChiSquared -> toExpr <$> chiSquaredDf cs + _ -> degreesOfFreedomFromTest test (chooseTypeOfTest cs) + let sub = getSubstitution cs + target' <- matchM doubleView (sub |-> target) + guard $ equalDouble df target' + +checkDfAnova :: Constraint ComponentSet +checkDfAnova = makeConstraint "check.df-anova" $ \cs -> do + relevance $ guard $ isDerived DfBetween cs || isDerived DfWithin cs + let cs' = substitute cs + (between, within) <- inferDfBetweenWithin cs' + -- check df-between + dfB <- getExpr DfBetween (derived cs) >>= matchM doubleView + unless (equalDouble dfB between) $ fail "df-between" + -- check df-within + dfW <- getExpr DfWithin (derived cs) >>= matchM doubleView + unless (equalDouble dfW within) $ fail "df-within" + +checkGroups :: Constraint ComponentSet +checkGroups = makeConstraint "check.groups" $ \cs -> do + relevance $ guard $ isDerived Groups cs + groups <- getExpr Groups (derived cs) >>= matchM doubleView + -- for now, only support for two groups + guard (groups == 2) + +checkCritical :: Constraint ComponentSet +checkCritical = makeConstraint "check.critical" $ \cs -> do + relevance $ guard $ isDerived Critical cs + unless (derived cs `contains` AlternativeHypothesis) $ + fail "alternative hypothesis missing" + let cs' = substitute cs + critical <- getRelation Critical cs + criticalValue <- matchM doubleView (rightHandSide critical) + let testTypes = inferTestChoices cs + testType <- + case leftHandSide critical of + Var "zcrit" | ZTest `elem` testTypes -> return ZTest + Var "rcrit" | RPearson `elem` testTypes -> return RPearson + Var "Fcrit" | Anova `elem` testTypes -> return Anova + Var "chicrit" | ChiSquared `elem` testTypes -> return ChiSquared + Var "tcrit" -> + case filter isTTest testTypes of + [] -> fail "test mismatch" + hd:_ -> return hd + this -> fail $ "test mismatch " ++ show this + sided <- inferSidedness cs + alpha <- matchM doubleView =<< getExpr SignificanceLevel cs + value' <- + case testType of + Anova -> do + (dfBetween, dfWithin) <- inferDfBetweenWithin cs + computeCriticalF dfBetween dfWithin alpha + ChiSquared -> do + df <- chiSquaredDf cs + computeCriticalChi sided alpha (fromIntegral df) + _ -> do + -- TO DO: use computeCriticalF for Anova + let cs'' = case inferDf cs' of + Just df -> substitute (append Df (CExpr df) cs') + Nothing -> cs' + let df = matchM doubleView =<< getExpr Df cs'' + computeCritical testType sided alpha df + testDoubles value' criticalValue + +checkRejectionCritical :: Constraint ComponentSet +checkRejectionCritical = makeConstraint "check.rejectioncritical" $ \cs -> do + relevance $ guard $ isDerived RejectionCritical cs + unless (derived cs `contains` AlternativeHypothesis) $ + fail "alternative hypothesis missing" + crit <- getRelation RejectionCritical cs + sided <- inferSidedness cs + + let testTypes = inferTestChoices cs + (testType, rel) <- + case withoutAbs (leftHandSide crit) of + Var "z" | ZTest `elem` testTypes -> + return (ZTest, sidedRelation sided (Var "z") (Var "zcrit")) + Var "r" | RPearson `elem` testTypes -> + return (RPearson, sidedRelation sided (Var "r") (Var "rcrit")) + Var "F" | Anova `elem` testTypes -> + return (Anova, sidedRelation sided (Var "F") (Var "Fcrit")) + expr | expr == chiSquared && ChiSquared `elem` testTypes -> + return (ChiSquared, sidedRelation sided chiSquared (Var "chicrit")) + Var "t" -> + case filter isTTest testTypes of + [] -> fail "test mismatch" + hd:_ -> return (hd, sidedRelation sided (Var "t") (Var "tcrit")) + expr -> fail $ "test mismatch" ++ show expr + + unless (withoutAbs (leftHandSide rel) == withoutAbs (leftHandSide crit) + && rightHandSide rel == rightHandSide crit) $ + fail "test mismatch" + unless (leftHandSide rel == leftHandSide crit + && relationType rel == relationType crit) $ + if testType == Anova + then fail "anova mismatch" + else if testType == ChiSquared + then fail "chi-squared mismatch" + else fail "sidedness mismatch" + +withoutAbs :: Expr -> Expr +withoutAbs (Sym s [a]) | s == newSymbol absSymbol = a +withoutAbs expr = expr + +checkPValue :: Constraint ComponentSet +checkPValue = makeConstraint "check.p-value" $ \cs -> do + relevance $ guard $ isDerived PValue cs + unless (derived cs `contains` AlternativeHypothesis) $ + fail "alternative hypothesis missing" + let cs' = substitute cs + t <- matchM doubleView =<< fmap rightHandSide (getRelation TestValue cs') + test <- inferTestChoice cs + pvalue <- matchM doubleView =<< getExpr PValue cs + when (pvalue < 0 || pvalue > 1) $ + fail "value not a probability" + sided <- inferSidedness cs + let cs'' = case inferDf cs' of + Just df -> substitute (append Df (CExpr df) cs') + Nothing -> cs' + let df = matchM doubleView =<< getExpr Df cs'' + pvalue' <- computePValue test sided t df + -- compare with 3 decimals + testDoublesWith 0.00055 pvalue' pvalue $ + if derived cs `contains` TestValue + then empty + else fail "TestValue missing" + +checkConclusionCritical :: Constraint ComponentSet +checkConclusionCritical = makeConstraint "check.conclusion-critical" $ \cs -> do + relevance $ guard $ isDerived ConclusionCritical cs + rej <- inferRejectionCritical cs + -- to do: rejection critical is added to component set only to get the substituted relation + let cs' = substitute (append RejectionCritical (CRelation rej) cs) + rejection <- getRelation RejectionCritical cs' + lhs <- matchM doubleView $ leftHandSide rejection + rhs <- matchM doubleView $ rightHandSide rejection + let result = eval (relationType rejection) lhs rhs + conclusion <- getConclusion ConclusionCritical cs + guard $ result == conclusion + +checkConclusionHypotheses :: Constraint ComponentSet +checkConclusionHypotheses = makeConstraint "check.conclusion-hypotheses" $ \cs -> do + relevance $ guard $ isDerived ConclusionHypotheses cs + concl <- getRejectionHypotheses ConclusionHypotheses cs + b1 <- case fmap relationType (inferConclusionPValue cs) of + Just LessThanOrEqualTo + | concl == RejectH0 -> return True + Just GreaterThan + | concl == DontRejectH0 -> return True + Just _ -> fail "conclusion mismatch pvalue" + Nothing -> return False + b2 <- case inferConclusionCritical cs of + Just False + | concl == DontRejectH0 -> return True + | concl `elem` [AcceptH0, RejectH1] -> fail "convention DontRejectH0" + Just True + | concl == RejectH0 -> return True + | concl `elem` [AcceptH1, DontRejectH1] -> fail "convention RejectH0" + Just _ -> fail "conclusion mismatch critical" + Nothing -> return False + guard (b1 || b2) + +checkStandardError :: Constraint ComponentSet +checkStandardError = makeConstraint "check.standard-error" $ \cs -> do + relevance $ guard $ isDerived StandardError cs + rel <- getRelation StandardError cs + n <- matchM doubleView =<< getExpr SampleSize cs + case (leftHandSide rel, getRhsExpr PopulationSdev cs, getRhsExpr SampleSdev cs) of + (Var "sigmaM", Just expr, Nothing) -> do + psdev <- matchM doubleView expr + se <- matchM doubleView (rightHandSide rel) + testDoubles se (psdev/sqrt n) + (Var "SEM", Nothing, Just expr) -> do + sdev <- matchM doubleView expr + se <- matchM doubleView (rightHandSide rel) + testDoubles se (sdev/sqrt n) + _ -> fail "standard error mismatch" + +checkObservedTotals :: Constraint ComponentSet +checkObservedTotals = makeConstraint "check.observed-totals" $ \cs -> do + relevance $ guard $ + all (derived cs `contains`) [ObservedRowTotals, ObservedColumnTotals, ObservedTotal] + rowTotals <- getExpr ObservedRowTotals cs >>= fromExpr + columnTotals <- getExpr ObservedColumnTotals cs >>= fromExpr + total <- getExpr ObservedTotal cs >>= fromExpr + table <- getTable ObservedFrequencies cs + let (expRowTotals, expColumnTotals, expTotal) = computeTotals table + unless (rowTotals == expRowTotals) $ fail "rows" + unless (columnTotals == expColumnTotals) $ fail "columns" + unless (total == expTotal) $ fail "total" + +checkExpectedFrequencies :: Constraint ComponentSet +checkExpectedFrequencies = makeConstraint "check.expected-frequencies" $ \cs -> do + relevance $ guard $ derived cs `contains` ExpectedFrequencies + expectedFrequencies <- getExpr ExpectedFrequencies cs >>= matchDoubleTable + observed <- getTable ObservedFrequencies cs + let totals = computeTotals observed + computedFrequencies = computeExpectedFrequencies totals + guard (equalTableDouble expectedFrequencies computedFrequencies) + +matchDoubleTable :: MonadPlus m => Expr -> m [[Double]] +matchDoubleTable = fromExpr >=> mapM (mapM (matchM doubleView)) + +equalTableDouble :: [[Double]] -> [[Double]] -> Bool +equalTableDouble = equalListBy (equalListBy equalDouble) + +equalListBy :: (a -> a -> Bool) -> [a] -> [a] -> Bool +equalListBy eq xs ys = length xs == length ys && and (zipWith eq xs ys) + +---------------------------------------------------------- + +testDoubles :: Double -> Double -> Result () +testDoubles value target = testDoublesWith defaultDelta value target empty + +testDoublesWith :: Double -> Double -> Double -> Result () -> Result () +testDoublesWith delta value target resultNotEqual = + unless (equalDoubleWith delta value target) $ + if equalDoubleWith (delta*10) value target + then fail "almost equal" + else resultNotEqual + +equalDouble :: Double -> Double -> Bool +equalDouble = equalDoubleWith defaultDelta + +equalDoubleWith :: Double -> Double -> Double -> Bool +equalDoubleWith delta x y = abs (x - y) < delta + +defaultDelta :: Double +defaultDelta = 0.0055 + +-- afronden op 3 cijfers achter de komma +-- rounded :: Double -> Double +-- rounded x = fromInteger (round (x * 1000)) / 1000.0 diff --git a/src/Domain/Hypothesis/Examples.hs b/src/Domain/Hypothesis/Examples.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Hypothesis/Examples.hs @@ -0,0 +1,450 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Hypothesis.Examples + {- ( hypothesisExamples + , ex_dwo2 + , ex_dwo3 + , ex_simple + ) -} where + +import Domain.Math.Data.Relation +import Domain.Math.Expr +import Domain.Statistics.ComponentSet +import Ideas.Common.Exercise +import Prelude hiding ((^)) + +hypothesisExamples :: Examples ComponentSet +hypothesisExamples = examplesWithDifficulty $ + [ (Easy, ex_simple) + , (Easy, ex_omi1) + , (Easy, ex_omi2) + , (Easy, ex_omi3) + , (Medium, ex_econ) + , (Medium, ex_dwo1) + , (Medium, ex_dwo2) + , (Medium, ex_dwo3) + , (Medium, ex_dwo4) + , (Medium, ex_dwo5) + , (Medium, ex_dwo6) + , (Medium, ex_dwo7) + , (Medium, sietske_1) + , (Medium, sietske_2) + , (Medium, sietske_3) + , (Difficult, ex_bio') + ] ++ + [ (Medium, cs) | cs <- [mayPilotA, mayPilotB, mayPilotC, mayPilotD] ] + +cExprDouble :: Double -> Component +cExprDouble = CExpr . fromDouble + +-------------------------------------------------------------------------------- +-- * Pilot exercises (May 2017) * -- +-------------------------------------------------------------------------------- + +{- Montarello and Martins (2005) discovered that children from the 7th group + were able to solve more difficult maths exercises if very simple exercises + were added to the exam. To research this phenomenon a researcher picks a + standardised exam of which the scores have a gaussian distribution with + + mu = 50, + sigma = 8, + + and adds very easy exercises to it. They give the exam to 40 children. The + mean outcome of the test is M = 52. Is this result enough to conclude that + adding the easy exercises improves the childrens prestation? Assume a + significance level of alpha = 0.05 -} + +-- Note: Claim assumed that `improvement' means mu_after - mu_before > 0 +mayPilotA :: ComponentSet +mayPilotA = initialSet + [ -- ( Claim , CRelation $ Var "mu" .>. fromDouble 50 ) + ( AlternativeHypothesis, CRelation $ Var "mu" .>. fromDouble 50) + , ( SampleMean , cExprDouble 52 ) + , ( PopulationMean , cExprDouble 50 ) -- I think that different means will + , ( PopulationSdev , cExprDouble 8 ) -- cause all sorts of trouble with + , ( SampleSize , CExpr 40 ) -- the strategy. + , ( SignificanceLevel , cExprDouble 0.05 ) ] + +{- It has been showed by Evans, Pelham, Smith et al. (2001) that the Ritalin + medicin improves the attention span of children with ADHD and that their + school results improve. + + A researcher draws a sample of N = 16 children that have been diagnosed with + ADHD to demonstrate that the medicine truly works. Their attention span is + measured before and after the medicine has been used. The results show an + average improvement of the attention span of M = 3.13 minutes with a standard + deviation of s = 5.39. Are the results of this experiment enough to conclude + that the the medicine improves the attention span of children with ADHD? + Assume a significance level of alpha = 0.01. -} + +-- Note: mu_d means mu_after - mu_before +mayPilotB :: ComponentSet +mayPilotB = initialSet + [ -- ( Claim , CRelation $ Var "mu" .>. fromDouble 0) + ( AlternativeHypothesis, CRelation $ Var "mu" .>. fromDouble 0) + , ( SampleSize , CExpr 16 ) + , ( SampleMean , cExprDouble 3.13 ) + , ( SampleSdev , cExprDouble 5.39 ) + , ( SignificanceLevel , cExprDouble 0.01 ) + , ( TestChoice , CChoice (TestType TTestPaired)) + ] + +{- A researcher compares two treatments during a design with independent groups. + Group A has 18 participants and group V has 16 participants. The means for the + groups are found to be M = 12.14 for group A and M = 16.17 for group B. The + pooled variance s^2 is approximated to be 7.74. The researcher wants to know + if the treatments will leaad to a different result. Assume a significance + level of alpha = 0.05. -} + +mayPilotC :: ComponentSet +mayPilotC = initialSet + [ -- ( Claim , CRelation $ Var "mu1" ./=. Var "mu2") + ( AlternativeHypothesis, CRelation $ Var "mu1" ./=. Var "mu2") + , ( One SampleSize , CExpr 18 ) + , ( Two SampleSize , CExpr 16 ) + , ( One SampleMean , cExprDouble 12.14 ) + , ( Two SampleMean , cExprDouble 16.70 ) + , ( PooledVariance , cExprDouble 7.74 ) + , ( SignificanceLevel , cExprDouble 0.05 ) + ] + +{- The mean annual income for a randomised sample of 51 househoulds turns out to + be 23 613 euros. The standard deviation of the sample is 658 euros. You may + assume that annual income is normally distributed. Research the claim that + mean annual income in the population is lower than 24000. Assume a + significance level of alpha = 0.10. -} + +mayPilotD :: ComponentSet +mayPilotD = initialSet + [ -- ( Claim , CRelation $ Var "mu" .<. fromDouble 24000) + ( AlternativeHypothesis, CRelation $ Var "mu" .<. fromDouble 24000) + , ( SampleSize , CExpr 51 ) + , ( SampleMean , cExprDouble 23613 ) + , ( PopulationMean , cExprDouble 24000 ) + , ( SampleSdev , cExprDouble 658 ) + , ( SignificanceLevel , cExprDouble 0.10 ) + ] + +-- Example from an exam +ex_econ :: ComponentSet +ex_econ = initialSet + [ (SampleSize, CExpr 61) + , (SampleMean, CExpr $ fromDouble 21023.0) + , (SampleSdev, CExpr $ fromDouble 634.0) + , (AlternativeHypothesis, CRelation $ Var "mu" .<=. 22000) + ] + +-- Example for an hypothesis test +ex_simple :: ComponentSet +ex_simple = initialSet [ + (AlternativeHypothesis, CRelation $ Var "mu" .>=. fromDouble 98.6), + (SampleSize, CExpr 101), + (SampleMean, CExpr $ fromDouble 98.9), + (SampleSdev, CExpr $ fromDouble 0.6)] + +-- Examples from the OMI course +ex_omi1 :: ComponentSet +ex_omi1 = initialSet [ + (AlternativeHypothesis, CRelation $ Var "mu" ./=. fromDouble 20.0), + (SampleSize, CExpr 10), + (SampleMean, CExpr $ fromDouble 22.3), + (SampleSdev, CExpr $ fromDouble 3.65) + ] + +ex_omi2 :: ComponentSet +ex_omi2 = initialSet [ + (AlternativeHypothesis, CRelation $ Var "mu" .>=. fromDouble 20.0), + (SampleSize, CExpr 10), + (SampleMean, CExpr $ fromDouble 22.3), + (SampleSdev, CExpr $ fromDouble 3.65) + ] + +ex_omi3 :: ComponentSet +ex_omi3 = initialSet [ + (AlternativeHypothesis, CRelation $ Var "mu1" .==. Var "mu2"), + (One SampleSize, CExpr 46), + (Two SampleSize, CExpr 56), + (One SampleMean, CExpr $ fromDouble 3.8698), + (Two SampleMean, CExpr $ fromDouble 4.5819), + (One SampleSdev, CExpr $ fromDouble 1.6714), + (Two SampleSdev, CExpr $ fromDouble 1.5216), + (PooledVariance, CExpr $ fromDouble 2.531) + ] + +-- Example from an exam simplified +ex_bio' :: ComponentSet +ex_bio' = initialSet [ + (AlternativeHypothesis, CRelation $ Var "mu1" .==. Var "mu2"), + (One SampleSize, CExpr 9), + (Two SampleSize, CExpr 9), + (One SampleMean, CExpr $ fromDouble 19.189), + (Two SampleMean, CExpr $ fromDouble 28.067), + (One SampleSdev, CExpr $ fromDouble 6.36), + (Two SampleSdev, CExpr $ fromDouble 6.12), + (PooledVariance, CExpr $ fromDouble 38.843)] + +-- Examples from the DWO +ex_dwo1 :: ComponentSet +ex_dwo1 = ex_econ + +ex_dwo2 :: ComponentSet +ex_dwo2 = initialSet [ + (SampleSize, CExpr 68), + (SampleMean, CExpr $ fromDouble 23035.0), + (SampleSdev, CExpr $ fromDouble 658.0), + (AlternativeHypothesis, CRelation $ Var "mu" .>=. 24000), + (SignificanceLevel, CExpr $ fromDouble 0.10)] + +ex_dwo3 :: ComponentSet +ex_dwo3 = initialSet [ + (SampleSize, CExpr 61), + (SampleMean, CExpr $ fromDouble 21023.0), + (SampleSdev, CExpr $ fromDouble 634.0), + (AlternativeHypothesis, CRelation $ Var "p" ./=. fromDouble 0.25), + (Proportion, CExpr $ fromDouble 0.21) +-- fix me (TestStatistic "t", CExpr $ (Var "p" - Var "p0") / sqrt ((Var "p0") * (1.0 - Var "p0") / (Var "n"))) + ] + +ex_dwo4 :: ComponentSet +ex_dwo4 = initialSet [ + (SampleSize, CExpr 72), + (SampleMean, CExpr $ fromDouble 24061.0), + (SampleSdev, CExpr $ fromDouble 663.0), + (AlternativeHypothesis, CRelation $ Var "p" ./=. fromDouble 0.3), + (Proportion, CExpr $ fromDouble 0.22), + --(TestStatistic "t", CExpr $ (Var "p" - Var "p0") / sqrt ((Var "p0") * (1.0 - Var "p0") / (Var "n"))), + (SignificanceLevel, CExpr $ fromDouble 0.01)] + +ex_dwo5 :: ComponentSet +ex_dwo5 = initialSet [ + (SampleSize, CExpr 100), + (SampleMean, CExpr $ fromDouble 24061.0), + (SampleSdev, CExpr $ fromDouble 663.0), + (AlternativeHypothesis, CRelation $ Var "p1" .>=. Var "p2"), + (One Proportion, CExpr $ fromDouble 0.3), + (Two Proportion, CExpr $ fromDouble 0.1), + (Other "p0", CExpr $ (Var "p1" + Var "p2") / 2), + (Other "d0", CExpr $ fromDouble 0.0) +-- fix me (TestStatistic "t", CExpr $ (Var "p1" - Var "p2" - Var "d0") / (sqrt $ (Var "p0") * (1.0 - Var "p0") / ((Var "n") / 2))) + ] + +ex_dwo6 :: ComponentSet +ex_dwo6 = initialSet [ + (SampleSize, CExpr 100), + (SampleMean, CExpr $ fromDouble 24061.0), + (SampleSdev, CExpr $ fromDouble 663.0), + (AlternativeHypothesis, CRelation $ Var "p1" ./=. Var "p2"), + (One Proportion, CExpr $ fromDouble 0.3), + (Two Proportion, CExpr $ fromDouble 0.1), + (Other "p0", CExpr $ (Var "p1" + Var "p2") / 2), + (Other "d0", CExpr $ fromDouble 0.0) +-- fix me (TestStatistic "t", CExpr $ (Var "p1" - Var "p2" - Var "d0") / (sqrt $ (Var "p0") * (1.0 - Var "p0") / ((Var "n") / 2))) + ] + +ex_dwo7 :: ComponentSet +ex_dwo7 = initialSet [ + (SampleSize, CExpr 490), + (SampleMean, CExpr $ fromDouble 83.0), + (Correlation, CExpr $ fromDouble (-0.557)), + (Df, CExpr $ Var "n" - 2), + (AlternativeHypothesis, CRelation $ Var "r" .>=. fromDouble 0.0) + ] + + +{- "Average income" example + - + - Steekproefgrootte n = 61 + - Steekproefgemiddelde m = 21023 + - Steekproefstandaardafwijking s = 634 + - “Je mag aannemen dat inkomen normaal is verdeeld” + - Significantieniveau α = 0,05. + - Claim: Het populatiegemiddelde μ is kleiner dan 22000. + -} + +sietske_1 :: ComponentSet +sietske_1 = initialSet [ + (SampleSize, CExpr 61), + (SampleMean, CExpr 21023), + (SampleSdev, CExpr 634), + (SignificanceLevel, CExpr $ fromDouble 0.05), + (AlternativeHypothesis, CRelation $ Var "mu" .<=. 22000) + ] + +{- "Car factory" example + - + - Steekproefgrootte: n_1 = 100, n_2 = 100 + - Steekproefproportie 1: p_A = 0.3 + - Steekproefproportie 2: p_B = 0.1 + - Significantieniveau: alpha = 0.05. + - p_0 = (p_A + p_B) / 2 + - Verschil tussen populatieproporties onder nulhypothese D_0 = 0 + - Toetsingsgrootheid + - Claim: Populatieproportie P_A is groter dan populatieproportie P_B + -} +sietske_2 :: ComponentSet +sietske_2 = initialSet + [ (SampleSize , CExpr 100) + , (One Proportion , CExpr $ fromDouble 0.3) + , (Two Proportion , CExpr $ fromDouble 0.1) + , (SignificanceLevel, CExpr $ fromDouble 0.05) + , (Other "p0" , CExpr $ (Var "p1" + Var "p2") / 2) + , (Other "d0" , CExpr $ fromDouble 0.0) +-- fix me +-- , (TestStatistic "t", +-- CExpr $ (Var "p1" - Var "p2" - Var "d0") / (sqrt $ (Var "p0") * +-- (1.0 - Var "p0") / ((Var "n") / 2))) + , (AlternativeHypothesis , CRelation $ Var "p1" .>=. Var "p2") + ] + +{- "Shop owner" example + - + - Steekproefgrootte n_x=20 en n_y=20 + - Steekproefgemiddelde 1: = 85 + - Steekproefgemiddelde 2: = 63 + - Steekproefstandaardafwijking 1: s_x = 11 + - Steekproefstandaardafwijking 2: s_y = 11 + - Aanname dat normaal verdeeld + - Verschil tussen populatiegrootheden onder nulhypothese D_0=0 + - Toetsingsgrootheid + - Significantieniveau α = 0,05. + - Aantal vrijheidsgraden DF=n_x+n_y-2 + - Claim: Populatiegemiddelde mu_A is ongelijk aan populatiegemiddelde mu_B. + -} +sietske_3 :: ComponentSet +sietske_3 = initialSet [ + (SampleSize, CExpr 100), + (One SampleMean, CExpr 85), + (Two SampleMean, CExpr 63), + (One SampleSdev, CExpr 11), + (Two SampleSdev, CExpr 11), +-- fix me(TestStatistic "t", CExpr $ (mean1 - mean2 - (Var "d0")) / (sqrt((sd1 ^ 2 / (Var "n")) + (sd2 ^ 2 / (Var "n")))) ), + (SignificanceLevel, CExpr $ fromDouble 0.05), + -- (Df, CExpr $ (Var "n1) + (Var "n2") - 2), + (AlternativeHypothesis, CRelation $ Var "muA" ./=. Var "muB") + ] + {- + where + mean1 = toExpr (One SampleMean) + mean2 = toExpr (Two SampleMean) + sd1 = toExpr (One SampleSdev) + sd2 = toExpr (Two SampleSdev) -} + +-------------------------------------------------------------------------------- +-- Opgaven Sociale Wetenschappen + +opgavenSW :: [ComponentSet] +opgavenSW = [ -- december 2017 + opgave_3_4, opgave_3_6, opgave_4_10, opgave_4_11, opgave_5_3, opgave_5_6 + -- februari 2018 + , opgave_1_5, opgave_3_5, opgave_5_4 + ] + +-------------------------------------------------------------------------------- +-- Opgaven Sociale Wetenschappen (december 2017) + +opgave_3_4 :: ComponentSet +opgave_3_4 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "mu" .>. 100) + , (PopulationSdev, CRelation $ Var "sigma" .==. 18) + , (SampleMean, CExpr 104) + , (SampleSize, CExpr 36) + , (SignificanceLevel, CExpr 0.01) + ] + +opgave_3_6 :: ComponentSet +opgave_3_6 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "mu" .<. 4.9) + , (PopulationSdev, CRelation $ Var "sigma" .==. 0.84) + , (SampleMean, CExpr 4.4) + , (SampleSize, CExpr 16) + , (SignificanceLevel, CExpr 0.05) + ] + +opgave_4_10 :: ComponentSet +opgave_4_10 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "mu" ./=. 0) + , (SampleSdev, CRelation $ Var "s" .==. 1.50) + , (SampleMean, CExpr 1.28) + , (SampleSize, CExpr 25) + , (TestChoice, CChoice $ TestType TTestPaired) + , (SignificanceLevel, CExpr 0.05) + ] + +opgave_4_11 :: ComponentSet +opgave_4_11 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "mu" ./=. 0) + , (SampleSdev, CRelation $ Var "s" .==. 2.45) + , (SampleMean, CExpr 2) + , (SampleSize, CExpr 6) + , (TestChoice, CChoice $ TestType TTestPaired) + , (SignificanceLevel, CExpr 0.05) + ] + +opgave_5_3 :: ComponentSet +opgave_5_3 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "mu1" ./=. Var "mu2") + , (Df, CExpr 30) + , (TestValue, CRelation $ Var "t" .==. 2.085) + , (TestChoice, CChoice $ TestType TTestTwo) + , (SignificanceLevel, CExpr 0.05) + ] + +opgave_5_6 :: ComponentSet +opgave_5_6 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "mu1" ./=. Var "mu2") + , (Df, CExpr 119.50) + , (TestValue, CRelation $ Var "t" .==. 3.379) + , (TestChoice, CChoice $ TestType TTestTwo) + , (SignificanceLevel, CExpr 0.05) + ] + +-------------------------------------------------------------------------------- +-- Opgaven Sociale Wetenschappen (februari 2018) + +opgave_1_5 :: ComponentSet +opgave_1_5 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "rho" ./=. 0) + , (SignificanceLevel, CExpr 0.01) + , (TestValue, CRelation $ Var "r" .==. 0.835) + , (SampleSize, CExpr 9) + , (TestChoice, CChoice $ TestType RPearson) + ] + +opgave_3_5 :: ComponentSet +opgave_3_5 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "mu1" ./=. Var "mu2") + , (TestValue, CRelation $ Var "F" .==. 4.00) + , (SignificanceLevel, CExpr 0.05) + , (SampleSize, CExpr 40) + , (TestChoice, CChoice $ TestType Anova) + ] + +opgave_5_4 :: ComponentSet +opgave_5_4 = initialSet + [ (SignificanceLevel, CExpr 0.01) + , (TestChoice, CChoice $ TestType ChiSquared) + , (ObservedFrequencies, CExpr $ toExpr [[18 :: Int, 4, 2], [4, 17, 15]]) + ] + +-------------------------------------------------------------------------------- +-- Opgaven Economie (februari 2018) + +opgavenEconomie :: [ComponentSet] +opgavenEconomie = [opgave_4_18] + +opgave_4_18 :: ComponentSet +opgave_4_18 = initialSet + [ (AlternativeHypothesis, CRelation $ Var "mu1" ./=. Var "mu2") + , (SignificanceLevel, CExpr 0.05) + , (One SampleSize, CExpr 20) + , (Two SampleSize, CExpr 20) + , (One SampleMean, CExpr 85) + , (Two SampleMean, CExpr 63) + , (One SampleSdev, CExpr 11) + , (Two SampleSdev, CExpr 11) + , (TestFormula, CRelation $ Var "t" .==. + (Var "samplemean1" - Var "samplemean2") / sqrt ((Var "samplesd1" ^ 2) / Var "n1" + (Var "samplesd2" ^ 2) / Var "n2")) + ] diff --git a/src/Domain/Hypothesis/Exercises.hs b/src/Domain/Hypothesis/Exercises.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Hypothesis/Exercises.hs @@ -0,0 +1,178 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Hypothesis.Exercises (hypothesisExercise) where + +import Control.Monad +import qualified Data.Map as M +import Data.List +import Data.Maybe +import Domain.Math.Data.Relation +import Domain.Hypothesis.Examples +import Domain.Hypothesis.Constraints +import Domain.Hypothesis.BuggyRules +import Domain.Hypothesis.Rules +import Domain.Hypothesis.Strategies +import Domain.Statistics.ComponentSet +import Domain.Statistics.Parser +import Ideas.Common.Library hiding (Predicate) +import Domain.Math.Expr.Data +import Ideas.Encoding.OpenMathSupport +import Ideas.Utils.Uniplate +import Ideas.Text.XML hiding (children) +import Prelude hiding (until) + +-- Debug functions +css :: [ComponentSet] +css = mapMaybe fromContext $ concatMap terms $ + mapMaybe (defaultDerivation hypothesisExercise) (examplesAsList hypothesisExercise) + +_testPP :: IO () +_testPP = forM_ css $ \x -> + let s = prettyPrinter hypothesisExercise x in + case parser hypothesisExercise s of + Right y | x == y -> putChar '.' + | otherwise -> error $ show (x, y) + Left msg -> error $ s ++ "\n" ++ msg + +_testOM :: IO () +_testOM = forM_ css $ \x -> + case toOpenMath hypothesisExercise x of + Right omobj -> + case fromOpenMath hypothesisExercise omobj of + Just y | x == y -> putChar '.' + | otherwise -> error $ "Not the same:\n" ++ show x ++ "\n" ++ show y + Nothing -> error $ show omobj + Left msg -> error $ show x ++ "\n" ++ msg + +_see :: Int -> IO () +_see n = printDerivation hypothesisExercise (examplesAsList hypothesisExercise !! n) + +_testje :: IO () +_testje = printDerivation hypothesisExercise opgave_5_4 + +_save :: IO () +_save = do + let pilot = [("A", mayPilotA), ("B", mayPilotB), ("C", mayPilotC), ("D", mayPilotD)] + ex = hypothesisExercise + forM_ pilot $ \(n, cs) -> do + writeFile ("Task" ++ n ++ ".txt") $ showDerivation ex cs + writeFile ("Task" ++ n ++ "-OM.txt") $ showDerivationOM ex cs + +showDerivationOM :: Exercise a -> a -> String +showDerivationOM ex a = + case defaultDerivation ex a of + Just d -> show $ biMap fst f d + Nothing -> "no derivation" + where + errorXML = makeXML "error" mempty + f ctx = fromMaybe errorXML $ do + x <- fromContext ctx + omobj <- toOpenMath ex x + return (toXML omobj) + +------------------------------ + +hypothesisExercise :: Exercise ComponentSet +hypothesisExercise = emptyExercise + { exerciseId = describe "Hypothesis testing" $ newId "hypothesis" + , prettyPrinter = show + , status = Experimental + , parser = parseComponentSet + , strategy = liftToContext hypothesisStrategy + , extraRules = fmap liftToContext buggyRules + , ruleOrdering = ruleOrderingWith highPriorityRules + , examples = examplesWithDifficulty [ (Medium, cs) | cs <- opgavenSW ++ opgavenEconomie ] -- hypothesisExamples + , equivalence = withoutContext eqComponentSet + , similarity = withoutContext similarComponentSet + , constraints = map liftToContext hypothesisConstraints + , hasTypeable = useTypeable + , ready = predicate (hasConclusionHypotheses <&&> checkConstraints) + , hasTermView = Just termView + } + +highPriorityRules :: [Rule ComponentSet] +highPriorityRules = + [ addHypothesesRule + , lookupTValueRule, lookupZValueRule + , computePValueTTest, computePValueZTest + , addTestValueRule + , hypothesesConclusionCriticalRule, hypothesesConclusionPValueRule + -- buggy rules + , buggyTValueSided, buggyTValueTestValue -- Sietske: higher priority than corresponding buggy alpha rule + , buggyRValueSided, buggyRValueTestValue -- idem + , buggyZValueSided, buggyZValueTestValue + , buggyChiValueSided, buggyChiValueTestValue + ] + +---------------------------------------------------------- +-- Equivalence +eqComponentSet :: ComponentSet -> ComponentSet -> Bool +eqComponentSet x y = compareIntials + && checkConstraints x + && checkConstraints y + where + compareIntials = + M.fromList (toList (initials x)) == M.fromList (toList (initials y)) + +checkConstraints :: ComponentSet -> Bool +checkConstraints cs = + all (`checkConstraint` cs) hypothesisConstraints + +checkConstraint :: Constraint ComponentSet -> ComponentSet -> Bool +checkConstraint p cs + | not (isRelevant p cs) = True + | otherwise = isSatisfied p cs + +---------------------------------------------------------- +-- Similarity + +-- To do: not all component ids compare doubles with the same precision +-- (e.g. pvalue uses 3 decimals, default is 2) +similarComponentSet :: ComponentSet -> ComponentSet -> Bool +similarComponentSet cs1 cs2 = sorted cs1 `eqList` sorted cs2 + where + sorted = map f . sortOn fst . toList + where + f (n, CExpr e) = (n, CExpr (normalizeExpr e)) + f (n, CRelation r) = (n, CRelation (fmap normalizeExpr r)) + f (n, c) = (n, c) + + eqList :: [(ComponentId, Component)] -> [(ComponentId, Component)] -> Bool + eqList xs ys = length xs == length ys && all (\((a, x), (b, y)) -> a == b && x `eqComponent` y) (zip xs ys) + + eqComponent :: Component -> Component -> Bool + eqComponent (CExpr x) (CExpr y) = eqExpr x y + eqComponent (CRelation x) (CRelation y) = + eqExpr (leftHandSide x) (leftHandSide y) && eqExpr (rightHandSide x) (rightHandSide y) + && relationType x == relationType y + eqComponent x y = x == y + + eqExprs :: [Expr] -> [Expr] -> Bool + eqExprs [] [] = True + eqExprs (x:xs) (y:ys) = eqExpr x y && eqExprs xs ys + eqExprs _ _ = False + + eqExpr :: Expr -> Expr -> Bool + eqExpr (Number x) (Number y) = equalDouble x y + eqExpr x y = + case (getFunction x, getFunction y) of + (Just (f, xs), Just (g, ys)) -> f == g && xs `eqExprs` ys + _ -> x == y + + +-- we vouwen hier de definities uit zodat we formules kunnen herkennen +-- afronden op 2 decimalen: beter is om expressies component-gewijs te +-- vergelijken en daarbij de hulpfunctie 'equal' te gebruiken +normalizeExpr :: Expr -> Expr +normalizeExpr (Var "sigmaM") = Var "sigma" / sqrt (Var "n") +normalizeExpr (Var "SEM") = Var "s" / sqrt (Var "n") +normalizeExpr (Nat n) = Number (fromInteger n) +normalizeExpr expr = descend normalizeExpr expr + +-- | Test whether the conclusionhypothesis is present in the componentset. If +-- the predicate holds then he exercise is considered /done/. +hasConclusionHypotheses :: ComponentSet -> Bool +hasConclusionHypotheses cs = contains cs ConclusionHypotheses diff --git a/src/Domain/Hypothesis/Rules.hs b/src/Domain/Hypothesis/Rules.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Hypothesis/Rules.hs @@ -0,0 +1,666 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Hypothesis.Rules + ( addAlphaRule, addDfRule, addDfBetweenWithin + , addH0FromHARule, addH0FromHAEqualSignRule, addHARule + , addHypothesesRule, addHypothesesChiSquaredRule + , addObservedTotals, addExpectedFrequencies + , addConclusionPValueRule, addRejectionRule + , addTestFormulaRule, addTestValueRule + , chooseTTestPairedRule, chooseTTestRule, chooseTTestTwoRule + , chooseZTestRule, chooseRPearsonRule, chooseAnovaRule, chooseChiSquaredRule + , computePValueTTest, computePValueZTest + , criticalConclusionRule + , hypothesesConclusionCriticalRule, hypothesesConclusionPValueRule + , addStandardErrorSigma, addStandardErrorSD + , determineSided + , lookupTValueRule, lookupZValueRule, lookupRValueRule, lookupFValueRule + , lookupChiValueRule + ------------- + , inferSidedness, inferTestChoice, inferTestChoices, inferRejectionCritical + , inferConclusionCritical, inferConclusionPValue, inferDf, inferVar + , inferTestFormula, inferCriticalZWith, inferCriticalTWith + , inferCriticalRWith, inferCriticalFWith, inferCriticalChiWith + , inferTestValue, inferDfBetweenWithin, chiSquaredDf, chiSquaredTestValue + , computeTotals, computeExpectedFrequencies, getTable + ) where + +import Control.Monad +import Data.List +import Data.Maybe +import Domain.Hypothesis.Common +import Domain.Math.Data.Relation +import Domain.Math.Expr hiding ((.*.), (./.), (^)) +import Domain.Math.Numeric.Views +import Domain.Statistics.ComponentSet +import Domain.Statistics.Views +import Ideas.Common.Library + +---------------------------------------------------------- +-- Rules for determining the confidence level + +addAlphaRule :: Rule ComponentSet +addAlphaRule = + describe "Rule for adding the alpha component" . + makeRule "component.alpha" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` SignificanceLevel) + let alpha = case get SignificanceLevel (initials cs) of + Just (CExpr a) -> a + _ -> toExpr pickAlpha + return $ append SignificanceLevel (CExpr alpha) cs + +---------------------------------------------------------- +-- Rules for constructing the hypotheses + +determineSided :: Rule ComponentSet +determineSided = + describe "Rule for determine one-/two-sided testing" . + makeRule "component.sided" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` Sidedness) + sided <- + case getTestType TestChoice (initials cs) of + Just Anova -> return RightSided + _ -> do + ha <- getRelation AlternativeHypothesis cs + return (sidedFromHA ha) + return $ append Sidedness (CChoice $ Sided sided) cs + +addHypothesesRule :: Rule ComponentSet +addHypothesesRule = + describe "Add null hypothesis and alternative hypothesis, in one step" $ + makeRule "component.hypotheses" f + where + f :: ComponentSet -> [ComponentSet] + f = applyAll $ (addH0FromHARule ./. addH0FromHAEqualSignRule) .*. addHARule + +addH0FromHARule :: Rule ComponentSet +addH0FromHARule = + describe "Rule for adding the H0 component based on HA" . + makeRule "component.h0-from-ha" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` NullHypothesis) + ha <- getRelation AlternativeHypothesis cs + return . append NullHypothesis + (CRelation $ h0FromHA ha) $ cs + +addH0FromHAEqualSignRule :: Rule ComponentSet +addH0FromHAEqualSignRule = + describe "Rule for adding the H0 component based on HA; use equal sign (by convention)" . + makeRule "component.h0-from-ha-eq" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` NullHypothesis) + ha <- getRelation AlternativeHypothesis cs + return . append NullHypothesis + (CRelation $ h0FromHAEqualSign ha) $ cs + +addHARule :: Rule ComponentSet +addHARule = + describe "Rule for adding the HA component" . + makeRule "component.ha" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` AlternativeHypothesis) + ha <- get AlternativeHypothesis (initials cs) + return $ append AlternativeHypothesis ha cs + +---------------------------------------------------------- +-- Rules for determining the properties of the data + +chooseTTestRule :: Rule ComponentSet +chooseTTestRule = + describe "Rule for choosing the type of test" . + makeRule "component.test.t-test" $ + addTestChoice TTestOne + +chooseTTestTwoRule :: Rule ComponentSet +chooseTTestTwoRule = + describe "Rule for choosing the type of test" . + makeRule "component.test.t-test-two" $ + addTestChoice TTestTwo + +chooseTTestPairedRule :: Rule ComponentSet +chooseTTestPairedRule = + describe "Rule for choosing the type of test" . + makeRule "component.test.t-test-paired" $ + addTestChoice TTestPaired + +-- | If the standard deviation of the population is known then the z-test is +-- /always/ chosen, thus a t-test is never chosen in that case. When this +-- information is unknown, then an estimation of the standard deviation of the +-- population has to be made based on the sample, in that case a t-test is +-- chosen. +-- +-- When the sample size becomes large enough* it is also possible to choose a +-- z-test because the z-distribution looks like the t-distribution for large +-- sample sizes. +-- +-- NOTE*: The threshold for what `large' means may vary. For now this threshold +-- is fixed at 100. +-- +-- Consisely: If PopulationSdev known: ZTest +-- else if SampleSize big enough: ZTest or TTest +-- else TTest +chooseZTestRule :: Rule ComponentSet +chooseZTestRule = + describe "Rule for choosing the type of test" . + makeRule "component.test.z-test" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = + let largeThreshold = 100 + in do + guard (derived cs `doesNotContain` TestChoice) + n <- match naturalView $ fromMaybe 0 (getExpr SampleSize cs) + if cs `contains` PopulationSdev || n >= largeThreshold + then return (append TestChoice (CChoice $ TestType ZTest) cs) + else do + -- This branch is similar to the previous rule body in revision 10549 + let tests = validTests cs + guard (ZTest `elem` tests) + return $ append TestChoice (CChoice $ TestType ZTest) cs + +addTestFormulaRule :: Rule ComponentSet +addTestFormulaRule = + describe "Rule for adding the test formula" . + makeRule "component.test-formula" $ f + where + f :: ComponentSet -> [ComponentSet] + f cs = do + guard (derived cs `doesNotContain` TestFormula) + rel <- inferTestFormula cs + return $ append TestFormula (CRelation rel) cs + +addTestValueRule :: Rule ComponentSet +addTestValueRule = + describe "Rule for adding the test value (from the formula)" . + makeRule "component.test-value" $ f + where + f :: ComponentSet -> [ComponentSet] + f cs = do + guard (derived cs `doesNotContain` TestValue) + tv <- inferTestValue cs + return $ append TestValue (CRelation tv) cs + +inferTestValue :: MonadPlus m => ComponentSet -> m (Relation Expr) +inferTestValue cs = + case getRelation TestValue (initials cs) of + Just initialTestValue -> + return initialTestValue + _ | inferTestChoice cs == Just ChiSquared -> + chiSquaredTestValue cs + _ -> do + new <- msum (map return $ applyAll addTestFormulaRule cs) + let cs' = substitute new + var <- leftHandSide <$> getRelation TestFormula new + expr <- rightHandSide <$> getRelation TestFormula cs' + val <- matchM doubleView expr + return $ var .==. toExpr val + +---------------------------------------------------------- +-- Rules for performing a T-Test + +addDfRule :: Rule ComponentSet +addDfRule = + describe "Rule for adding the degrees of freedom component" . + makeRule "component.df" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` Df) + case get Df (initials cs) of + Just initialDf -> + return $ append Df initialDf cs + _ | inferTestChoice cs == Just ChiSquared -> do + df <- chiSquaredDf cs + return $ append Df (CExpr (toExpr df)) cs + _ -> do + test <- inferTestChoice cs + df <- degreesOfFreedomFromTest test (chooseTypeOfTest cs) + val <- matchM doubleView (getSubstitution cs |-> df) + return $ append Df (CExpr (toExpr val)) cs + +lookupTValueRule :: Rule ComponentSet +lookupTValueRule = + describe "Rule for looking up a t-value" . + makeRule "component.critical.t-value" $ + inferCriticalTWith $ \sided alpha df -> + [computeCriticalT sided alpha df] + +-- shared function for computing t-value and r-value +inferCriticalGenericWith :: (TestType -> Bool) -> Expr -> (Sided -> Double -> Double -> [Double]) + -> ComponentSet -> [ComponentSet] +inferCriticalGenericWith forTestType var compute cs = do + guard (cs `doesNotContain` Critical) + guard (derived cs `contains` AlternativeHypothesis) + guard (any forTestType (inferTestChoices cs)) + cs' <- matchM substitutedView cs + alpha <- matchM doubleView <=< getExpr SignificanceLevel $ cs + let cs'' = case inferDf cs' of + Just df -> substitute (append Df (CExpr df) cs') + Nothing -> cs' + df <- matchM doubleView =<< getExpr Df cs'' + sided <- inferSidedness cs + value <- compute sided alpha df + return $ append Critical (CRelation $ var .==. fromDouble value) cs + +inferCriticalTWith :: (Sided -> Double -> Double -> [Double]) -> ComponentSet -> [ComponentSet] +inferCriticalTWith = inferCriticalGenericWith isTTest (Var "tcrit") + +computePValueTTest :: Rule ComponentSet +computePValueTTest = + describe "Rule for computing the p-value for a t-test" . + makeRule "component.p-value.t-test" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` PValue) + guard (derived cs `contains` AlternativeHypothesis) + guard (maybe False isTTest (inferTestChoice cs)) + let cs' = substitute cs + testStatistic <- match doubleView =<< fmap rightHandSide (getRelation TestValue cs') + + let cs'' = case inferDf cs' of + Just df -> substitute (append Df (CExpr df) cs') + Nothing -> cs' + df <- matchM doubleView =<< getExpr Df cs'' + sided <- inferSidedness cs + let value = computePValueT sided testStatistic df + return $ append PValue (CExpr $ fromDouble value) cs + +---------------------------------------------------------- +-- Rules for performing a Z-Test + +lookupZValueRule :: Rule ComponentSet +lookupZValueRule = + describe "Rule for looking up a z-value" . + makeRule "component.critical.z-value" $ + inferCriticalZWith $ \sided alpha -> + [computeCriticalZ sided alpha] + +inferCriticalZWith :: (Sided -> Double -> [Double]) -> ComponentSet -> [ComponentSet] +inferCriticalZWith compute cs = do + guard (cs `doesNotContain` Critical) + guard (derived cs `contains` AlternativeHypothesis) + guard (maybe False (ZTest ==) (inferTestChoice cs)) + alpha <- matchM doubleView <=< getExpr SignificanceLevel $ cs + sided <- inferSidedness cs + value <- compute sided alpha + return $ append Critical (CRelation $ Var "zcrit" .==. fromDouble value) cs + +computePValueZTest :: Rule ComponentSet +computePValueZTest = + describe "Rule for computing the p-value for a z-test" . + makeRule "component.p-value.z-test" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` PValue) + guard (maybe False (ZTest ==) (inferTestChoice cs)) + guard (derived cs `contains` AlternativeHypothesis) + let cs' = substitute cs + -- TestStatistic was renamed to TestFormula, which is now a CRelation + testStatistic <- match doubleView =<< fmap rightHandSide (getRelation TestValue cs') + sided <- inferSidedness cs + let value = computePValueZ sided testStatistic + return $ append PValue (CExpr $ fromDouble value) cs + +---------------------------------------------------------- +-- Rules for performing R-Pearson + +lookupRValueRule :: Rule ComponentSet +lookupRValueRule = + describe "Rule for looking up a r-value" . + makeRule "component.critical.r-value" $ + inferCriticalRWith $ \sided alpha df -> + [computeCriticalR sided alpha df] + +chooseRPearsonRule :: Rule ComponentSet +chooseRPearsonRule = + describe "Rule for choosing the type of test" . + makeRule "component.test.r-pearson" $ + addTestChoice RPearson + +addTestChoice :: TestType -> ComponentSet -> Maybe ComponentSet +addTestChoice testType cs = do + guard (derived cs `doesNotContain` TestChoice) + let tests = validTests cs + guard (testType `elem` tests) + return $ append TestChoice (CChoice $ TestType testType) cs + +inferCriticalRWith :: (Sided -> Double -> Double -> [Double]) -> ComponentSet -> [ComponentSet] +inferCriticalRWith = inferCriticalGenericWith (== RPearson) (Var "rcrit") + +---------------------------------------------------------- +-- Rules for performing Anova + +lookupFValueRule :: Rule ComponentSet +lookupFValueRule = + describe "Rule for looking up a F-value" . + makeRule "component.critical.f-value" $ inferCriticalFWith computeCriticalF + +chooseAnovaRule :: Rule ComponentSet +chooseAnovaRule = + describe "Rule for choosing the type of test" . + makeRule "component.test.anova" $ + addTestChoice Anova + +addDfBetweenWithin :: Rule ComponentSet +addDfBetweenWithin = describe "Add df between and within (for Anova)" $ + makeRule "component.df-anova" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` DfBetween) + guard (derived cs `doesNotContain` DfWithin) + (between, within) <- inferDfBetweenWithin cs + return $ append DfBetween (CExpr $ toExpr between) $ + append DfWithin (CExpr $ toExpr within) cs + +inferDfBetweenWithin :: Monad m => ComponentSet -> m (Double, Double) +inferDfBetweenWithin cs = do + n <- matchM doubleView <=< getExpr SampleSize $ cs + let nrOfGroups = 2 -- always 2, for now + dfBetween = nrOfGroups - 1 + dfWithin = n - nrOfGroups + return (dfBetween, dfWithin) + +inferCriticalFWith :: (Double -> Double -> Double -> [Double]) -> ComponentSet -> [ComponentSet] +inferCriticalFWith compute cs = do + guard (cs `doesNotContain` Critical) + guard (derived cs `contains` AlternativeHypothesis) + guard (maybe False (Anova ==) (inferTestChoice cs)) + alpha <- matchM doubleView <=< getExpr SignificanceLevel $ cs + (dfBetween, dfWithin) <- inferDfBetweenWithin cs + value <- compute dfBetween dfWithin alpha + return $ append Critical (CRelation $ Var "Fcrit" .==. fromDouble value) cs + +---------------------------------------------------------- +-- Rules for performing Chi-Squared + +chooseChiSquaredRule :: Rule ComponentSet +chooseChiSquaredRule = describe "Rule for choosing the type of test" . + makeRule "component.test.chi-squared" $ + addTestChoice ChiSquared + +addHypothesesChiSquaredRule :: Rule ComponentSet +addHypothesesChiSquaredRule = describe "Add hypotheses (null and alternative) for chi-squared" $ + makeRule "component.hypotheses-chi-squared" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` NullHypothesis) + guard (derived cs `doesNotContain` AlternativeHypothesis) + guard (maybe False (ChiSquared ==) (inferTestChoice cs)) + return $ append NullHypothesis (CExpr independent) + $ append AlternativeHypothesis (CExpr dependent) cs + +addObservedTotals :: Rule ComponentSet +addObservedTotals = describe "Add totals (rows and columns) for observed frequencies" $ + makeRule "component.observed-totals" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` ObservedColumnTotals) + guard (derived cs `doesNotContain` ObservedRowTotals) + guard (derived cs `doesNotContain` ObservedTotal) + table <- getTable ObservedFrequencies cs + let (rowTotals, columnTotals, total) = computeTotals table + return $ append ObservedRowTotals (CExpr $ toExpr rowTotals) + $ append ObservedColumnTotals (CExpr $ toExpr columnTotals) + $ append ObservedTotal (CExpr $ toExpr total) cs + +type ChiSquaredTotals = ([Int], [Int], Int) + +computeTotals :: [[Int]] -> ChiSquaredTotals +computeTotals table = + let rowTotals = map sum table + columnTotals = map sum (transpose table) + total = sum rowTotals + in (rowTotals, columnTotals, total) + +computeExpectedFrequencies :: ChiSquaredTotals -> [[Double]] +computeExpectedFrequencies (rowTotals, columnTotals, total) = table + where + xss = map (replicate (length columnTotals)) rowTotals + yss = replicate (length rowTotals) columnTotals + table = zipWith (zipWith f) xss yss + f x y = fromIntegral (x*y) / fromIntegral total + +addExpectedFrequencies :: Rule ComponentSet +addExpectedFrequencies = describe "Add expected frequencies" $ + makeRule "component.expected-frequencies" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (derived cs `doesNotContain` ExpectedFrequencies) + observed <- getTable ObservedFrequencies cs + let totals = computeTotals observed + expected = computeExpectedFrequencies totals + return $ append ExpectedFrequencies (CExpr $ toExpr expected) cs + +lookupChiValueRule :: Rule ComponentSet +lookupChiValueRule = describe "Rule for looking up a chi^2-value" . + makeRule "component.critical.chi-value" $ + inferCriticalChiWith $ \sided alpha df -> + computeCriticalChi sided alpha df + +inferCriticalChiWith :: (Sided -> Double -> Double -> [Double]) -> ComponentSet -> [ComponentSet] +inferCriticalChiWith = inferCriticalGenericWith (== ChiSquared) (Var "chicrit") + +chiSquaredTestValue :: MonadPlus m => ComponentSet -> m (Relation Expr) +chiSquaredTestValue cs = do + observed <- getTable ObservedFrequencies cs + expected <- getDoubleTable ExpectedFrequencies cs + let table = zipWith (zipWith f) observed expected + f o e = (fromIntegral o-e)^(2 :: Int) / e + value = sum (map sum table) + return (chiSquared .==. toExpr value) + +chiSquaredDf :: MonadPlus m => ComponentSet -> m Int +chiSquaredDf cs = do + observed <- getTable ObservedFrequencies cs + guard (not $ null observed) + let r = length observed + c = length $ head observed + return ((r-1)*(c-1)) + +getTable :: MonadPlus m => ComponentId -> ComponentSet -> m [[Int]] +getTable n cs = getExpr n cs >>= fromExpr + +getDoubleTable :: MonadPlus m => ComponentId -> ComponentSet -> m [[Double]] +getDoubleTable n cs = getExpr n cs >>= fromExpr + +---------------------------------------------------------- +-- Rules for constructing the rejection region + +-- lookup the Sidedness component; if not present, try to infer (using the rule) +inferSidedness :: MonadPlus m => ComponentSet -> m Sided +inferSidedness cs = + getSided Sidedness cs `mplus` do + testType <- inferTestChoice cs + guard (testType `elem` [Anova, ChiSquared]) + return RightSided + `mplus` do + cs' <- applyM determineSided cs + getSided Sidedness cs' + +inferVar :: MonadPlus m => ComponentSet -> m Expr +inferVar cs = do + testType <- inferTestChoice cs + return (varForTestType testType) + +varForTestType :: TestType -> Expr +varForTestType testType = + case testType of + ZTest -> Var "z" + RPearson -> Var "r" + Anova -> Var "F" + ChiSquared -> chiSquared + _ -> Var "t" + +inferTestChoice :: MonadPlus m => ComponentSet -> m TestType +inferTestChoice cs = + case validTests cs of + hd:_ -> return hd + _ -> fail "no valid test choice" + +inferTestChoices :: ComponentSet -> [TestType] +inferTestChoices cs = + case getTestType TestChoice cs of + Just test -> [test] + Nothing -> validTests cs + +inferTestFormula :: ComponentSet -> [Relation Expr] +inferTestFormula cs = + case getRelation TestFormula cs of + Just formula -> [formula] + Nothing -> do + testType <- inferTestChoices cs + let var = varForTestType testType + let test = chooseTypeOfTest cs + t <- testFormulaFromTest testType test + return (var .==. t) + +inferRejectionCritical :: MonadPlus m => ComponentSet -> m (Relation Expr) +inferRejectionCritical cs = + getRelation RejectionCritical cs `mplus` do + cs' <- applyM addRejectionRule cs + getRelation RejectionCritical cs' + +inferConclusionCritical :: MonadPlus m => ComponentSet -> m Bool +inferConclusionCritical cs = + getConclusion ConclusionCritical cs `mplus` do + cs' <- applyM criticalConclusionRule cs + getConclusion ConclusionCritical cs' + +inferConclusionPValue :: MonadPlus m => ComponentSet -> m (Relation Expr) +inferConclusionPValue cs = + getRelation ConclusionPValue cs `mplus` do + cs' <- applyM addConclusionPValueRule cs + getRelation ConclusionPValue cs' + +inferDf :: MonadPlus m => ComponentSet -> m Expr +inferDf cs = + getExpr Df cs `mplus` do + cs' <- applyM addDfRule cs + getExpr Df cs' + +addRejectionRule :: Rule ComponentSet +addRejectionRule = + describe "Rule for constructing the rejection critical component" . + makeRule "component.rejection.critical" $ f + where + f :: ComponentSet -> [ComponentSet] + f cs = do + guard (cs `doesNotContain` RejectionCritical) + guard (derived cs `contains` AlternativeHypothesis) + let cs' = substitute cs + sided <- inferSidedness cs' + testType <- inferTestChoices cs + let rel = case testType of + ZTest -> sidedRelation sided (Var "z") (Var "zcrit") + RPearson -> sidedRelation sided (Var "r") (Var "rcrit") + Anova -> sidedRelation sided (Var "F") (Var "Fcrit") + ChiSquared -> sidedRelation sided chiSquared (Var "chicrit") + _ -> sidedRelation sided (Var "t") (Var "tcrit") + return . append RejectionCritical + (CRelation rel) $ cs + +addConclusionPValueRule :: Rule ComponentSet +addConclusionPValueRule = + describe "Rule for constructing the conclusion p-value component" . + makeRule "component.conclusion.p-value" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` ConclusionPValue) + pv <- matchM doubleView =<< getExpr PValue cs + alpha <- matchM doubleView =<< getExpr SignificanceLevel cs + let relType = if pv > alpha then GreaterThan else LessThanOrEqualTo + return . append ConclusionPValue + (CRelation $ makeType relType (Var "p") (Var "alpha")) $ cs + +---------------------------------------------------------- +-- Rules for making a conclusion + +criticalConclusionRule :: Rule ComponentSet +criticalConclusionRule = + describe ("Rule for creating a conclusion based on a critical value " ++ + "and test statistic") . + makeRule "component.critical-conclusion" $ f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + rej <- inferRejectionCritical cs + -- to do: rejection critical is added to component set only to get the substituted relation + let cs' = substitute (append RejectionCritical (CRelation rej) cs) + rejection <- getRelation RejectionCritical cs' + guard (cs `doesNotContain` ConclusionCritical) + lhs <- match doubleView $ leftHandSide rejection + rhs <- match doubleView $ rightHandSide rejection + let result = eval (relationType rejection) lhs rhs + return $ append ConclusionCritical (CChoice $ Conclusion result) cs + +hypothesesConclusionCriticalRule :: Rule ComponentSet +hypothesesConclusionCriticalRule = + describe "derive the hypotheses conclusion from the critical conclusion" $ + makeRule "component.hypotheses-conclusion-critical" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` ConclusionHypotheses) + concl <- inferConclusionCritical cs + let rejhyp = if concl then RejectH0 else DontRejectH0 + return $ append ConclusionHypotheses (CChoice (RejectionHypotheses rejhyp)) cs + +hypothesesConclusionPValueRule :: Rule ComponentSet +hypothesesConclusionPValueRule = + describe "derive the hypotheses conclusion from the p-value" $ + makeRule "component.hypotheses-conclusion-pvalue" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` ConclusionHypotheses) + rel <- inferConclusionPValue cs + let rejhyp = if relationType rel == LessThanOrEqualTo then RejectH0 else DontRejectH0 + return $ append ConclusionHypotheses (CChoice (RejectionHypotheses rejhyp)) cs + +addStandardErrorSigma :: Rule ComponentSet +addStandardErrorSigma = + describe "derive standard error from population standard deviation and sample size" $ + makeRule "component.standard-error-sigma" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` StandardError) + n <- matchM doubleView =<< getExpr SampleSize cs + psdev <- matchM doubleView =<< getRhsExpr PopulationSdev cs + let se = Var "sigmaM" .==. toExpr (psdev / sqrt n) + return $ append StandardError (CRelation se) cs + +addStandardErrorSD :: Rule ComponentSet +addStandardErrorSD = + describe "derive standard error from sample standard deviation and sample size" $ + makeRule "component.standard-error-sd" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do + guard (cs `doesNotContain` StandardError) + n <- matchM doubleView =<< getExpr SampleSize cs + sdev <- matchM doubleView =<< getRhsExpr SampleSdev cs + let se = Var "SEM" .==. toExpr (sdev / sqrt n) + return $ append StandardError (CRelation se) cs diff --git a/src/Domain/Hypothesis/Strategies.hs b/src/Domain/Hypothesis/Strategies.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Hypothesis/Strategies.hs @@ -0,0 +1,74 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Hypothesis.Strategies (hypothesisStrategy) where + +import Domain.Hypothesis.Rules +import Domain.Statistics.Rules +import Domain.Statistics.ComponentSet +import Ideas.Common.Strategy.Combinators hiding (not, while) +import Ideas.Common.Library hiding (while) +import Prelude hiding (until, sequence, repeat) + + +---------------------------------------------------------- +-- Template for the hypothesis testing strategy + +hypothesisStrategy :: LabeledStrategy ComponentSet +hypothesisStrategy = label "Hypothesis testing" $ + label "Preparation" (whileNotReady $ choice $ + [ addHypothesesRule, addH0FromHARule, addH0FromHAEqualSignRule, addHARule + , addHypothesesChiSquaredRule + , addAlphaRule, determineSided, chooseTTestRule + , chooseTTestTwoRule, chooseTTestPairedRule, chooseZTestRule + , chooseRPearsonRule, chooseAnovaRule, chooseChiSquaredRule + ] ++ sampleStatistics) + .*. + check (\cs -> all (derived cs `contains`) [NullHypothesis, AlternativeHypothesis]) + .*. + label "Computation" (whileNotReady $ + (check (\cs -> derived cs `doesNotContain` TestValue) .*. + (addTestFormulaRule .|. choice sampleStatistics)) + .|. (check allowCriticalRoute .*. choice + [ addTestValueRule, addRejectionRule + , lookupZValueRule, lookupTValueRule, lookupRValueRule, lookupFValueRule, lookupChiValueRule + ]) + .|. (check allowPValueRoute .*. choice + [ computePValueZTest, computePValueTTest + ]) + ) + .*. + check (\cs -> derived cs `contains` TestValue && + derived cs `contains` Critical || derived cs `contains` PValue) + .*. + label "Conclusion" ( + whileNotReady (criticalConclusionRule .|. addConclusionPValueRule) + .*. + (hypothesesConclusionCriticalRule .|. hypothesesConclusionPValueRule)) + where + sampleStatistics = + [ addNRule, addAverageRule, addVarianceRule, addStandardDeviationRule + , addStandardErrorRule, addStandardErrorSD, addDfRule, addDfBetweenWithin + , addStandardErrorSigma, addObservedTotals, addExpectedFrequencies + ] + +-- customized while (combinator in library is greedy and uses repeat instead of many) +while :: IsStrategy f => (a -> Bool) -> f a -> Strategy a +while p s = many (check p .*. s) + +whileNotReady :: Strategy ComponentSet -> Strategy ComponentSet +whileNotReady = while (`doesNotContain` ConclusionHypotheses) + +allowCriticalRoute :: ComponentSet -> Bool +allowCriticalRoute cs = criticalRoute cs || not (pvalueRoute cs) + +allowPValueRoute :: ComponentSet -> Bool +allowPValueRoute cs = pvalueRoute cs || not (criticalRoute cs) + +criticalRoute :: ComponentSet -> Bool +criticalRoute cs = derived cs `contains` RejectionCritical || derived cs `contains` Critical + +pvalueRoute :: ComponentSet -> Bool +pvalueRoute cs = derived cs `contains` PValue diff --git a/src/Domain/Hypothesis/Tables.hs b/src/Domain/Hypothesis/Tables.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Hypothesis/Tables.hs @@ -0,0 +1,88 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Hypothesis.Tables + ( zTableS, zTable, + tTableS, tTable, + fTable, chiTable + ) where + +zTableS :: Bool -> Double -> Maybe Double +zTableS True alpha = zTable alpha +zTableS False alpha = zTable (alpha * 2) + +zTable :: Double -> Maybe Double +zTable 0.05 = Just 1.645 +zTable 0.025 = Just 1.960 +zTable 0.01 = Just 2.326 +zTable 0.005 = Just 2.576 +zTable 0.0025 = Just 2.807 +zTable 0.001 = Just 3.090 +zTable _ = Nothing + +tTableS :: Bool -> Double -> Int -> Maybe Double +tTableS True alpha df = tTable alpha df +tTableS False alpha df = tTable (alpha * 2) df + +tTable :: Double -> Int -> Maybe Double +tTable alpha df | df > 60 = lookup alpha [(0.1, 1.282), (0.05, 1.645), (0.025, 1.960), (0.01, 2.326), (0.005, 2.576)] + | otherwise = tTable' alpha df + +tTable' :: Double -> Int -> Maybe Double +tTable' 0.1 df = lookup df [(1, 3.078), (2, 1.886), (3, 1.638), (4, 1.533), (5, 1.476), (6, 1.440), (7, 1.415), (8, 1.397), (9, 1.383), (10, 1.372), + (11, 1.363), (12, 1.356), (13, 1.350), (14, 1.345), (15, 1.341), (16, 1.337), (17, 1.333), (18, 1.330), (19, 1.328), (20, 1.325), + (21, 1.323), (22, 1.321), (23, 1.319), (24, 1.318), (25, 1.316), (26, 1.315), (27, 1.314), (28, 1.313), (29, 1.311), (30, 1.310), + (40, 1.303), (60, 1.296)] +tTable' 0.05 df = lookup df [(1, 6.314), (2, 2.920), (3, 2.353), (4, 2.132), (5, 2.015), (6, 1.943), (7, 1.895), (8, 1.860), (9, 1.833), (10, 1.812), + (11, 1.796), (12, 1.782), (13, 1.771), (14, 1.761), (15, 1.753), (16, 1.746), (17, 1.740), (18, 1.734), (19, 1.729), (20, 1.725), + (21, 1.721), (22, 1.717), (23, 1.714), (24, 1.711), (25, 1.708), (26, 1.706), (27, 1.703), (28, 1.701), (29, 1.699), (30, 1.697), + (40, 1.684), (60, 1.671)] +tTable' 0.025 df = lookup df [(1, 12.706), (2, 4.303), (3, 3.182), (4, 2.776), (5, 2.571), (6, 2.447), (7, 2.365), (8, 2.306), (9, 2.262), (10, 2.228), + (11, 2.201), (12, 2.179), (13, 2.160), (14, 2.145), (15, 2.131), (16, 2.120), (17, 2.110), (18, 2.101), (19, 2.093), (20, 2.086), + (21, 2.080), (22, 2.074), (23, 2.069), (24, 2.064), (25, 2.060), (26, 2.056), (27, 2.052), (28, 2.048), (29, 2.045), (30, 2.042), + (40, 2.021), (60, 2.000)] +tTable' 0.01 df = lookup df [(1, 31.821), (2, 6.965), (3, 4.541), (4, 3.747), (5, 3.365), (6, 3.143), (7, 2.998), (8, 2.896), (9, 2.821), (10, 2.764), + (11, 2.718), (12, 2.681), (13, 2.650), (14, 2.624), (15, 2.602), (16, 2.583), (17, 2.567), (18, 2.552), (19, 2.539), (20, 2.528), + (21, 2.518), (22, 2.508), (23, 2.500), (24, 2.492), (25, 2.485), (26, 2.479), (27, 2.473), (28, 2.467), (29, 2.462), (30, 2.457), + (40, 2.423), (60, 2.390)] +tTable' 0.005 df = lookup df [(1, 63.657), (2, 9.925), (3, 5.841), (4, 4.604), (5, 4.032), (6, 3.707), (7, 3.499), (8, 3.355), (9, 3.250), (10, 3.169), + (11, 3.106), (12, 3.055), (13, 3.012), (14, 2.997), (15, 2.947), (16, 2.921), (17, 2.898), (18, 2.878), (19, 2.861), (20, 2.845), + (21, 2.831), (22, 2.819), (23, 2.807), (24, 2.797), (25, 2.787), (26, 2.779), (27, 2.771), (28, 2.763), (29, 2.756), (30, 2.750), + (40, 2.704), (60, 2.660)] +tTable' _ _ = Nothing + +-- for Anova +fTable :: Double -> Double -> Double -> Maybe Double +fTable dfBetween dfWithin alpha + | alpha == 0.05 = lookup (dfBetween, dfWithin) + [((1, 38), 4.10), ((1, 40), 4.08), ((2, 38), 3.25), ((2, 40), 3.23) ] + | alpha == 0.01 = lookup (dfBetween, dfWithin) + [((1, 38), 7.35), ((1, 40), 7.31), ((2, 38), 5.21), ((2, 40), 5.18) ] + | otherwise = Nothing + +-- for Chi squared +chiTable :: Double -> Int -> Maybe Double +chiTable alpha df + | alpha == 0.10 && df > 0 && df <= 20 = Just $ + [ 2.71, 4.61, 6.25, 7.78, 9.24, 10.64, 12.02, 13.36, 14.68, 15.99 + , 17.28, 18.55, 19.81, 21.06, 22.31, 23.54, 24.77, 25.99, 27.20, 28.41 + ] !! (df - 1) + | alpha == 0.05 && df > 0 && df <= 20 = Just $ + [ 3.84, 5.99, 7.81, 9.49, 11.07, 12.59, 14.07, 15.51, 16.92, 18.31 + , 19.68, 21.03, 22.36, 23.68, 25.00, 26.30, 27.59, 28.87, 30.14, 31.41 + ] !! (df - 1) + | alpha == 0.025 && df > 0 && df <= 20 = Just $ + [ 5.02, 7.38, 9.35, 11.14, 12.83, 14.45, 16.01, 17.53, 19.02, 20.48 + , 21.92, 23.34, 24.74, 26.12, 27.49, 28.85, 30.19, 31.53, 32.85, 34.17 + ] !! (df - 1) + | alpha == 0.01 && df > 0 && df <= 20 = Just $ + [ 6.63, 9.21, 11.34, 13.28, 15.09, 16.81, 18.48, 20.09, 21.67, 23.21 + , 24.72, 26.22, 27.69, 29.14, 30.58, 32.00, 33.41, 34.81, 36.19, 37.57 + ] !! (df - 1) + | alpha == 0.005 && df > 0 && df <= 20 = Just $ + [ 7.88, 10.60, 12.84, 14.86, 16.75, 18.55, 20.28, 21.96, 23.59, 25.19 + , 26.76, 28.30, 29.82, 31.32, 32.80, 34.27, 35.72, 37.16, 38.58, 40.00 + ] !! (df - 1) + | otherwise = Nothing diff --git a/src/Domain/Statistics/Component.hs b/src/Domain/Statistics/Component.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Statistics/Component.hs @@ -0,0 +1,328 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +{-# LANGUAGE FlexibleInstances #-} + +module Domain.Statistics.Component + ( Choice (..), Component (..), ComponentId (..) + , Sided (..), TestType (..), RejectionHypotheses(..) + , choices, getCId, fromCId, sidedRelation + , isData, isExpr, isRelation + , isTestType, isSided, isRejectionHypotheses, isConclusion + ) where + +import Control.Applicative ((<|>)) +import Control.Monad +import Data.Char +import Data.Function +import Data.List +import Domain.Statistics.Symbols +import Domain.Math.Data.Relation +import Domain.Math.Expr +import Domain.Statistics.Data +import Ideas.Common.Rewriting hiding (trueSymbol,falseSymbol) +import Prelude hiding (lookup) + +data Component + = CData Data -- Data component type (sample data) + | CExpr Expr -- Expression component type + | CRelation (Relation Expr) + | CChoice Choice -- Expressing various choices + deriving Eq + +instance Show Component where + show (CRelation r) = show r + show (CExpr e) = show e + show (CChoice c) = show c + show (CData d) = show d + +instance IsTerm Component where + toTerm (CData x) = unary cDataSymbol (toTerm x) + toTerm (CChoice x) = toTerm x + toTerm (CRelation x) = toTerm x + toTerm (CExpr x) = toTerm x + + -- Check to see if the whole x can be parsed as an expression + fromTerm (TCon s [val]) + | s == cDataSymbol = CData <$> fromTerm val + -- First try to convert it to a relation because relation is more specific + -- than an Expression (IMPORTANT) + fromTerm t = CChoice <$> fromTerm t <|> + CRelation <$> fromTerm t <|> + CExpr <$> fromTerm t + +data Choice + = TestType TestType + | Sided Sided + | RejectionHypotheses RejectionHypotheses + | Conclusion Bool + deriving Eq + +instance Show Choice where + show (TestType tt) = show tt + show (Sided sd) = show sd + show (RejectionHypotheses rh) = show rh + show (Conclusion c) = map toLower (show c) + +instance IsTerm Choice where + toTerm (TestType x) = toTerm x + toTerm (Sided x) = toTerm x + toTerm (RejectionHypotheses rh) = toTerm rh + toTerm (Conclusion x) = toTermBool x + + fromTerm t = + Sided <$> fromTerm t <|> + TestType <$> fromTerm t <|> + RejectionHypotheses <$> fromTerm t <|> + Conclusion <$> fromTermBool t + +choices :: [Choice] +choices = [ TestType tt | tt <- [minBound .. maxBound] ] ++ + [ Sided sd | sd <- [minBound .. maxBound] ] ++ + [ RejectionHypotheses rj | rj <- [minBound .. maxBound] ] ++ + [ Conclusion c | c <- [minBound .. maxBound] ] + +data TestType + = TTestOne + | TTestTwo + | TTestPaired + | ZTest + | RPearson + | Anova + | ChiSquared + deriving (Eq, Enum, Bounded) + +instance Show TestType where + -- to do: use same string as corresponding symbols (without dashes) + show TTestOne = "t-test-one-sample" + show TTestTwo = "t-test-two-sample" + show TTestPaired = "t-test-paired" + show ZTest = "z-test" + show RPearson = "r-pearson" + show Anova = "anova" + show ChiSquared = "chi-squared" + +data Sided + = TwoSided + | LeftSided + | RightSided + deriving (Eq, Enum, Bounded) + +instance Show Sided where + show TwoSided = "two" + show LeftSided = "left" + show RightSided = "right" + +sidedRelation :: Sided -> Expr -> Expr -> Relation Expr +sidedRelation LeftSided lhs rhs = makeType LessThan lhs rhs +sidedRelation RightSided lhs rhs = makeType GreaterThan lhs rhs +sidedRelation TwoSided lhs rhs = makeType GreaterThan (abs lhs) rhs + +data RejectionHypotheses + = RejectH0 + | AcceptH0 + | DontRejectH0 + | RejectH1 + | AcceptH1 + | DontRejectH1 + deriving (Eq, Enum, Bounded) + +instance Show RejectionHypotheses where + show RejectH0 = "rejecth0" + show AcceptH0 = "accepth0" + show DontRejectH0 = "dontrejecth0" + show RejectH1 = "rejecth1" + show AcceptH1 = "accepth1" + show DontRejectH1 = "dontrejecth1" + +data ComponentId + = DataSet + | SampleSize -- n + | SampleMean -- M + | SampleVariance -- s^2 + | SampleSdev -- s + | ObservedFrequencies -- fo (two-dimensional matrix), for chi squared + | ObservedRowTotals + | ObservedColumnTotals + | ObservedTotal + | ExpectedFrequencies -- fe (two-dimensional matrix), for chi squared + | PopulationMean -- mu + | PopulationSdev -- sigma + | StandardError -- sigmaM or SEM + | SignificanceLevel -- alpha + | NullHypothesis + | AlternativeHypothesis + | PooledVariance -- S_p + | Df -- df + | DfBetween -- for Anova + | DfWithin -- for Anova + | Groups + | Sidedness + | TestChoice + | Critical + | PValue + | RejectionCritical + | ConclusionPValue + | ConclusionCritical + | ConclusionHypotheses + | TestFormula + | TestValue + | Proportion -- P + | Correlation -- r + | One ComponentId -- Belonging to the first sample (in case of multiple samples) + | Two ComponentId -- Belonging to the second sample + | Other String + deriving (Eq, Show) + +instance Ord ComponentId where + compare = compare `on` getCId + +instance IsTerm ComponentId where + toTerm = toTerm . getCId + fromTerm t = fromCId <$> fromTerm t + +getCId :: ComponentId -> String +getCId (Other nm) = nm +getCId (One nm) = getCId nm ++ "1" +getCId (Two nm) = getCId nm ++ "2" +getCId cId = case lookup cId cIdTable of + Just nm -> nm + _ -> error "Missing name for ComponentId" +fromCId :: String -> ComponentId +fromCId s + | suffix == "1" = One (fromCId $ init s) + | suffix == "2" = Two (fromCId $ init s) + | otherwise = + case [ x | (x, y) <- cIdTable, y == s ] of + [x] -> x + _ -> Other s + where + suffix = take 1 (reverse s) + +cIdTable :: [(ComponentId, String)] +cIdTable = + [ (DataSet, "data") + , (SampleSize, "n") + , (SampleMean, "samplemean") -- was: "mean" + , (SampleVariance, "samplevariance") -- was: "variance" + , (SampleSdev, "samplesd") -- was: "sd" + , (ObservedFrequencies, "fo") + , (ObservedRowTotals, "forowtotal") + , (ObservedColumnTotals, "focolumntotal") + , (ObservedTotal, "fototal") + , (ExpectedFrequencies, "fe") + , (PopulationMean, "populationmean") -- was: "pmean" + , (PopulationSdev, "populationsd") -- was "psd" + , (StandardError, "standarderror") -- was: "se" + , (SignificanceLevel, "alpha") + , (NullHypothesis, "h0") + , (AlternativeHypothesis, "ha") + , (PooledVariance, "pooledvariance") -- was "sp2" + , (Df, "df") + , (DfWithin, "dfwithin") + , (DfBetween, "dfbetween") + , (Groups, "groups") + , (Sidedness, "sided") + , (TestChoice, "test") + , (Critical, "critical") + , (PValue, "pvalue") + , (RejectionCritical, "rejectioncritical") -- was: "rejection" + , (ConclusionPValue, "conclusionpvalue") + , (ConclusionCritical, "conclusioncritical") -- was: "conclusion" + , (ConclusionHypotheses, "conclusionhypotheses") + , (TestFormula, "testformula") + , (TestValue, "testvalue") + , (Proportion, "p") + , (Correlation, "r") + ] + +---------------------------------------------------------- +-- Component utilities + +isData :: Monad m => Component -> m Data +isData (CData d) = return d +isData _ = failMsg "data" + +isExpr :: Monad m => Component -> m Expr +isExpr (CExpr e) = return e +isExpr _ = failMsg "expr" + +isRelation :: Monad m => Component -> m (Relation Expr) +isRelation (CRelation r) = return r +isRelation _ = failMsg "relation" + +isTestType :: Monad m => Component -> m TestType +isTestType (CChoice (TestType tt)) = return tt +isTestType _ = failMsg "testtype" + +isSided :: Monad m => Component -> m Sided +isSided (CChoice (Sided s)) = return s +isSided _ = failMsg "sided" + +isRejectionHypotheses :: Monad m => Component -> m RejectionHypotheses +isRejectionHypotheses (CChoice (RejectionHypotheses h)) = return h +isRejectionHypotheses _ = failMsg "rejectionhypotheses" + +isConclusion :: Monad m => Component -> m Bool +isConclusion (CChoice (Conclusion b)) = return b +isConclusion _ = failMsg "conclusion" + +failMsg :: Monad m => String -> m a +failMsg s = fail ("component is not of type " ++ s) + +---------------------------------------------------------- +-- Component utilities + +-- | Types that have only atomic constructors are easy and follow a pattern for +-- their implementation +instance IsTerm TestType where + toTerm TTestOne = symbol tTestOneSymbol + toTerm TTestTwo = symbol tTestTwoSymbol + toTerm TTestPaired = symbol tTestPairedSymbol + toTerm ZTest = symbol zTestSymbol + toTerm RPearson = symbol rPearsonSymbol + toTerm Anova = symbol anovaSymbol + toTerm ChiSquared = symbol chiSquaredSymbol + + fromTerm (TCon s []) + | s == tTestOneSymbol = return TTestOne + | s == tTestTwoSymbol = return TTestTwo + | s == tTestPairedSymbol = return TTestPaired + | s == zTestSymbol = return ZTest + | s == rPearsonSymbol = return RPearson + | s == anovaSymbol = return Anova + | s == chiSquaredSymbol = return ChiSquared + fromTerm (TVar "rpearson") = return RPearson -- temporary fix!!! + fromTerm (TVar "anova") = return Anova -- temporary fix!!! + fromTerm (TVar "chisquared") = return ChiSquared -- temporary fix!!! + fromTerm t = fail $ "Invalid term: " ++ show t ++ " not TestType" + +instance IsTerm Sided where + toTerm TwoSided = symbol twoSidedSymbol + toTerm LeftSided = symbol leftSidedSymbol + toTerm RightSided = symbol rightSidedSymbol + + fromTerm (TCon s []) + | s == twoSidedSymbol = return TwoSided + | s == leftSidedSymbol = return LeftSided + | s == rightSidedSymbol = return RightSided + fromTerm t = fail $ "Invalid term: " ++ show t ++ " not Sided" + +instance IsTerm RejectionHypotheses where + toTerm RejectH0 = symbol rejectH0Symbol + toTerm AcceptH0 = symbol acceptH0Symbol + toTerm DontRejectH0 = symbol dontRejectH0Symbol + toTerm RejectH1 = symbol rejectH1Symbol + toTerm AcceptH1 = symbol acceptH1Symbol + toTerm DontRejectH1 = symbol dontRejectH1Symbol + + fromTerm (TCon s []) + | s == rejectH0Symbol = return RejectH0 + | s == acceptH0Symbol = return AcceptH0 + | s == dontRejectH0Symbol = return DontRejectH0 + | s == rejectH1Symbol = return RejectH1 + | s == acceptH1Symbol = return AcceptH1 + | s == dontRejectH1Symbol = return DontRejectH1 + fromTerm t = fail $ "Invalid term: " ++ show t ++ " not RejectionHypotheses" diff --git a/src/Domain/Statistics/ComponentSet.hs b/src/Domain/Statistics/ComponentSet.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Statistics/ComponentSet.hs @@ -0,0 +1,163 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Statistics.ComponentSet + ( -- * Types + ComponentSet, module Domain.Statistics.Component + -- * Constructors + , initialSet, derivedSet + -- * Combining and transforming + , append, delete, mapComponent, initials, derived, toList + -- * Getters + , get, getData, getExpr, getRelation, getRhsExpr + , getTestType, getSided, getRejectionHypotheses, getConclusion + -- * Membership + , contains, doesNotContain + , isInitial, isDerived + ) where + +import Control.Monad +import Data.List hiding (delete) +import Data.Maybe +import Domain.Statistics.Symbols +import Domain.Math.Data.Relation +import Domain.Math.Expr +import Domain.Statistics.Data +import Domain.Statistics.Component +import Ideas.Common.Rewriting hiding (trueSymbol,falseSymbol) +import Ideas.Common.Classes + +------------------------------------------------------------------------------ +-- Types + +data ComponentSet = CS + { initialsList :: [(ComponentId, Component)] + , derivedList :: [(ComponentId, Component)] + } + +instance Show ComponentSet where + show cs = unlines (map (f True) (initialsList cs) ++ map (f False) (derivedList cs)) + where + f b (k, v) = bracketsIf b (getCId k) ++ ": " ++ show v ++ ";" + + bracketsIf True s = "[" ++ s ++ "]" + bracketsIf False s = s + +instance Eq ComponentSet where + CS xs1 xs2 == CS ys1 ys2 = f xs1 == f ys1 && f xs2 == f ys2 + where + f = sortBy (\(x, _) (y, _) -> x `compare` y) + +instance Semigroup ComponentSet where + CS xs1 xs2 <> CS ys1 ys2 = CS (xs1 ++ ys1) (xs2 ++ ys2) + +instance Monoid ComponentSet where + mempty = CS [] [] + mappend = (<>) + +instance IsTerm ComponentSet where + toTerm cs = TList $ + [ ternary componentSymbol (toTerm k) (symbol initialSymbol) (toTerm c) + | (k, c) <- initialsList cs + ] ++ + [ ternary componentSymbol (toTerm k) (symbol derivedSymbol) (toTerm c) + | (k, c) <- derivedList cs + ] + + fromTerm (TList xs) = + let f (TCon s [t1, TCon sv [], t3]) | s == componentSymbol = do + k <- fromTerm t1 + c <- fromTerm t3 + if sv == initialSymbol + then return $ CS [(k, c)] [] + else if sv == derivedSymbol + then return $ CS [] [(k, c)] + else fail "expected initial/derived" + f _ = fail "Term is not a component" + in mconcat <$> mapM f xs + fromTerm _ = fail "Term is not a component (not a list)" + +------------------------------------------------------------------------------ +-- Constructors + +initialSet :: [(ComponentId, Component)] -> ComponentSet +initialSet xs = CS xs [] + +derivedSet :: [(ComponentId, Component)] -> ComponentSet +derivedSet = CS [] + +------------------------------------------------------------------------------ +-- Combining and transforming + +append :: ComponentId -> Component -> ComponentSet -> ComponentSet +append k v cs = cs <> derivedSet [(k, v)] + +delete :: ComponentId -> ComponentSet -> ComponentSet +delete k (CS xs ys) = CS (f xs) (f ys) + where + f = filter ((/= k) . fst) + +mapComponent :: (Component -> Component) -> ComponentSet -> ComponentSet +mapComponent f (CS xs ys) = CS (change xs) (change ys) + where + change = map (mapSecond f) + +initials :: ComponentSet -> ComponentSet +initials cs = mempty { initialsList = initialsList cs } + +derived :: ComponentSet -> ComponentSet +derived cs = mempty { derivedList = derivedList cs } + +toList :: ComponentSet -> [(ComponentId, Component)] +toList cs = initialsList cs ++ derivedList cs + +------------------------------------------------------------------------------ +-- Getters + +get :: Monad m => ComponentId -> ComponentSet -> m Component +get n cs = + -- first consider derived components + case lookup n (derivedList cs ++ initialsList cs) of + Just c -> return c + Nothing -> fail $ show n ++ " missing" + +getData :: Monad m => ComponentId -> ComponentSet -> m Data +getData n = get n >=> isData + +getExpr :: Monad m => ComponentId -> ComponentSet -> m Expr +getExpr n = get n >=> isExpr + +getRelation :: Monad m => ComponentId -> ComponentSet -> m (Relation Expr) +getRelation n = get n >=> isRelation + +getRhsExpr :: MonadPlus m => ComponentId -> ComponentSet -> m Expr +getRhsExpr n cs = fmap rightHandSide (getRelation n cs) `mplus` getExpr n cs + +getTestType :: Monad m => ComponentId -> ComponentSet -> m TestType +getTestType n = get n >=> isTestType + +getSided :: Monad m => ComponentId -> ComponentSet -> m Sided +getSided n = get n >=> isSided + +getConclusion :: Monad m => ComponentId -> ComponentSet -> m Bool +getConclusion n = get n >=> isConclusion + +getRejectionHypotheses :: Monad m => ComponentId -> ComponentSet -> m RejectionHypotheses +getRejectionHypotheses n = get n >=> isRejectionHypotheses + +------------------------------------------------------------------------------ +-- Membership + +contains :: ComponentSet -> ComponentId -> Bool +contains cs x = isJust (get x cs) + +doesNotContain :: ComponentSet -> ComponentId -> Bool +doesNotContain cs = not . contains cs + +isInitial :: ComponentId -> ComponentSet -> Bool +isInitial n cs = initials cs `contains` n + +isDerived :: ComponentId -> ComponentSet -> Bool +isDerived n cs = derived cs `contains` n diff --git a/src/Domain/Statistics/Data.hs b/src/Domain/Statistics/Data.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Statistics/Data.hs @@ -0,0 +1,86 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Statistics.Data + ( Data (..) + , isSingleSample + , sampleSize + , sampleSum + , sampleMean + , sampleVariance + ) where + +import Ideas.Common.Rewriting +import Domain.Statistics.Symbols + +-- | Data component +data Data + = SingleSample [Double] -- Single sample + | TwoSample [Double] [Double] Bool -- Two samples, either paired or not + deriving (Eq, Show, Read) + +instance IsTerm Data where + toTerm (SingleSample xs) = unary singleSampleSymbol (toTerm xs) + toTerm (TwoSample xs ys b) = + ternary twoSampleSymbol (toTerm xs) (toTerm ys) (toTermBool b) + + fromTerm (TCon s [x]) + | s == singleSampleSymbol = + SingleSample <$> fromTerm x + fromTerm (TCon s [xs, ys, b]) + | s == twoSampleSymbol = + TwoSample <$> fromTerm xs <*> fromTerm ys <*> fromTermBool b + fromTerm t = fail $ "Invalid term: " ++ show t ++ " not Data" + +-- | Utility functions +isSingleSample :: Data -> Bool +isSingleSample (SingleSample _) = True +isSingleSample _ = False + +-- | +sampleSize :: Data -> Int +sampleSize (SingleSample xs) = length xs +sampleSize (TwoSample xs _ True) = length xs -- Both samples should have the same sample size +sampleSize _ = error "two sample sizes for unpaired samples" + +-- | +sampleSum :: Data -> [Double] +sampleSum (SingleSample xs) = [sum xs] +sampleSum (TwoSample xs ys _) = [sum xs, sum ys] + +-- | +sampleMean :: Data -> [Double] +sampleMean (SingleSample xs) = [ sum xs / fromIntegral (length xs) ] +sampleMean (TwoSample xs ys _) = [ sum xs / fromIntegral (length xs) + , sum ys / fromIntegral (length ys) + ] + +-- | +sampleVariance :: Bool -> Data -> [Double] +sampleVariance True = sampleVarianceT +sampleVariance False = sampleVarianceZ + +-- | +sampleVarianceZ :: Data -> [Double] +sampleVarianceZ (SingleSample xs) = + let mean = head $ sampleMean (SingleSample xs) + in [foldr (\x r -> r + (x - mean) ** 2) 0 xs / fromIntegral (length xs)] +sampleVarianceZ (TwoSample xs ys _) = + let (mean1:mean2:_) = sampleMean (TwoSample xs ys False) + in [ foldr (\x r -> r + (x - mean1) ** 2) 0 xs / fromIntegral (length xs) + , foldr (\x r -> r + (x - mean2) ** 2) 0 ys / fromIntegral (length ys) + ] + +-- | +sampleVarianceT :: Data -> [Double] +sampleVarianceT (SingleSample xs) = + let mean = head $ sampleMean (SingleSample xs) + in [foldr (\x r -> r + (x - mean) ** 2) 0 xs / fromIntegral (length xs - 1)] +sampleVarianceT (TwoSample xs ys _) = + let (mean1:mean2:_) = sampleMean (TwoSample xs ys False) + in [ foldr (\x r -> r + (x - mean1) ** 2) 0 xs / fromIntegral (length xs - 1) + , foldr (\x r -> r + (x - mean2) ** 2) 0 ys / fromIntegral (length ys - 1) + ] + diff --git a/src/Domain/Statistics/Parser.hs b/src/Domain/Statistics/Parser.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Statistics/Parser.hs @@ -0,0 +1,71 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Statistics.Parser (parseComponentSet) where + +import Domain.Math.Expr +import Domain.Statistics.ComponentSet +import Ideas.Utils.Parsing hiding (char) +import qualified Ideas.Utils.Parsing as P + +parseComponentSet :: String -> Either String ComponentSet +parseComponentSet = parseSimple pComponentSet + +pComponentSet :: Parser ComponentSet +pComponentSet = mconcat <$> many pComponent + +pComponent :: Parser ComponentSet +pComponent = do + (k, b) <- pKey + char ':' + v <- pValue + char ';' + return $ if b then initialSet [(k, v)] else derivedSet [(k, v)] + +pKey :: Parser (ComponentId, Bool) +pKey = (\n -> (n, True)) <$> brackets pComponentId + <|> (\n -> (n, False)) <$> pComponentId + +pComponentId :: Parser ComponentId +pComponentId = fromCId <$> parseName + +pValue :: Parser Component +pValue = CChoice <$> try parseChoice + <|> CRelation <$> try pRelExpr + <|> CExpr <$> try pExpr + +brackets :: Parser a -> Parser a +brackets p = do + char '[' + a <- p + char ']' + return a + +parseChoice :: Parser Choice +parseChoice = choice [ x <$ try (string (show x)) | x <- choices ] + +parseName :: Parser String +parseName = + let + firstChar = oneOf $ ['a'..'z'] ++ ['A'..'Z'] + nonFirstChar = oneOf $ ['a'..'z'] ++ ['A'..'Z'] ++ ['0'..'9'] ++ "-" + in do + fc <- firstChar + rest <- many nonFirstChar + return (fc:rest) + +---------------------------------------------------------------------------- + +char :: Char -> Parser () +char c = lexeme (P.char c) >> return () + +whitespace :: Parser () +whitespace = many (oneOf " \n\t") >> return () + +lexeme :: Parser a -> Parser a +lexeme p = do + x <- p + whitespace + return x diff --git a/src/Domain/Statistics/Rules.hs b/src/Domain/Statistics/Rules.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Statistics/Rules.hs @@ -0,0 +1,94 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Statistics.Rules + ( addNRule + , addAverageRule + , addVarianceRule + , addStandardDeviationRule + , addStandardErrorRule + ) where + +import Control.Monad +import Domain.Math.Expr +import Domain.Statistics.ComponentSet +import Ideas.Common.Library +import Domain.Statistics.Data + +---------------------------------------------------------- +-- Rules on the data component + +addNRule :: Rule ComponentSet +addNRule = describe "Rule for computing the n of the dataset" $ makeRule "component.n" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do guard (not $ contains cs SampleSize) + sample <- getData DataSet cs + let nCs = nComponents sample + return $ foldr (\(nm, c) cs' -> append nm c cs') cs nCs + nComponents (SingleSample xs) = [(SampleSize, CExpr $ Nat $ toInteger (length xs))] + nComponents (TwoSample xs ys _) = [(One SampleSize, CExpr $ Nat $ toInteger (length xs)), + (Two SampleSize, CExpr $ Nat $ toInteger (length ys))] + +addAverageRule :: Rule ComponentSet +addAverageRule = describe "Rule for adding the formula for computing the mean" $ makeRule "component.mean" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do guard (not $ contains cs SampleMean) + guard (not $ contains cs (One SampleMean)) + guard (not $ contains cs (Two SampleMean)) + sample <- getData DataSet cs + -- n <- getExpr SampleMean cs + case sampleMean sample of + [mean] -> return $ append SampleMean (CExpr $ fromDouble mean) cs + [mean0, mean1] -> return $ append (Two SampleMean) (CExpr $ fromDouble mean1) (append (One SampleMean) (CExpr $ fromDouble mean0) cs) + _ -> Nothing + + +addVarianceRule :: Rule ComponentSet +addVarianceRule = describe "Rule for adding the formula for computing the variance" $ makeRule "component.variance" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do guard (not $ contains cs SampleVariance) + guard (not $ contains cs (One SampleVariance)) + guard (not $ contains cs (Two SampleVariance)) + sample <- getData DataSet cs + -- n <- getExpr SampleSize cs + -- mean <- getExpr SampleMean cs + test <- getTestType TestChoice cs + -- s = sum((X_i - avg(X))^2) / n + case sampleVariance (test /= ZTest) sample of + [var] -> return $ append SampleVariance (CExpr $ fromDouble var) cs + [var0, var1] -> return $ append (Two SampleVariance) (CExpr $ fromDouble var1) (append (One SampleVariance) (CExpr $ fromDouble var0) cs) + _ -> Nothing + +addStandardDeviationRule :: Rule ComponentSet +addStandardDeviationRule = describe "Rule for adding the formula for computing the deviation" $ makeRule "component.standard-deviation" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do guard (not $ contains cs SampleSdev) + guard (not $ contains cs (One SampleSdev)) + guard (not $ contains cs (Two SampleSdev)) + sample <- getData DataSet cs + -- n <- getExpr SampleSize cs + -- mean <- getExpr SampleMean cs + test <- getTestType TestChoice cs + -- s = sqrt(sum((X_i - avg(X))^2) / n) + case map sqrt $ sampleVariance (test /= ZTest) sample of + [sd] -> return $ append SampleSdev (CExpr $ fromDouble sd) cs + [sd0, sd1] -> return $ append (Two SampleSdev) (CExpr $ fromDouble sd1) (append (One SampleSdev) (CExpr $ fromDouble sd0) cs) + _ -> Nothing + +addStandardErrorRule :: Rule ComponentSet +addStandardErrorRule = describe "Rule for adding the formula for computing the standard error" $ makeRule "component.standard-error" f + where + f :: ComponentSet -> Maybe ComponentSet + f cs = do guard (not $ contains cs StandardError) + guard (not $ contains cs (One StandardError)) + guard (not $ contains cs (Two StandardError)) + sd <- getExpr SampleSdev cs + -- n <- getExpr SampleSize cs + -- SE = s / sqrt(n) + return $ append StandardError (CExpr $ sd / sqrt(Var "n")) cs diff --git a/src/Domain/Statistics/Symbols.hs b/src/Domain/Statistics/Symbols.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Statistics/Symbols.hs @@ -0,0 +1,133 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +{- Implements symbols for the statistics domain. The only symbols that are + - re-used are the true and the false symbols from the Logic1 module. + -} +module Domain.Statistics.Symbols + ( -- cChoiceSymbol + cConstSymbol + , cDataSymbol + , componentSymbol + -- , conclusionSymbol + , derivedSymbol + , falseSymbol + , initialSymbol + , leftSidedSymbol + , rightSidedSymbol + -- , sidedSymbol + , singleSampleSymbol + , tTestOneSymbol + , tTestPairedSymbol + , tTestTwoSymbol + , rPearsonSymbol + , anovaSymbol + , chiSquaredSymbol + , dependentSymbol, independentSymbol + -- , testTypeSymbol + , toOMSymbol + , trueSymbol + , twoSampleSymbol + , twoSidedSymbol + , zTestSymbol + , rejectH0Symbol, acceptH0Symbol, dontRejectH0Symbol, rejectH1Symbol + , acceptH1Symbol, dontRejectH1Symbol + , toTermBool, fromTermBool + ) where + +import Ideas.Common.Rewriting.Term hiding (trueSymbol, falseSymbol) +import qualified Ideas.Text.OpenMath.Dictionary.Logic1 as L +import Ideas.Common.Id +import qualified Ideas.Text.OpenMath.Symbol as OM + +componentSymbol :: Symbol +componentSymbol = newSymbol "stats.component" + +-- cChoiceSymbol -- component data symbol +cConstSymbol -- component expression symbol, filtered in OM output + , cDataSymbol -- component relation symbol, filtered in OM output +-- , conclusionSymbol -- All of these symbols up to intialSymbol are tied to a + , derivedSymbol -- constructor from Domain.Statistics.ComponentSet + , initialSymbol -- + , leftSidedSymbol -- + , rightSidedSymbol -- +-- , sidedSymbol -- + , singleSampleSymbol -- + , tTestOneSymbol -- + , tTestPairedSymbol -- + , tTestTwoSymbol -- when the step is in the initial set +-- , testTypeSymbol -- when the step is derived + , twoSampleSymbol + , twoSidedSymbol + , zTestSymbol + , rPearsonSymbol + , anovaSymbol + , chiSquaredSymbol + , dependentSymbol + , independentSymbol :: Symbol + +cDataSymbol = newSymbol "stats.cdata" +-- cChoiceSymbol = newSymbol "stats.cchoice" +cConstSymbol = newSymbol "stats.cconst" + +--testTypeSymbol = newSymbol "stats.testtype" +--sidedSymbol = newSymbol "stats.sided" +--conclusionSymbol = newSymbol "stats.conclusion" + +tTestOneSymbol = newSymbol "stats.ttestone" +tTestTwoSymbol = newSymbol "stats.ttesttwo" +tTestPairedSymbol = newSymbol "stats.ttestpaired" +zTestSymbol = newSymbol "stats.ztest" +rPearsonSymbol = newSymbol "stats.rpearson" +anovaSymbol = newSymbol "stats.anova" +chiSquaredSymbol = newSymbol "stats.chisquared" + +twoSidedSymbol = newSymbol "stats.twosided" +leftSidedSymbol = newSymbol "stats.leftsided" +rightSidedSymbol = newSymbol "stats.rightsided" + +singleSampleSymbol = newSymbol "stats.singlesample" +twoSampleSymbol = newSymbol "stats.twosample" + +dependentSymbol = newSymbol "stats.dependent" +independentSymbol = newSymbol "stats.independent" + +initialSymbol = newSymbol "stats.initial" +derivedSymbol = newSymbol "stats.derived" + +trueSymbol, falseSymbol :: Symbol +trueSymbol = newSymbol L.trueSymbol +falseSymbol = newSymbol L.falseSymbol + +toOMSymbol :: Symbol -> OM.Symbol +toOMSymbol s = idToSymbol (getId s) + +rejectH0Symbol, acceptH0Symbol, dontRejectH0Symbol, rejectH1Symbol, + acceptH1Symbol, dontRejectH1Symbol :: Symbol +rejectH0Symbol = newSymbol "stats.rejecth0" +acceptH0Symbol = newSymbol "stats.accepth0" +dontRejectH0Symbol = newSymbol "stats.dontrejecth0" +rejectH1Symbol = newSymbol "stats.rejecth1" +acceptH1Symbol = newSymbol "stats.accepth1" +dontRejectH1Symbol = newSymbol "stats.dontrejecth1" + +-- Copied from Ideas.Encoding.OpenMathSuppoer +idToSymbol :: Id -> OM.Symbol +idToSymbol a + | null (qualifiers a) = + OM.extraSymbol (unqualified a) + | otherwise = + OM.makeSymbol (qualification a) (unqualified a) + +toTermBool :: Bool -> Term +toTermBool b = symbol s + where + s = if b then trueSymbol else falseSymbol + +fromTermBool :: Monad m => Term -> m Bool +fromTermBool (TCon s []) + | s == trueSymbol = return True + | s == falseSymbol = return False +fromTermBool _ = fail "Term is not a boolean" diff --git a/src/Domain/Statistics/Views.hs b/src/Domain/Statistics/Views.hs new file mode 100644 --- /dev/null +++ b/src/Domain/Statistics/Views.hs @@ -0,0 +1,183 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Domain.Statistics.Views + ( Substitution, Substitutable(..), single + , getSubstitution, substitute + , checkSubstitution + , substitutedView, evaluatedView + ) where + +import qualified Data.Map as M +import qualified Data.Set as S +import Data.List (intercalate) +import Data.Maybe +import Domain.Math.Data.Relation +import Domain.Math.Expr +import Domain.Math.Numeric.Views +import Domain.Statistics.ComponentSet +import Ideas.Common.View +import Ideas.Utils.Uniplate + +-- idempotent, non-recursive substitutions +newtype Substitution = S (M.Map String Expr) + +instance Show Substitution where + show (S m) = "{" ++ intercalate ", " (map f (M.toList m)) ++ "}" + where + f (s, expr) = s ++ "=" ++ show expr + +instance Semigroup Substitution where + (<>) (S m) = flip (M.foldrWithKey extend) m + +instance Monoid Substitution where + mempty = S M.empty + mappend = (<>) + +member :: String -> Substitution -> Bool +member s (S m) = s `M.member` m + +single :: String -> Expr -> Substitution +single s expr + | s `elem` vars expr = mempty + | otherwise = S (M.singleton s expr) + +extend :: String -> Expr -> Substitution -> Substitution +extend s expr sub@(S m) + | s `elem` vars expr = mempty -- recursive expression + | otherwise = + case M.lookup s m of + Nothing -> S (M.insert s (sub |-> expr) (M.map (single s expr |->) m)) + Just e + | e == (sub |-> expr) -> sub + -- inconsistent extension; should not happen. + -- deleting s in m breaks monoid property for Substitution + | otherwise -> S (M.delete s m) + +-- throws an error for non-idempotent substitutions +checkSubstitution :: Substitution -> Substitution +checkSubstitution sub@(S m) + | disjoint (M.keysSet m) (vars (M.elems m)) = sub + | otherwise = error $ "Non-idempotent substitution: " ++ show sub + where + -- in Data.Set, since 0.5.11 + disjoint x y = S.null (S.intersection x y) + +class Substitutable a where + (|->) :: Substitution -> a -> a + vars :: a -> S.Set String + +instance Substitutable a => Substitutable [a] where + (|->) = map . (|->) + vars = S.unions . map vars + +instance Substitutable Expr where + S xs |-> Var s = fromMaybe (Var s) (M.lookup s xs) + sub |-> expr = descend (sub |->) expr + + vars (Var s) = S.singleton s + vars expr = S.unions (map vars (children expr)) + +instance Substitutable a => Substitutable (Relation a) where + sub |-> r = fmap (sub |->) r + vars r = vars (leftHandSide r) `S.union` vars (rightHandSide r) + +instance Substitutable Component where + sub |-> CExpr e = CExpr $ sub |-> e + sub |-> CRelation r = CRelation $ sub |-> r + _ |-> c = c + + vars (CExpr e) = vars e + vars (CRelation r) = vars r + vars _ = S.empty + +instance Substitutable ComponentSet where + sub |-> cs = mapComponent (sub |->) cs + vars = S.unions . map (vars . snd) . toList + +getSubstitution :: ComponentSet -> Substitution +getSubstitution = addSigmaM . addSEM . make . resolveTestFormulaValue . addComponentVars + where + -- if test value and test formula are both present (but differ because of + -- rounding), prefer the test value + resolveTestFormulaValue cs + | cs `contains` TestValue && cs `contains` TestFormula = + delete TestFormula cs + | otherwise = cs + + make cs = mconcat (map f (xs ++ ys)) + where + -- prefer initials in substitution + xs = toList (initials cs) + ys = filter notInital (toList (derived cs)) + notInital = (`notElem` map fst xs) . fst + + f (cid, CExpr e) = single (getCId cid) e + f (_, CRelation r) | relationType r == EqualTo = + case leftHandSide r of + Var v -> single v (rightHandSide r) + _ -> mempty + f _ = mempty + + -- add formulas, but only if not already present in the substitution + addSigmaM sub + | "sigmaM" `member` sub = sub + | otherwise = sub <> single "sigmaM" (Var "sigma" / sqrt (Var "n")) + addSEM sub + | "SEM" `member` sub = sub + | otherwise = sub <> single "SEM" (Var "s" / sqrt (Var "n")) + +substitute :: ComponentSet -> ComponentSet +substitute cs = getSubstitution cs |-> cs + +addComponentVars :: ComponentSet -> ComponentSet +addComponentVars cs = foldr (uncurry addComponentVar) (addPopulationMean cs) table + where + table = + [ (PopulationSdev, "sigma") + , (SampleSdev, "s") + , (SampleMean, "M") + , (One SampleMean, "M1") + , (Two SampleMean, "M2") + , (One SampleSdev, "s1") + , (Two SampleSdev, "s2") + , (PopulationMean, "mu") + ] + +-- extract population mean from alternative hypothesis +addPopulationMean :: ComponentSet -> ComponentSet +addPopulationMean cs = + case (get PopulationMean cs, getRelation AlternativeHypothesis cs) of + (Nothing, Just rel) -> append PopulationMean (CExpr (rightHandSide rel)) cs + _ -> cs + +addComponentVar :: ComponentId -> String -> ComponentSet -> ComponentSet +addComponentVar cid var cs = + case getRhsExpr cid cs of + Just e | cs `doesNotContain` Other var -> + -- to do: is it necessary to add (Other var) to the component set?? + append (Other var) (CExpr e) (single var e |-> cs) + _ -> cs + +---------------------------------------------------------- +-- Substituted view + +substitutedView :: View ComponentSet ComponentSet +substitutedView = makeView (Just . substitute) id + +---------------------------------------------------------- +-- Evaluated view + +evaluatedView :: View ComponentSet ComponentSet +evaluatedView = makeView f id + where + f :: ComponentSet -> Maybe ComponentSet + f cs = Just $ mapComponent evaluateComponent cs + +evaluateComponent :: Component -> Component +evaluateComponent (CExpr x) = CExpr $ simplify doubleView x +evaluateComponent (CRelation r) = CRelation $ fmap (simplify doubleView) r +evaluateComponent c = c + diff --git a/src/Main.hs b/src/Main.hs new file mode 100644 --- /dev/null +++ b/src/Main.hs @@ -0,0 +1,38 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Main where + +-- Include the Ideas framework +import Ideas.Common.Library +import Ideas.Main.Default + +-- Include the hypothesis exercise and predicate service +import Domain.Hypothesis.Exercises +import Service.ConstraintServices + +main :: IO () +main = defaultMain ideasStatistics + +ideasStatistics :: DomainReasoner +ideasStatistics = + describe "Domain reasoner for statistics" $ + (newDomainReasoner "ideas.statistics") + { exercises = myExercises + , services = myServices + } + +myExercises :: [Some Exercise] +myExercises = + [ -- hypothesis testing + Some hypothesisExercise + ] + +myServices :: [Service] +myServices = + concat [ metaServiceList ideasStatistics + , filter ((/= newId "basic.diagnose") . getId) serviceList + , constraintServiceList + ] diff --git a/src/Service/ConstraintServices.hs b/src/Service/ConstraintServices.hs new file mode 100644 --- /dev/null +++ b/src/Service/ConstraintServices.hs @@ -0,0 +1,93 @@ +----------------------------------------------------------------------------- +-- Copyright 2020, Ideas project team. This file is distributed under the +-- terms of the Apache License 2.0. For more information, see the files +-- "LICENSE.txt" and "NOTICE.txt", which are included in the distribution. +----------------------------------------------------------------------------- +module Service.ConstraintServices + ( constraintServiceList + ) where + +import Ideas.Common.Library +import Ideas.Service.Diagnose +import Ideas.Service.BasicServices +import Ideas.Service.State +import Ideas.Service.Types + +-- Constraint services that are available +constraintServiceList :: [Service] +constraintServiceList = [diagnoseS, constraintsS, onehintS, allhintsS] + +diagnoseS :: Service +diagnoseS = makeService "basic.diagnose" + "Diagnose an expression submitted by a student. Possible diagnosis are \ + \Buggy (a common misconception was detected), NotEquivalent (something is \ + \wrong, but we don't know what), Similar (the expression is pretty similar \ + \to the last expression in the derivation), Expected (the submitted \ + \expression was anticipated by the strategy), Detour (the submitted \ + \expression was not expected by the strategy, but the applied rule was \ + \detected), and Correct (it is correct, but we don't know which rule was \ + \applied). Extended version for statistics domain: check predicates for NotEquiv." $ + diagnoseWithConstraints ::: tState .-> tContext .-> tMaybe tId .-> tDiagnosis + +constraintsS :: Service +constraintsS = makeService "basic.constraints" + "Check all constraints" $ + checkConstraints ::: tState .-> tList (tPair tConstraint tString) + +onehintS :: Service +onehintS = makeService "basic.onehint" + "Check constraints, and provide a hint" $ + onehint ::: tState .-> tError (Tag "elem" tHint) + +allhintsS :: Service +allhintsS = makeService "basic.allhints" + "Check constraints, and provide (multiple) hints" $ + allhints ::: tState .-> tError (tList tHint) + +-------------------------------------------------------------------------------- + +diagnoseWithConstraints :: State a -> Context a -> Maybe Id -> Diagnosis a +diagnoseWithConstraints st ctx mid = f (diagnose st ctx mid) + where + f (NotEquivalent "") = + case violations (exercise st) ctx of + (n, msg):_ -> NotEquivalent $ violationMsg n msg + [] | null (violations (exercise st) (stateContext st)) -> + NotEquivalent [] + | otherwise -> -- special case: previous state is invalid (constraint violated) + Correct (finished st) st + f d = d + +data Hint a = Violation (Constraint (Context a)) String + | HintStep (StepInfo a) (State a) + +tHint :: Type a (Hint a) +tHint = Iso (f <-> g) tp + where + tp = tPair tConstraint tString :|: tPair tStepInfo tState + + f (Left (c, msg)) = Violation c msg + f (Right (stp, st)) = HintStep stp st + g (Violation c msg) = Left (c, msg) + g (HintStep stp st) = Right (stp, st) + +allhints :: State a -> Either String [Hint a] +allhints st = + case violations (exercise st) (stateContext st) of + [] -> fmap (map (uncurry HintStep)) (allfirsts st) + xs -> Right (map (uncurry Violation) xs) + +onehint :: State a -> Either String (Hint a) +onehint st = + case allhints st of + Left msg -> Left msg + Right [] -> Left "no hint available" + Right (h:_) -> Right h + +checkConstraints :: State a -> [(Constraint (Context a), String)] +checkConstraints st = map f (constraints (exercise st)) + where + f c = (c, show $ getResult c (stateContext st)) + +violationMsg :: Constraint a -> String -> String +violationMsg n msg = show n ++ if null msg then "" else ": " ++ msg