Skip to content

bug found in predict.locfit in locfit package (PR#8057)

4 messages · apipatta@colorado.edu, Erich Neuwirth, Martin Maechler

Original
apipatta@colorado.edu
# 
Full_Name: Somkiat Apipattanavis
Version: 2.1.1
OS: Windows
Submission from: (NULL) (128.138.44.123)


Bug found in predict.locfit for density estimation

# Example of bug found in prdict.locfit (Locfit)
library('locfit')

# generate data
y =c(4281,2497,4346,5588,5593,3474,4291,2542,5195,4056,
     3114,2864,4904,7625,3377,4001,4999,7191,8062,5668)
x1=c( 0.258729, 1.460156, 0.192323, 0.067062,-0.304291,
     -0.420917, 0.214729, 0.239979,-0.421938,-0.571229,
      1.310990, 2.043032, 0.449906,-0.951917,-0.077104,
     -0.356833,-0.286042, 0.065750, 0.159677,-0.075792)
x2=c(-0.3050, 1.0125, 0.2050, 0.1025, 0.9550,
      0.6975, 1.5550, 0.0225, 0.2575, 0.3725,
      2.0075, 2.1275, 0.7200, 0.2950, 0.2875,
     -0.2800,-0.6050, 0.2125,-0.5525,-1.7850)

ndat=length(y)
ybk=y
x1bk=x1
x2bk=x2
######## Joint probability function of y, x1 and x2
# fit joint probability function
fityxv=locfit(~y+x1+x2,alpha=1,deg=1)
fyxv=predict(fityxv,where="data")

######## Marginal distribution of gxv
# fit marginal distribution of y
fitxv=locfit(~x1+x2,alpha=0.5,deg=1)
gxv=predict(fitxv,where="data")

######## Prediction of fyxv and gxv
# new data
vx1=0.2
vx2=0.7
x1new=rep(vx1,ndat)
x2new=rep(vx2,ndat)
ynew=y

# marginal distribution of gxv for new data
newdata=data.frame(x1new,x2new)
gxvnew=predict(fitxv,newdata)   #bug!!! gave the same values as gxv

# This bug can be avoid by setting new values into old variables
# then, we will get the new predicted values
# for example
x1=x1new
x2=x2new
gxvnew2=predict(fitxv,where="data")

# predict joint probability function of fyxv for new data
newdat2=data.frame(ynew,x1new,x2new)
fyxvnew=predict(fityxv,newdat2) #bug! same as 2D density

# but setting new values into old variables DOES NOT
# work for solving this kind of problem for 3D density
# for example
x1=x1bk
x2=x2bk
fyxvnew2=predict(fityxv,where="data")
12 days later
Erich Neuwirth
# 
In R 2.2.0 density now can work with weighted obesrvations.
It would be nice if boxplot also would
accept a weight parameter,
then one could produce consistent density estimators and boxplots.

Could the developers consider adding this feature?
1 day later
Martin Maechler
# 
Erich> In R 2.2.0 density now can work with weighted
    Erich> obesrvations.  It would be nice if boxplot also would
    Erich> accept a weight parameter, then one could produce
    Erich> consistent density estimators and boxplots.

    Erich> Could the developers consider adding this feature?

The first thing I'd want is  quantile() with weights --- which I
personally find quite interesting and have wanted several times
in the past --- not wanted enough to implement though.

I'm interested to hear of (or even see C or R implementations of)
fast algorithms for "weight quantiles".  
Code contributions are welcome too..

(And yes, I do know that boxplots are base on "hinges" rather than
 quartiles but that's less interesting here.)

Martin Maechler <maechler at stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: +41-44-632-3408		fax: ...-1228			<><
2 days later
Erich Neuwirth
# 
Martin,

Frank Harrell's Hmisc has weigthed variants of quite a few statistical
estimators, including quantiles. I never have look at
how efficiently this is implemented, but as far as I know
it works.

Erich
Martin Maechler wrote: