diff --git a/README.md b/README.md index a9e55825..618bc7a4 100644 --- a/README.md +++ b/README.md @@ -279,3 +279,6 @@ _____________________________________________________________ Last changed on: Source code: https://github.com/smartell/iSCAM _____________________________________________________________ + +## Code Changes + - March 30, 2015. Added time varying catchbility options (random walk). diff --git a/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl b/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl index 94eaed55..b92d97ab 100644 --- a/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl +++ b/fba/HALIBUT/DATA/DMVL1988/Halibut2014.ctl @@ -83,7 +83,7 @@ 1888 1888 1888 -## +## ## ———————————————————————————————————————————————————————————————————————————————————— ## ## TIME VARYING NATURAL MORTALIIY RATES ## ## ———————————————————————————————————————————————————————————————————————————————————— ## @@ -104,18 +104,20 @@ ## ———————————————————————————————————————————————————————————————————————————————————— ## ## ## -## ------------------------------------------------------------------------- ## -## PRIORS FOR SURVEY Q ## -## Prior type: ## -## 0 - uninformative prior ## -## 1 - normal prior density for log(q) ## -## 2 - random walk in q ## -## ------------------------------------------------------------------------- ## -2 # -number of surveys (nits) -0 0 # -prior type (see legend above -0 0 # -prior log(mean) -0 0 # -prior sd -## ------------------------------------------------------------------------- ## +## ———————————————————————————————————————————————————————————————————————————————————— ## +## ABUNDANCE OBSERVATION MODELS +## ———————————————————————————————————————————————————————————————————————————————————— ## +## QTYPE: +## 0 = FIXED SURVEY Q (specify log(mean) for prior log(mean)) +## 1 = CONSTANT Q (use MLE for q and optional informative prior) +## 2 = RANDOM WALK Q (use prior mean & sd for penalized random walk) +## ———————————————————————————————————————————————————————————————————————————————————— ## +2 # -number of surveys (n_it_nobs) + 2 2 # -QTYPE (see legend above) + 0 0 # -prior log(mean) + 0 0 # -prior sd (set to 0 for uniformative prior) + 1 1 # -Estimation phase +## ———————————————————————————————————————————————————————————————————————————————————— ## ## ## ------------------------------------------------------------------------- ## diff --git a/fba/HALIBUT/DATA/Halibut2014.dat b/fba/HALIBUT/DATA/Halibut2014.dat index 42cd972d..9715bca7 100644 --- a/fba/HALIBUT/DATA/Halibut2014.dat +++ b/fba/HALIBUT/DATA/Halibut2014.dat @@ -384,123 +384,123 @@ ## commercial_data: ## survey ## iyr index(it) gear area group sex log(SE) log(PE) timing - 1907 280 1 1 1 0 0.100 0.010 0.5 - 1910 271 1 1 1 0 0.100 0.010 0.5 - 1911 237 1 1 1 0 0.100 0.010 0.5 - 1912 176 1 1 1 0 0.100 0.010 0.5 - 1913 129 1 1 1 0 0.100 0.010 0.5 - 1914 124 1 1 1 0 0.100 0.010 0.5 - 1915 118 1 1 1 0 0.100 0.010 0.5 - 1916 137 1 1 1 0 0.100 0.010 0.5 - 1917 98 1 1 1 0 0.100 0.010 0.5 - 1918 96 1 1 1 0 0.100 0.010 0.5 - 1919 93 1 1 1 0 0.100 0.010 0.5 - 1920 96 1 1 1 0 0.100 0.010 0.5 - 1921 88 1 1 1 0 0.100 0.010 0.5 - 1922 73 1 1 1 0 0.100 0.010 0.5 - 1923 78 1 1 1 0 0.100 0.010 0.5 - 1924 74 1 1 1 0 0.100 0.010 0.5 - 1925 68 1 1 1 0 0.100 0.010 0.5 - 1926 67 1 1 1 0 0.100 0.010 0.5 - 1927 65 1 1 1 0 0.100 0.010 0.5 - 1928 58 1 1 1 0 0.100 0.010 0.5 - 1929 51 1 1 1 0 0.100 0.010 0.5 - 1930 46 1 1 1 0 0.100 0.010 0.5 - 1931 50 1 1 1 0 0.100 0.010 0.5 - 1932 60 1 1 1 0 0.100 0.010 0.5 - 1933 63 1 1 1 0 0.100 0.010 0.5 - 1934 62 1 1 1 0 0.100 0.010 0.5 - 1935 76 1 1 1 0 0.100 0.010 0.5 - 1936 71 1 1 1 0 0.100 0.010 0.5 - 1937 80 1 1 1 0 0.100 0.010 0.5 - 1938 88 1 1 1 0 0.100 0.010 0.5 - 1939 80 1 1 1 0 0.100 0.010 0.5 - 1940 81 1 1 1 0 0.100 0.010 0.5 - 1941 85 1 1 1 0 0.100 0.010 0.5 - 1942 90 1 1 1 0 0.100 0.010 0.5 - 1943 95 1 1 1 0 0.100 0.010 0.5 - 1944 110 1 1 1 0 0.100 0.010 0.5 - 1945 102 1 1 1 0 0.100 0.010 0.5 - 1946 101 1 1 1 0 0.100 0.010 0.5 - 1947 99 1 1 1 0 0.100 0.010 0.5 - 1948 99 1 1 1 0 0.100 0.010 0.5 - 1949 95 1 1 1 0 0.100 0.010 0.5 - 1950 95 1 1 1 0 0.100 0.010 0.5 - 1951 96 1 1 1 0 0.100 0.010 0.5 - 1952 110 1 1 1 0 0.100 0.010 0.5 - 1953 131 1 1 1 0 0.100 0.010 0.5 - 1954 133 1 1 1 0 0.100 0.010 0.5 - 1955 119 1 1 1 0 0.100 0.010 0.5 - 1956 129 1 1 1 0 0.100 0.010 0.5 - 1957 110 1 1 1 0 0.100 0.010 0.5 - 1958 121 1 1 1 0 0.100 0.010 0.5 - 1959 129 1 1 1 0 0.100 0.010 0.5 - 1960 132 1 1 1 0 0.100 0.010 0.5 - 1961 127 1 1 1 0 0.100 0.010 0.5 - 1962 115 1 1 1 0 0.100 0.010 0.5 - 1963 105 1 1 1 0 0.100 0.010 0.5 - 1964 100 1 1 1 0 0.100 0.010 0.5 - 1965 99 1 1 1 0 0.100 0.010 0.5 - 1966 100 1 1 1 0 0.100 0.010 0.5 - 1967 101 1 1 1 0 0.100 0.010 0.5 - 1968 103 1 1 1 0 0.100 0.010 0.5 - 1969 95 1 1 1 0 0.100 0.010 0.5 - 1970 91 1 1 1 0 0.100 0.010 0.5 - 1971 89 1 1 1 0 0.100 0.010 0.5 - 1972 78 1 1 1 0 0.100 0.010 0.5 - 1973 63 1 1 1 0 0.100 0.010 0.5 - 1974 61 1 1 1 0 0.100 0.010 0.5 - 1975 61 1 1 1 0 0.100 0.010 0.5 - 1976 55 1 1 1 0 0.100 0.010 0.5 - 1977 63 1 1 1 0 0.100 0.010 0.5 - 1978 71 1 1 1 0 0.100 0.010 0.5 - 1979 75 1 1 1 0 0.100 0.010 0.5 - 1980 94 1 1 1 0 0.100 0.010 0.5 - 1981 111 1 1 1 0 0.100 0.010 0.5 - 1982 127 1 1 1 0 0.100 0.010 0.5 - 1984 291 1 1 1 0 0.100 0.010 0.5 - 1985 357 1 1 1 0 0.036 0.010 0.5 - 1986 320 1 1 1 0 0.041 0.010 0.5 - 1987 321 1 1 1 0 0.040 0.010 0.5 - 1988 368 1 1 1 0 0.035 0.010 0.5 - 1989 358 1 1 1 0 0.025 0.010 0.5 - 1990 319 1 1 1 0 0.028 0.010 0.5 - 1991 318 1 1 1 0 0.038 0.010 0.5 - 1992 319 1 1 1 0 0.034 0.010 0.5 - 1993 373 1 1 1 0 0.099 0.010 0.5 - 1994 306 1 1 1 0 0.072 0.010 0.5 - 1995 330 1 1 1 0 0.036 0.010 0.5 - 1996 392 1 1 1 0 0.038 0.010 0.5 - 1997 405 1 1 1 0 0.025 0.010 0.5 - 1998 408 1 1 1 0 0.025 0.010 0.5 - 1999 393 1 1 1 0 0.023 0.010 0.5 - 2000 403 1 1 1 0 0.022 0.010 0.5 - 2001 362 1 1 1 0 0.041 0.010 0.5 - 2002 360 1 1 1 0 0.019 0.010 0.5 - 2003 329 1 1 1 0 0.018 0.010 0.5 - 2004 319 1 1 1 0 0.019 0.010 0.5 - 2005 296 1 1 1 0 0.017 0.010 0.5 - 2006 270 1 1 1 0 0.019 0.010 0.5 - 2007 252 1 1 1 0 0.020 0.010 0.5 - 2008 232 1 1 1 0 0.017 0.010 0.5 - 2009 222 1 1 1 0 0.018 0.010 0.5 - 2010 204 1 1 1 0 0.020 0.010 0.5 - 2011 198 1 1 1 0 0.015 0.010 0.5 - 2012 195 1 1 1 0 0.021 0.010 0.5 - 2013 179 1 1 1 0 0.017 0.010 0.5 - 2014 191 1 1 1 0 0.052 0.010 0.5 + 1907 280 1 1 1 0 0.100 0.050 0.5 + 1910 271 1 1 1 0 0.100 0.050 0.5 + 1911 237 1 1 1 0 0.100 0.050 0.5 + 1912 176 1 1 1 0 0.100 0.050 0.5 + 1913 129 1 1 1 0 0.100 0.050 0.5 + 1914 124 1 1 1 0 0.100 0.050 0.5 + 1915 118 1 1 1 0 0.100 0.050 0.5 + 1916 137 1 1 1 0 0.100 0.050 0.5 + 1917 98 1 1 1 0 0.100 0.050 0.5 + 1918 96 1 1 1 0 0.100 0.050 0.5 + 1919 93 1 1 1 0 0.100 0.050 0.5 + 1920 96 1 1 1 0 0.100 0.050 0.5 + 1921 88 1 1 1 0 0.100 0.050 0.5 + 1922 73 1 1 1 0 0.100 0.050 0.5 + 1923 78 1 1 1 0 0.100 0.050 0.5 + 1924 74 1 1 1 0 0.100 0.050 0.5 + 1925 68 1 1 1 0 0.100 0.050 0.5 + 1926 67 1 1 1 0 0.100 0.050 0.5 + 1927 65 1 1 1 0 0.100 0.050 0.5 + 1928 58 1 1 1 0 0.100 0.050 0.5 + 1929 51 1 1 1 0 0.100 0.050 0.5 + 1930 46 1 1 1 0 0.100 0.050 0.5 + 1931 50 1 1 1 0 0.100 0.050 0.5 + 1932 60 1 1 1 0 0.100 0.050 0.5 + 1933 63 1 1 1 0 0.100 0.050 0.5 + 1934 62 1 1 1 0 0.100 0.050 0.5 + 1935 76 1 1 1 0 0.100 0.050 0.5 + 1936 71 1 1 1 0 0.100 0.050 0.5 + 1937 80 1 1 1 0 0.100 0.050 0.5 + 1938 88 1 1 1 0 0.100 0.050 0.5 + 1939 80 1 1 1 0 0.100 0.050 0.5 + 1940 81 1 1 1 0 0.100 0.050 0.5 + 1941 85 1 1 1 0 0.100 0.050 0.5 + 1942 90 1 1 1 0 0.100 0.050 0.5 + 1943 95 1 1 1 0 0.100 0.050 0.5 + 1944 110 1 1 1 0 0.100 0.050 0.5 + 1945 102 1 1 1 0 0.100 0.050 0.5 + 1946 101 1 1 1 0 0.100 0.050 0.5 + 1947 99 1 1 1 0 0.100 0.050 0.5 + 1948 99 1 1 1 0 0.100 0.050 0.5 + 1949 95 1 1 1 0 0.100 0.050 0.5 + 1950 95 1 1 1 0 0.100 0.050 0.5 + 1951 96 1 1 1 0 0.100 0.050 0.5 + 1952 110 1 1 1 0 0.100 0.050 0.5 + 1953 131 1 1 1 0 0.100 0.050 0.5 + 1954 133 1 1 1 0 0.100 0.050 0.5 + 1955 119 1 1 1 0 0.100 0.050 0.5 + 1956 129 1 1 1 0 0.100 0.050 0.5 + 1957 110 1 1 1 0 0.100 0.050 0.5 + 1958 121 1 1 1 0 0.100 0.050 0.5 + 1959 129 1 1 1 0 0.100 0.050 0.5 + 1960 132 1 1 1 0 0.100 0.050 0.5 + 1961 127 1 1 1 0 0.100 0.050 0.5 + 1962 115 1 1 1 0 0.100 0.050 0.5 + 1963 105 1 1 1 0 0.100 0.050 0.5 + 1964 100 1 1 1 0 0.100 0.050 0.5 + 1965 99 1 1 1 0 0.100 0.050 0.5 + 1966 100 1 1 1 0 0.100 0.050 0.5 + 1967 101 1 1 1 0 0.100 0.050 0.5 + 1968 103 1 1 1 0 0.100 0.050 0.5 + 1969 95 1 1 1 0 0.100 0.050 0.5 + 1970 91 1 1 1 0 0.100 0.050 0.5 + 1971 89 1 1 1 0 0.100 0.050 0.5 + 1972 78 1 1 1 0 0.100 0.050 0.5 + 1973 63 1 1 1 0 0.100 0.050 0.5 + 1974 61 1 1 1 0 0.100 0.050 0.5 + 1975 61 1 1 1 0 0.100 0.050 0.5 + 1976 55 1 1 1 0 0.100 0.050 0.5 + 1977 63 1 1 1 0 0.100 0.050 0.5 + 1978 71 1 1 1 0 0.100 0.050 0.5 + 1979 75 1 1 1 0 0.100 0.050 0.5 + 1980 94 1 1 1 0 0.100 0.050 0.5 + 1981 111 1 1 1 0 0.100 0.050 0.5 + 1982 127 1 1 1 0 0.100 0.050 0.5 + 1984 291 1 1 1 0 0.100 1.000 0.5 + 1985 357 1 1 1 0 0.036 0.050 0.5 + 1986 320 1 1 1 0 0.041 0.050 0.5 + 1987 321 1 1 1 0 0.040 0.050 0.5 + 1988 368 1 1 1 0 0.035 0.050 0.5 + 1989 358 1 1 1 0 0.025 0.050 0.5 + 1990 319 1 1 1 0 0.028 0.050 0.5 + 1991 318 1 1 1 0 0.038 0.050 0.5 + 1992 319 1 1 1 0 0.034 0.050 0.5 + 1993 373 1 1 1 0 0.099 0.050 0.5 + 1994 306 1 1 1 0 0.072 0.050 0.5 + 1995 330 1 1 1 0 0.036 0.050 0.5 + 1996 392 1 1 1 0 0.038 0.050 0.5 + 1997 405 1 1 1 0 0.025 0.050 0.5 + 1998 408 1 1 1 0 0.025 0.050 0.5 + 1999 393 1 1 1 0 0.023 0.050 0.5 + 2000 403 1 1 1 0 0.022 0.050 0.5 + 2001 362 1 1 1 0 0.041 0.050 0.5 + 2002 360 1 1 1 0 0.019 0.050 0.5 + 2003 329 1 1 1 0 0.018 0.050 0.5 + 2004 319 1 1 1 0 0.019 0.050 0.5 + 2005 296 1 1 1 0 0.017 0.050 0.5 + 2006 270 1 1 1 0 0.019 0.050 0.5 + 2007 252 1 1 1 0 0.020 0.050 0.5 + 2008 232 1 1 1 0 0.017 0.050 0.5 + 2009 222 1 1 1 0 0.018 0.050 0.5 + 2010 204 1 1 1 0 0.020 0.050 0.5 + 2011 198 1 1 1 0 0.015 0.050 0.5 + 2012 195 1 1 1 0 0.021 0.050 0.5 + 2013 179 1 1 1 0 0.017 0.050 0.5 + 2014 191 1 1 1 0 0.052 0.050 0.5 ## survey_data: 1977 1.4713 6 1 1 0 0.100 0.010 0.5 1978 1.1112 6 1 1 0 0.100 0.010 0.5 1980 2.0090 6 1 1 0 0.100 0.010 0.5 1981 2.6670 6 1 1 0 0.100 0.010 0.5 - 1982 2.8739 6 1 1 0 0.100 0.010 0.5 + 1982 2.8739 6 1 1 0 0.100 1.000 0.5 1983 2.8816 6 1 1 0 0.100 0.010 0.5 - 1984 9.3014 6 1 1 0 0.100 0.010 0.5 + 1984 9.3014 6 1 1 0 0.100 1.000 0.5 1985 8.9419 6 1 1 0 0.100 0.010 0.5 1986 6.2568 6 1 1 0 0.100 0.010 0.5 1996 12.8858 6 1 1 0 0.100 0.010 0.5 - 1997 7.7841 6 1 1 0 0.034 0.010 0.5 + 1997 7.7841 6 1 1 0 0.034 1.000 0.5 1998 6.9727 6 1 1 0 0.034 0.010 0.5 1999 6.2649 6 1 1 0 0.034 0.010 0.5 2000 6.6265 6 1 1 0 0.034 0.010 0.5 diff --git a/scripts/Template.ctl b/scripts/Template.ctl index bb28e38d..f3ff2643 100644 --- a/scripts/Template.ctl +++ b/scripts/Template.ctl @@ -100,20 +100,23 @@ ## ———————————————————————————————————————————————————————————————————————————————————— ## ## ## -## ------------------------------------------------------------------------- ## -## PRIORS FOR SURVEY Q ## -## Prior type: ## -## 0 - uninformative prior ## -## 1 - normal prior density for log(q) ## -## 2 - random walk in q ## -## ------------------------------------------------------------------------- ## -1 # -number of surveys (nits) -0 # -prior type (see legend above) -0 # -prior log(mean) -0 # -prior sd -## ------------------------------------------------------------------------- ## +## ———————————————————————————————————————————————————————————————————————————————————— ## +## ABUNDANCE OBSERVATION MODELS +## ———————————————————————————————————————————————————————————————————————————————————— ## +## QTYPE: +## 0 = FIXED SURVEY Q (specify log(mean) for prior log(mean)) +## 1 = CONSTANT Q (use MLE for q and optional informative prior) +## 2 = RANDOM WALK Q (use prior mean & sd for penalized random walk) +## ———————————————————————————————————————————————————————————————————————————————————— ## +1 # -number of surveys (n_it_nobs) + 2 # -QTYPE (see legend above) + 0 # -prior log(mean) + 0 # -prior sd (set to 0 for uniformative prior) + 1 # -Estimation phase +## ———————————————————————————————————————————————————————————————————————————————————— ## ## + ## ------------------------------------------------------------------------- ## ## OTHER MISCELANEOUS CONTROLS ## ## ------------------------------------------------------------------------- ## diff --git a/src/admb-code/include/compositionLikelihoods.hpp b/src/admb-code/include/compositionLikelihoods.hpp new file mode 100644 index 00000000..185a7e82 --- /dev/null +++ b/src/admb-code/include/compositionLikelihoods.hpp @@ -0,0 +1,292 @@ +// negLogLikelihood.hpp +#include + +#ifndef negLogLikelihood_H +#define negLogLikelihood_H + +/** + * @file negLogLikelihood.hpp + * @defgroup Likelihoods + * @author Steven Martell + * @namespace acl + * @date Feb 10, 2014 + * @title Selectivity functions + * @details Uses abstract base class for computing negative loglikelihoods + */ +namespace acl +{ + + /** + * @brief Base class for composition likelihoods. + * @details Virtual methods for negative loglikelihood and standardized residuals. + * + * @tparam DATA Data type argument. + * @tparam DVAR Variable type argument. + */ + template + class negLogLikelihood + { + private: + int r1,r2; /// index for first and last row. + int c1,c2; /// index for rist and last column. + ivector m_jmin; /// index for ragged start columns. + ivector m_jmax; /// index for ragged end columns. + imatrix m_jagg; /// index for aggregation among bins. + + DATA m_O; /// observed composition object + DATA m_rO; /// observed ragged composition object + + DVAR m_P; /// predicted composition object + DVAR m_rP; /// predicted ragged composition object + + + public: + // constructors + negLogLikelihood(){} + negLogLikelihood(const DATA& _O, const DVAR& _P) + :m_O(_O), m_P(_P) + { + r1 = m_O.rowmin(); + r2 = m_O.rowmax(); + c1 = m_O.colmin(); + c2 = m_O.colmax(); + + getRaggedVectors(); + } + + // pure virtual methods + virtual ~negLogLikelihood(){}; + virtual const dvariable nloglike() const = 0; + virtual const DVAR residual() const = 0; + + // virtual methods + inline + const DATA tailCompression(const DATA &_M) const; + + // methods + inline + void getRaggedVectors(double pmin = 0.0); + + inline + void aggregate(const double pmin = 0.0); + + // getters & setters + void set_O(DATA & O) { this -> m_O = O; } + DATA get_O() const{ return m_O; } + + DVAR get_P() const {return m_P; } + void set_P(DVAR &P) {this->m_P=P;} + + DATA get_rO() const {return m_rO; } + void set_rO(DATA &X){this->m_rO=X;} + + DVAR get_rP() const {return m_rP; } + void set_rP(DVAR &X){this->m_rP=X;} + + + template + inline + const T compress(const T& _M) const + { + T R; + T M = _M; + R.allocate(r1,r2,m_jmin,m_jmax); + R.initialize(); + // fill ragged array R + for(int i = r1; i <= r2; i++ ) + { + M(i) /= sum(M(i)); + R(i)(m_jmin(i),m_jmax(i)) = M(i)(m_jmin(i),m_jmax(i)); + + // add cumulative sum to tails. + R(i)(m_jmin(i)) = sum(M(i)(c1,m_jmin(i))); + R(i)(m_jmax(i)) = sum(M(i)(m_jmax(i),c2)); + } + + return R; + } + + }; + + template + inline + void acl::negLogLikelihood::aggregate(const double pmin) + { + m_jmin.allocate(r1,r2); + m_jmax.allocate(r1,r2); + ivector n(r1,r2); + n.initialize(); + + // get number of observations > pmin each year. + for(int i = r1; i <= r2; i++ ) + { + dvector oo = m_O(i)/sum(m_O(i)); + for(int j = c1; j <= c2; j++ ) + { + if( oo(j) > pmin ) n(i)++; + } + } + + m_rO.allocate(r1,r2,1,n); + m_rP.allocate(r1,r2,1,n); + m_rO.initialize(); + m_rP.initialize(); + + for(int i = r1; i <= r2; i++ ) + { + int k = 1; + dvector oo = m_O(i)/sum(m_O(i)); + dvar_vector pp = m_P(i)/sum(m_P(i)); + for(int j = c1; j <= c2; j++ ) + { + m_rO(i)(k) += oo(j); + m_rP(i)(k) += pp(j); + if( oo(j)>pmin && k + inline + void acl::negLogLikelihood::getRaggedVectors(double pmin) + { + m_jmin.allocate(r1,r2); + m_jmax.allocate(r1,r2); + + for(int i = r1; i <= r2; i++ ) + { + dvector o = m_O(i); + + m_jmin(i) = c1; // index for min column + m_jmax(i) = c2; // index for max column + dvector cumsum = o/sum(o); + for(int j = c1; j < c2; j++ ) + { + // min bin and plus bin column index. + cumsum(j) <= pmin ? m_jmin(i)++ : NULL; + 1.0 - cumsum(j) <= pmin ? m_jmax(i)-- : NULL; + cumsum(j+1) += cumsum(j); + + // o(j) > pmin ? n(i)++ : NULL; + } + } + + } + + + + + + + + + + + + + + + + + // |--------------------------------------------------------------------------------| + // | MULTIVARIATE LOGISTIC NEGATIVE LOGLIKELIHOOD derived class | + // |--------------------------------------------------------------------------------| + + /** + * @brief multivariate logistic negative log likelihood + * @details The negative log likeihood for the multivariate + * logistic distribution based on the MLE of the variance. + * + * @param NU Matrix of residuals + * @return the negative loglikelihood + */ + template + inline + const T dmvlogistic(const DVAR& NU) + { + int n = size_count(NU) - (NU.rowmax()-NU.rowmin()+1); + T var = 1.0/n * norm2(NU); + return n * log(var); + } + + template + inline + const DVAR dmvlogisticResidual(const DATA& O, const DVAR& P) + { + DVAR nu; + nu.allocate(P); + int i; + for( i = O.rowmin(); i<= O.rowmax(); i++) + { + nu(i) = log(O(i)) - log(P(i)); + nu(i) -= mean(nu(i)); + } + return nu; + } + + /** + * @brief Negative loglikelihood for compostion data using multivaritae logistic + * @author Steve Martell + * @param T usually a dmatrix + */ + template + class multivariteLogistic: public negLogLikelihood + { + private: + DATA m_rO; /// ragged observed composition object. + DVAR m_rP; /// ragged predicted composition object. + DVAR m_nu; /// logistic residuals. + + public: + // constructor + // todo: add eps value to constructor. + multivariteLogistic(const DATA &_O, const DVAR &_P, const double eps=1.e-3) + :negLogLikelihood(_O,_P) + { + // tail compression + DATA tmp = this->compress(this->get_O()) + eps; + this->set_rO(tmp); + DVAR vmp = this->compress(this->get_P()) + eps; + this->set_rP(vmp); + + + + // residuals + this->aggregate(eps); + DVAR tnu = acl::dmvlogisticResidual(this->get_rO(),this->get_rP()); + set_nu(tnu); + + } + + + const T nloglike() const + { + return acl::dmvlogistic(this->get_nu()); + } + + const DVAR residual() const + { + return this->get_nu(); + } + + DVAR get_nu() const {return m_nu; } + + void set_nu(DVAR &R){this->m_nu=R;} + + }; + + +} + + + + +#endif \ No newline at end of file diff --git a/src/admb-code/include/lib_iscam.h b/src/admb-code/include/lib_iscam.h index 650b961a..600579d5 100644 --- a/src/admb-code/include/lib_iscam.h +++ b/src/admb-code/include/lib_iscam.h @@ -8,6 +8,7 @@ #include "growthModels.hpp" #include "multinomial.h" #include "selex.hpp" +#include "compositionLikelihoods.hpp" #endif \ No newline at end of file diff --git a/src/admb-code/iscam.tpl b/src/admb-code/iscam.tpl index ad4a550f..04566ac8 100644 --- a/src/admb-code/iscam.tpl +++ b/src/admb-code/iscam.tpl @@ -1233,6 +1233,7 @@ DATA_SECTION init_ivector q_prior(1,nits); init_vector mu_log_q(1,nits); init_vector sd_log_q(1,nits); + init_ivector q_phz(1,nits); // |---------------------------------------------------------------------------------| // | Miscellaneous controls | @@ -1291,12 +1292,13 @@ DATA_SECTION // | ilvec[2] -> number of surveys (nItNobs) // | ilvec[3] -> number of age-compisition data sets (nAgears) // | ilvec[4] -> container for recruitment deviations. - ivector ilvec(1,7); + ivector ilvec(1,8); !! ilvec = ngear; !! ilvec(1) = 1; !! ilvec(2) = nItNobs; - !! ilvec(3) = nAgears; - !! ilvec(4) = ngroup; + !! ilvec(3) = nItNobs; + !! ilvec(4) = nAgears; + !! ilvec(5) = ngroup; // |---------------------------------------------------------------------------------| @@ -1408,29 +1410,16 @@ DATA_SECTION // | MANAGEMENT STRATEGY EVALUATION INPUTS // |---------------------------------------------------------------------------------| // | - - - //LOC_CALCS - // ifstream ifile("Halibut2012.mse"); - // if(ifile) - // { - // cout<<"Vader is happy"< residuals between estimated R and R from S-R curve (process err) // | matrix log_rt(1,n_ag,syr-nage+sage,nyr); - matrix nlvec(1,7,1,ilvec); + matrix nlvec(1,8,1,ilvec); matrix epsilon(1,nItNobs,1,n_it_nobs); matrix xi(1,nItNobs,1,n_it_nobs); matrix it_hat(1,nItNobs,1,n_it_nobs); @@ -1659,7 +1653,7 @@ PARAMETER_SECTION matrix rt(1,ngroup,syr+sage,nyr); matrix delta(1,ngroup,syr+sage,nyr); - + // |---------------------------------------------------------------------------------| // | THREE DIMENSIONAL ARRAYS // |---------------------------------------------------------------------------------| @@ -1721,7 +1715,7 @@ PRELIMINARY_CALCS_SECTION // | - nf is a function evaluation counter. // | - SimFlag comes from the -sim command line argument to simulate fake data. // | - + nf=0; if( testMSY ) { @@ -1870,7 +1864,7 @@ FUNCTION void initParameters() steepness = theta(2); m = mfexp(theta(3)); rho = theta(6); - sigma_r = theta(7); + sigma_r = theta(7); //varphi = sqrt(1.0/theta(7)); sig = elem_prod(sqrt(rho) , varphi); //tau = elem_prod(sqrt(1.0-rho) , varphi); @@ -1960,23 +1954,7 @@ FUNCTION dvector cubic_spline(const dvector& spline_coffs, const dvector& la) //return(1.0*la); } -// FUNCTION dvar_matrix cubic_spline_matrix(const dvar_matrix& spline_coffs) -// { -// RETURN_ARRAYS_INCREMENT(); -// int nodes= spline_coffs.colmax()-spline_coffs.colmin()+1; -// int rmin = spline_coffs.rowmin(); -// int rmax = spline_coffs.rowmax(); - -// dvector ia(1,nodes); -// dvector fa(sage,nage); -// ia.fill_seqadd(0,1./(nodes-1)); -// //fa.fill_seqadd(sage,1); -// fa.fill_seqadd(0,1./(nage-sage)); -// vcubic_spline_function_array fna(rmin,rmax,ia,spline_coffs); -// RETURN_ARRAYS_DECREMENT(); -// return(fna(fa)); - -// } + /** @@ -2580,9 +2558,9 @@ FUNCTION calcComposition dvar_vector pred_ca = ca * age_age(e); A_hat(kk)(ii) = pred_ca(n_A_sage(kk),n_A_nage(kk)); if( n_A_nage(kk) < nage ) - { - A_hat(kk)(ii)(n_A_nage(kk)) += sum( pred_ca(n_A_nage(kk)+1,nage) ); - } + { + A_hat(kk)(ii)(n_A_nage(kk)) += sum( pred_ca(n_A_nage(kk)+1,nage) ); + } } else { @@ -2761,6 +2739,7 @@ FUNCTION calcTotalCatch fishery that would remove potential spawn that would not be surveyed. - d3_survey_data: (iyr index(it) gear area group sex wt timing) - for MLE of survey q, using weighted mean of zt to calculate q. + - March 30, 2015. Added deviation in q for random walk. TODO list: [] - add capability to accomodate priors for survey q's. @@ -2773,11 +2752,12 @@ FUNCTION calcSurveyObservations int ii,kk,ig,nz; double di; - dvariable ftmp; + dvariable ftmp,zbar; dvar_vector Na(sage,nage); dvar_vector va(sage,nage); dvar_vector sa(sage,nage); epsilon.initialize(); + qt.initialize(); xi.initialize(); it_hat.initialize(); @@ -2839,42 +2819,85 @@ FUNCTION calcSurveyObservations dvar_vector t1 = rowsum(V); - dvar_vector zt = log(it) - log(t1(iz,nz)); - dvariable zbar = sum(elem_prod(zt,wt)); - q(kk) = mfexp(zbar); - - // | survey residuals - epsilon(kk).sub(iz,nz) = elem_div(zt - zbar,it_log_se(kk)(iz,nz)); - it_hat(kk).sub(iz,nz) = q(kk) * t1(iz,nz); + dvar_vector zt = log(it); + - // | SPECIAL CASE: penalized random walk in q process error only. - if( q_prior(kk)==2 ) + // | March 30, 2015. Issue #37. + // | Added switch for 3 q_prior(kk) options. + // | q_prior(kk) = 0 = constant fixed Q + // | q_prior(kk) = 1 = constant MLE Q + // | q_prior(kk) = 2 = penalized random walk in Q + switch( q_prior(kk) ) { - epsilon(kk).initialize(); - dvar_vector fd_zt = first_difference(zt); - dvariable zw_bar = sum(elem_prod(fd_zt,wt(iz,nz-1))); - epsilon(kk).sub(iz,nz-1) = elem_div(fd_zt - zw_bar,it_log_se(kk)(iz,nz-1)); - qt(kk)(iz) = exp(zt(iz)); - for(ii=iz+1;ii<=nz;ii++) - { - qt(kk)(ii) = qt(kk)(ii-1) * exp(fd_zt(ii-1)); - } - it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); + case 0: // Constant fixed Q + q(kk) = exp(mu_log_q(kk)); + it_hat(kk).sub(iz,nz) = q(kk) * t1(iz,nz); + zt -= log(it_hat(kk).sub(iz,nz)); + break; + + case 1: // Constant MLE Q + zt -= log(t1(iz,nz)); + zbar = sum(elem_prod(zt,wt)); + q(kk) = mfexp(zbar); + + // | survey residuals + it_hat(kk).sub(iz,nz) = q(kk) * t1(iz,nz); + zt -= zbar; + //epsilon(kk).sub(iz,nz) = elem_div(zt,it_log_se(kk)(iz,nz)); + break; + + case 2: // Penalized random walk in Q + zt -= log(t1(iz,nz)); + qt(kk)(iz) = exp( zt(iz) + log_q_devs(kk)(iz) ); + for(ii=iz+1; ii<=nz; ii++) + { + qt(kk)(ii) = qt(kk)(ii-1) * exp(log_q_devs(kk)(ii)); + } + it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); + zt -= log(qt(kk)(iz,nz)); + //epsilon(kk).sub(iz,nz)= elem_div(zt,it_log_se(kk)(iz,nz)); + break; } - // | MIXED ERROR MODEL for random walk in q - if( q_prior(kk)==3 ) + // Standardized observation error residuals. + epsilon(kk).sub(iz,nz) = elem_div(zt,it_log_se(kk)(iz,nz)); + + // Standardized process error residuals. + if(active(log_q_devs(kk))) { - dvar_vector proerr = zt - zbar; - qt(kk)(ii) = exp(zbar + proerr(iz)); - for(ii=iz+1;ii<=nz;ii++) - { - proerr(ii) = zt(ii) - zt(ii-1); - qt(kk)(ii) = qt(kk)(ii-1) * exp(proerr(ii)); - } - xi(kk).sub(iz,nz) = elem_div(proerr,it_log_pe(kk)(iz,nz)); - it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); + dvar_vector fd_qt = first_difference( log_q_devs(kk)(iz,nz) ); + xi(kk).sub(iz,nz-1) = elem_div(fd_qt,it_log_pe(kk)(iz,nz-1)); } + +// TO BE DEPRECATED +// // | SPECIAL CASE: penalized random walk in q process error only. +// if( q_prior(kk)==2 ) +// { +// epsilon(kk).initialize(); +// dvar_vector fd_zt = first_difference(zt); +// dvariable zw_bar = sum(elem_prod(fd_zt,wt(iz,nz-1))); +// epsilon(kk).sub(iz,nz-1) = elem_div(fd_zt - zw_bar,it_log_se(kk)(iz,nz-1)); +// qt(kk)(iz) = exp(zt(iz)); +// for(ii=iz+1;ii<=nz;ii++) +// { +// qt(kk)(ii) = qt(kk)(ii-1) * exp(fd_zt(ii-1)); +// } +// it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); +// } +// +// // | MIXED ERROR MODEL for random walk in q +// if( q_prior(kk)==3 ) +// { +// dvar_vector proerr = zt - zbar; +// qt(kk)(ii) = exp(zbar + proerr(iz)); +// for(ii=iz+1;ii<=nz;ii++) +// { +// proerr(ii) = zt(ii) - zt(ii-1); +// qt(kk)(ii) = qt(kk)(ii-1) * exp(proerr(ii)); +// } +// xi(kk).sub(iz,nz) = elem_div(proerr,it_log_pe(kk)(iz,nz)); +// it_hat(kk).sub(iz,nz) = elem_prod(qt(kk)(iz,nz),t1(iz,nz)); +// } } if(verbose)cout<<"**** Ok after calcSurveyObservations ****"< *ptr_AgeCompLike; + + + // delete *ptr_AgeComeLike; + + A_nu.initialize(); for(k=1;k<=nAgears;k++) { @@ -3148,6 +3181,10 @@ FUNCTION calcObjectiveFunction //if( iyr <= nyr ) naa++; //if( iyr < syr ) iaa++; } + ptr_AgeCompLike = new acl::multivariteLogistic(O,P,dMinP(k)); + dvariable ell = ptr_AgeCompLike->nloglike(); + //COUT(ell); + //dmatrix O = trans(trans(d3_A_obs(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); //dvar_matrix P = trans(trans(A_hat(k)).sub(n_A_sage(k),n_A_nage(k))).sub(iaa,naa); @@ -3160,19 +3197,20 @@ FUNCTION calcObjectiveFunction switch( int(nCompLikelihood(k)) ) { case 1: - nlvec(3,k) = dmvlogistic(O,P,nu,age_tau2(k),dMinP(k)); + nlvec(4,k) = dmvlogistic(O,P,nu,age_tau2(k),dMinP(k)); + cout<<"like: "<log_sel(k,ig,i,j+1)) - nlvec(6,k)+=lambda_2(k) + nlvec(7,k)+=lambda_2(k) *square( log_sel(k,ig,i,j)-log_sel(k,ig,i,j+1) ); } } @@ -3312,7 +3350,7 @@ FUNCTION calcObjectiveFunction for(j=sage;j<=nage;j++) { dvar_vector df2 = first_difference(first_difference(trans_log_sel(j))); - nlvec(7,k) += lambda_3(k)/(nage-sage+1)*norm2(df2); + nlvec(8,k) += lambda_3(k)/(nage-sage+1)*norm2(df2); } } } @@ -3413,7 +3451,7 @@ FUNCTION calcObjectiveFunction qvec.initialize(); for(k=1;k<=nits;k++) { - if(q_prior(k) == 1 ) + if(q_prior(k) == 1 && sd_log_q(k) != 0 ) { qvec(k) = dnorm( log(q(k)), mu_log_q(k), sd_log_q(k) ); }