Nothing
#' Local polynomial estimator of a psychometric function
#'
#' Local polynomial estimator for the psychometric function and eta function (psychometric function
#' transformed by link) for binomial data; also returns the hat matrix H.
#'
#' @usage locglmfit( xfit, r, m, x, h, returnH = FALSE, link = "logit",
#' guessing = 0, lapsing = 0, K = 2, p = 1,
#' ker = "dnorm", maxiter = 50, tol = 1e-6 )
#'
#' @param xfit points at which to calculate the estimate pfit
#' @param r number of successes at points x
#' @param m number of trials at points x
#' @param x stimulus levels
#' @param h bandwidth(s)
#' @param returnH (optional) logical, if TRUE then hat matrix is calculated; default is FALSE
#' @param link (optional) name of the link function; default is 'logit'
#' @param guessing (optional) guessing rate; default is 0
#' @param lapsing (optional) lapsing rate; default is 0
#' @param K (optional) power parameter for Weibull and reverse Weibull link; default is 2
#' @param p (optional) degree of the polynomial; default is 1
#' @param ker (optional) kernel function for weights; default is 'dnorm'
#' @param maxiter (optional) maximum number of iterations in Fisher scoring; default is 50
#' @param tol (optional) tolerance level at which to stop Fisher scoring; default is 1e-6
#'
#' @returns \verb{pfit } value of the local polynomial estimate at points xfit
#' @returns \verb{etafit } estimate of eta (link of pfit)
#' @returns \verb{H } hat matrix (OPTIONAL)
#'
#' @examples
#' data("Miranda_Henson")
#' x = Miranda_Henson$x
#' r = Miranda_Henson$r
#' m = Miranda_Henson$m
#' numxfit <- 199; # Number of new points to be generated minus 1
#' xfit <- (max(x)-min(x)) * (0:numxfit) / numxfit + min(x)
#' # Find a plug-in bandwidth
#' bwd <- bandwidth_plugin( r, m, x)
#' pfit <- locglmfit( xfit, r, m, x, bwd )$pfit
#' # Plot the fitted curve
#' plot( x, r / m, xlim = c( 0.1, 1.302 ), ylim = c( 0.0165, 0.965 ), type = "p", pch="*" )
#' lines(xfit, pfit )
#'
#' @export
locglmfit<-function( xfit, r, m, x, h, returnH = FALSE, link = "logit",
guessing = 0, lapsing = 0, K = 2, p = 1,
ker = "dnorm", maxiter = 50, tol = 1e-6 ) {
#
# Local polynomial estimator for the psychometric function and eta function (psychometric function
# transformed by link) for binomial data; also returns the hat matrix H. Actual calculations are
# done in LOCGLMFIT_PRIVATE or LOCGLMFIT_SPARSE_PRIVATE depending on the size of the data set.
# Here, the data are split into several parts to speed up the calculations.
#
#
#INPUT
#
# xfit - points at which to calculate the estimate pfit
# r - number of successes at points x
# m - number of trials at points x
# x - stimulus levels
# h - bandwidth(s)
#
# OPTIONAL INPUT
#
# returnH - logical, if TRUE then hat matrix is calculated; default is FALSE
# link - name of the link function; default is 'logit'
# guessing - guessing rate; default is 0
# lapsing - lapsing rate; default is 0
# K - power parameter for Weibull and reverse Weibull link; default is 2
# p - degree of the polynomial; default is 1
# ker - kernel function for weights; default is 'dnorm'
# maxiter - maximum number of iterations in Fisher scoring; default is 50
# tol - tolerance level at which to stop Fisher scoring; default is 1e-6
#
# OUTPUT
#
# pfit - value of the local polynomial estimate at points xfit
# etafit - estimate of eta (link of pfit)
# H - hat matrix (OPTIONAL)
# First 5 arguments are mandatory
if( missing("xfit") || missing("r") || missing("m") || missing("x") || missing("h") ) {
stop("Check input. First 5 arguments are mandatory");
}
# CHECK ROBUSTNESS OF INPUT PARAMETERS
if ( !is.vector( xfit ) ) {
stop("Vector xfit should be a column vector");
}
checkdata<-list();
checkdata[[1]] <- x;
checkdata[[2]] <- r;
checkdata[[3]] <- m;
checkinput("checkdesignpoints",x)
checkinput( "psychometricdata", checkdata );
rm( checkdata )
if ( !is.vector( h ) ) {
stop("Bandwidths h should be a vector");
}
if( !( length( h ) == 1 || length( h ) == length( xfit ) ) ) {
stop( "Bandwidth h must be either a scalar or a vector with the same number of elements as xfit" );
}
checkinput( "linkfunction", link );
if( length( guessing ) > 1 ) {
stop( "Guessing rate must be a scalar" );
}
if( length( lapsing ) > 1 ) {
stop( "Lapsing rate must be a scalar" );
}
checkinput( 'guessingandlapsing', c(guessing, lapsing) );
if (link == "weibull" || link == "revweibull"){
checkinput( "exponentk", K );
}
pn <- list()
pn[[1]] <- p
pn[[2]] <- x
checkinput( "degreepolynomial", pn );
checkinput( 'kernel', ker );
checkinput( 'maxiter', maxiter );
checkinput( 'tolerance', tol );
# INITIAL VALUES
split <- 20;
Lxfit <- length(xfit);
Lx <- length(x);
value <- NULL
pfit <- NULL;
etafit <- NULL;
if( returnH ) H <- NULL;
if( link == "logit" ||
link == "probit" ||
link == "loglog" ||
link == "comploglog" ||
link == "weibull" ||
link == "revweibull" ) {
link <- paste( link, "_link_private", sep = "" );
}
if( Lx > 15 ) {
# big data
# First try to load package SparseM
if (requireNamespace("SparseM", quietly = TRUE)) {
fun_estim <- locglmfit_sparse_private;
} else {
fun_estim <- locglmfit_private;
message("Package SparseM not installed. No sparse matrices used");
message("SparseM can be found at CRAN web site http://cran.r-project.org/");
}
# options(warn = -1);
# existSparseM <- library( SparseM, logical.return = TRUE );
# options(warn = 0);
# if( existSparseM ) {
# fun_estim <- locglmfit_sparse_private;
# }
# else {
# fun_estim <- locglmfit_private;
# message("Package SparseM not installed. No sparse matrices used");
# message("SparseM can be found at CRAN web site http://cran.r-project.org/");
#}
}
else {
# small data
fun_estim <- locglmfit_private;
}
# SPLIT AND EVALUATION
################################################################ SCALAR h
if( length( h ) == 1 ) {
# with Hat matrix
if( returnH ) {
if( Lxfit <= split ) {
# small x
value <- fun_estim( xfit, r, m, x, h, returnH, link, guessing,
lapsing, K, p, ker, maxiter, tol);
}
else {
# large x
# number of parts into which the fitting is divided
fLx = floor( Lxfit / split );
# initialise output
for( i in 0:(fLx-1) ) {
# part of the fit
value1 <- fun_estim( xfit[i*split+c(1:split)], r, m, x, h,
returnH, link, guessing, lapsing, K, p, ker,
maxiter, tol );
# put the fits together
pfit <- c( pfit, value1$pfit );
etafit <- c( etafit, value1$etafit );
H <- rbind(H, value1$H );
}
# final part of the fit
if( ( split * fLx ) < Lxfit ) {
value1 <- fun_estim( xfit[c(1+split*fLx):Lxfit], r, m, x, h,
returnH, link, guessing, lapsing, K, p, ker,
maxiter, tol );
# put the fits together
pfit <- c( pfit, value1$pfit );
etafit <- c( etafit, value1$etafit );
H <- rbind(H, value1$H );
}
# values to return
value$pfit <- pfit
value$etafit <- etafit
value$H <- H
}
}
else { # no Hat matrix
if( Lxfit <= split ) {
# small x
value <- fun_estim( xfit, r, m, x, h,returnH, link, guessing,
lapsing, K, p, ker, maxiter, tol );
}
else {
# large x
# number of parts into which the fitting is divided
fLx = floor( Lxfit / split );
# initialise output
for( i in 0:(fLx-1) ) {
# part of the fit
value1 <- fun_estim( xfit[i*split + c(1:split)], r, m, x, h,
returnH, link, guessing, lapsing, K, p, ker,
maxiter, tol );
# put the fits together
pfit <- c( pfit, value1$pfit );
etafit <- c( etafit, value1$etafit );
}
# final part of the fit
if( ( split * fLx ) < Lxfit ) {
value1 <- fun_estim( xfit[c(1+split*fLx):Lxfit], r, m, x, h,
returnH, link, guessing, lapsing, K, p, ker,
maxiter, tol );
# put the fits together
pfit <- c( pfit, value1$pfit );
etafit <- c( etafit, value1$etafit );
}
# values to return
value$pfit <- pfit
value$etafit <- etafit
}
} # if( returnH )
}
################################################################ VECTOR h
else { # if( length( h ) == 1 )
# with Hat matrix
if( returnH ) {
if ( Lxfit <= split ) {
# small x
value <- fun_estim( xfit, r, m, x, h, returnH, link, guessing,
lapsing, K, p, ker, maxiter, tol );
}
else {
# large x
# number of parts into which the fitting is divided
fLx = floor( Lxfit / split );
for( i in 0:(fLx-1) ) {
# part of the fit
value1 <- fun_estim( xfit[i*split + c(1:split)], r, m, x,
h[i*split + c(1:split)], returnH, link,
guessing, lapsing, K, p, ker, maxiter,
tol );
# put the fits together
pfit <- c( pfit, value1$pfit );
etafit <- c( etafit, value1$etafit );
H <- rbind(H, value1$H );
}
# final part of the fit
if( ( split * fLx ) < Lxfit ) {
value1 <- fun_estim( xfit[c(1+split*fLx):Lxfit], r, m, x,
h[c(1+split*fLx):Lxfit], returnH, link,
guessing, lapsing, K, p, ker, maxiter,
tol );
# put the fits together
pfit <- c( pfit, value1$pfit );
etafit <- c( etafit, value1$etafit );
H <- rbind(H, value1$H );
}
# values to return
value$pfit <- pfit
value$etafit <- etafit
value$H <- H
}
}# end if(retunH)
else { # no Hat matrix
if( Lxfit <= split ) {
# small x
value <- fun_estim( xfit, r, m, x, h, returnH, link, guessing,
lapsing, K, p, ker, maxiter, tol );
}
else {
# large x
# number of parts into which the fitting is divided
fLx = floor( Lxfit / split );
# initialise output
for( i in 0:(fLx-1) ) {
# part of the fit
value1 <- fun_estim( xfit[i*split + c(1:split)], r, m, x,
h[i*split + c(1:split)], returnH, link,
guessing, lapsing, K, p, ker, maxiter, tol );
# put the fits together
pfit <- c( pfit, value1$pfit );
etafit <- c( etafit, value1$etafit );
}
# final part of the fit
if( ( split * fLx ) < Lxfit ) {
value1 <- fun_estim( xfit[c(1+split*fLx):Lxfit], r, m, x,
h[c(1+split*fLx):Lxfit], returnH, link,
guessing, lapsing, K, p, ker, maxiter, tol );
# put the fits together
pfit <- c( pfit, value1$pfit );
etafit <- c( etafit, value1$etafit );
}
# values to return
value$pfit <- pfit
value$etafit <- etafit
}
} # if( returnH )
} # if( length( h ) == 1 )
return( value )
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.