> Date: Mon, 20 Jun 2011 11:25:54 +0200
> From: lig...@statistik.tu-dortmund.de
> To: p.e.b...@dunelm.org.uk
> CC: r-help@r-project.org; pb...@astro.uni-bonn.de
> Subject: Re: [R] different results from nls in 2.10.1 and 2.11.1
>
> Since one is a 32-bit and the other one a 64-bit, and therefore the
> compiler is also different, you can get different numerical results
> easily, even with identical versions of R (and most of us do not have
> outdated R installations around).
>
> I just tried your example on 3 different systems with R-2.13.0 and all
> told me "singular convergence"...
>
> Uwe Ligges
>
>
>
>
>
> On 18.06.2011 15:44, Philip Bett wrote:
> > Hi,
> > I've noticed I get different results fitting a function to some data on
> > my laptop to when I do it on my computer at work.
I guess the best all around approach is to dump the results from your
FisherAvgPdf function and get some idea what trajectory it takes
in the different cases. This is presumably not just an issue with R
and could have something to tell you about your data vs the model.
Even if it converged without incident, you'd probably want to
look at more details. You apparently are trying to fit a histogram to
a pdf and it is not too hard to just plot fit over histo and
spot check parameter space.
Consider for example,
plot(h$mids,h$density,type="b")
points(h$mids,FisherAvgdPdf(acos(h$mids),3.527e-8,3.198),pch=20)
points(h$mids,FisherAvgdPdf(acos(h$mids),3.527e-8,3.298),pch=13)
points(h$mids,FisherAvgdPdf(acos(h$mids),4.527e-8,3.198),pch=14)
and the find various measures for how good/bad these are etc.
e<-h$density-FisherAvgdPdf(acos(h$mids),3.527e-8,3.198)
sum(e*e)
e<-h$density-FisherAvgdPdf(acos(h$mids),3.527e-8,3.298)
sum(e*e)
the error surface as function of parameters seems smooth etc,
foo <- function (x,y) { e<-h$density-FisherAvgdPdf(acos(h$mids),x,y) ; sum(e*e)
}
x=(2+(1:100)*.02)/1e8;
y=3+((1:100)/100);
df<-data.frame(x=0,y=0,z=0);
for ( i in x ) { for ( j in y ) { gg<-data.frame(i,j,foo(i,j));str(gg);
colnames(gg)<-colnames(df); df<-rbind(df,gg); }}
sel=which(df$x!=0)
scatterplot3d(df$x[sel],df$y[sel],df$z[sel])
> >
> > Here's a code snippet of what I do:
> >
> > ##------------------------------------------------------------------
> > require(circular) ## for Bessel function I.0
> >
> > ## Data:
> > dd <- c(0.9975948929787, 0.9093316197395, 0.7838819026947,
> > 0.9096108675003, 0.8901804089546, 0.2995955049992, 0.9461286067963,
> > 0.8248071670532, 0.2442084848881, 0.2836948633194, 0.7353935241699,
> > 0.5812761187553, 0.8705610632896, 0.8744471669197, 0.7490273118019,
> > 0.9947383403778, 0.9154829382896, 0.8659985661507, 0.6448246836662,
> > 0.8588128685951, 0.7347437739372, -0.1645197421312, 0.970999121666,
> > 0.8038327097893, 0.9558997154236, 0.6846113204956, 0.6286814808846,
> > 0.9201356172562, 0.9422197341919, 0.3470877110958, 0.4154576957226,
> > 0.0721184238791, 0.14151956141, -0.6142936348915, -0.4688512086868,
> > 0.6805665493011, 0.3594025671482, 0.8991097211838, 0.7656877636909,
> > 0.9282909035683, 0.9454715847969, 0.9766132831573, 0.4316343963146,
> > 0.62679708004, 0.2093886137009, 0.3937581181526, 0.4254160523415,
> > 0.8684504628181, 0.3844584524632, 0.9578431844711, 0.956972181797,
> > 0.4456568360329, 0.9793710708618, 0.5825698971748, 0.929228246212,
> > 0.9211971759796, 0.9407976865768, 0.821156680584, 0.2048042863607,
> > 0.6473184227943, 0.9456319212914, 0.7021154165268, 0.9761978387833,
> > 0.1485801786184, 0.2195029109716, 0.5378785729408, 0.8304615020752,
> > 0.8596342802048, 0.950027525425, 0.9102076888084, 0.5108731985092,
> > 0.7200184464455, 0.3571084141731, 0.9765330553055, -0.143017962575,
> > 0.8576183915138, 0.1283493340015, -0.3226098418236, 0.7031792402267,
> > 0.8708582520485, 0.56754809618, 0.060470353812, 0.8015220761299,
> > 0.7363410592079, 0.671902179718, 0.8082517385483, 0.9468197822571,
> > 0.9729647636414, 0.7919752597809, 0.9539568424225, 0.4840737581253,
> > 0.850653231144, 0.5909016132355, 0.8414449691772, 0.9699150323868)
> >
> > xlims <- c(-1,1)
> > bw <- 0.05
> > b <- seq(xlims[1],xlims[2],by=bw) ; nb <- length(b)
> > h <- hist( dd, breaks=b, plot=FALSE)
> >
> > FisherAvgdPdf <- function(theta,theta0,kappa){
> > A <- kappa/(2*sinh(kappa))
> > A * I.0( kappa*sin(theta)*sin(theta0) ) * exp(
> > kappa*cos(theta)*cos(theta0) )
> > }
> >
> >
> > nls(dens ~ FisherAvgdPdf(theta,theta0,kappa),
> > data = data.frame( theta=acos(h$mids), dens=h$density ),
> > start=c( theta0=0.5, kappa=4.0 ),
> > algorithm="port", lower=c(0,0), upper=c(acos(xlims[1]),500),
> > control=list(warnOnly=TRUE) )
> >
> > ##------------------------------------------------------------------
> >
> > On one machine, nls converges, and on the other it doesn't. Any ideas
> > why, and which is "right"? I can't see anything in R News that could be
> > relevant.
> >
> >
> > The different R versions (and computers) are:
> > > R.version
> > _
> > platform i686-pc-linux-gnu
> > arch i686
> > os linux-gnu
> > system i686, linux-gnu
> > status
> > major 2
> > minor 11.1
> > year 2010
> > month 05
> > day 31
> > svn rev 52157
> > language R
> > version.string R version 2.11.1 (2010-05-31)
> >
> >
> > and
> >
> > > R.version
> > _
> > platform x86_64-pc-linux-gnu
> > arch x86_64
> > os linux-gnu
> > system x86_64, linux-gnu
> > status
> > major 2
> > minor 10.1
> > year 2009
> > month 12
> > day 14
> > svn rev 50720
> > language R
> > version.string R version 2.10.1 (2009-12-14)
> >
> >
> >
> >
> > Thanks,
> >
> > Phil Bett
> >
> >
> >
>
> ______________________________________________
> 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.
[[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.
