Dear Colleagues, I have x, y data (pollen and seed dispersal from oaks !) that I would like to fit with a function which look like this: p(a,b,x,y)=b/(2*pi*a?gamma(2/b))*exp(-(square_root(x?+y?)/a)power(b)) I am looking for a and b values that fit my data at best. Can someone give me hints to perform such an analysis with R ? Thanks a lot Sophie Sophie Gerber gerber at pierroton.inra.fr INRA - UMR BIOGECO 69 route d'Arcachon tel (33) (0)5 57 12 28 30 33612 Cestas cedex fax (33) (0)5 57 12 28 81 http://www.pierroton.inra.fr/genetics/Perso/Sophie/ page_sophie_english.html
Have you considered "nls" or "optim"? My way of handling this kind of problem is to assume that x and / or y follow some probability distribution with parameters a and b. Then I write a function to compute deviance = (-2*log(likelihood)) = (-2*log(probability density for x and / or y given a and b)) and feed it to "optim" to minimize this "deviance". Or if the response variable is a nonlinear model plus independent normal errors with constant variance, I may use "nls". hope this helps. spencer graves Sophie Gerber wrote:> Dear Colleagues, > > I have x, y data (pollen and seed dispersal from oaks !) that I would > like to fit with a function which look like this: > > p(a,b,x,y)=b/(2*pi*a?gamma(2/b))*exp(-(square_root(x?+y?)/a)power(b)) > > I am looking for a and b values that fit my data at best. > Can someone give me hints to perform such an analysis with R ? > > Thanks a lot > > Sophie > > > Sophie Gerber gerber at pierroton.inra.fr > INRA - UMR BIOGECO > 69 route d'Arcachon tel (33) (0)5 57 12 28 30 > 33612 Cestas cedex fax (33) (0)5 57 12 28 81 > http://www.pierroton.inra.fr/genetics/Perso/Sophie/ > page_sophie_english.html > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Your function is not quite complete -- you need to specify an error distribution such as Poisson or negative binomial (see literature from the past decade by J. Clark, Muller-Landau, Ribbens and others) Assuming Poisson errors, if x and y are locations and d is observed density: f <- function(p) { a <- p[1] b <- p[2] r <- sqrt(x^2+y^2) exp.dens = b/(2*pi*a^2*gamma(2/b))*exp(-(r/a)^b) lik = -sum(dpois(d,exp.dens,log=TRUE)) } optim(fn=f,par=(a vector of reasonable starting values)) I see from my e-mail that Spencer Graves has just replied to this message too so I'll go see what he wrote ... On Wed, 22 Oct 2003, Sophie Gerber wrote:> Dear Colleagues, > > I have x, y data (pollen and seed dispersal from oaks !) that I would > like to fit with a function which look like this: > > p(a,b,x,y)=b/(2*pi*a?gamma(2/b))*exp(-(square_root(x?+y?)/a)power(b)) > > I am looking for a and b values that fit my data at best. > Can someone give me hints to perform such an analysis with R ? > > Thanks a lot > > Sophie > > > Sophie Gerber gerber at pierroton.inra.fr > INRA - UMR BIOGECO > 69 route d'Arcachon tel (33) (0)5 57 12 28 30 > 33612 Cestas cedex fax (33) (0)5 57 12 28 81 > http://www.pierroton.inra.fr/genetics/Perso/Sophie/ > page_sophie_english.html > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help >-- 620B Bartram Hall bolker at zoo.ufl.edu Zoology Department, University of Florida http://www.zoo.ufl.edu/bolker Box 118525 (ph) 352-392-5697 Gainesville, FL 32611-8525 (fax) 352-392-3704