Hi
You indeed cannot use conventional optimisation algorithms as the indicator
function is discontinuous, with no derivative. The search hence has to be
done with a grid search over each potential threshold value, which is
termed "conditional/concentrated" least square, or also "profile
likelihood".
Another approach is to make the indicator function smooth, replacing it
with a logistic or cumulative normal, in which case you can then use
conventional methods. This means that instead of having a clear jump from
one regime to the other, you have a smooth transition.
Package tsDyn implements both these models for time series
auto-regressions, the first being called SETAR and the second star. You can
have a look at the code and adapt it easily for your purpose (although
there the y is continuous)!
Regarding litterature, the paper that fits most what you look for is
probably:
-Samia, N. and Chan, K. (2011). Maximum likelihood estimation of a
generalized threshold stochastic regression model. Biometrika,
98(2):433448.
Best
Matthieu
Message: 13
Date: Sat, 20 Oct 2012 09:44:00 -0400
From: V?ronique Boucher Lalonde <veronique.boucher.lalo...@gmail.com
>
To: Ravi Varadhan <ravi.varad...@jhu.edu>
Cc: "r-help@r-project.org" <r-help@r-project.org>, Bert Gunter
<gunter.ber...@gene.com>
Subject: Re: [R] fit a "threshold" function with nls
Message-ID:
<cabkqxbzspkfe6yqve85phkdp8ando7ri6gcmzk8mafrs2cb...@mail.gmail.com>
Content-Type: text/plain
Thank you all for your help. I think I now understand the issue.
I tried to write a likelihood function for my binomial model.
Please excuse my ignorance if I am not doing this right; I do not have any
statistical background.
#Example data
x <- seq(0, 1000)
y <- ifelse(x < 300, 0, ifelse(x>700, 0, 1))
#Function assuming binomial errors
fn <- function(par) {
y.pred = ifelse(x<par[1], 0, ifelse(x>par[2], 0, 1))
-sum(dbinom(1-abs(y-y.pred), prob = mean(y.pred == y), size = 1, log = T))
}
optim(par=c(300,700), fn)
Would the "Nelder-Mead" optimization method be appropriate?
On Tue, Oct 16, 2012 at 11:28 AM, Ravi Varadhan <ravi.varad...@jhu.edu
>wrote:
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
