[R] Error "singular gradient matrix at initial parameter estimates" in nls
Greetings,
I am struggling a bit with a non-linear regression. The problem is
described below with the known values r and D inidcated.
I tried to alter the start values but get always following error
message:
Error in nlsModel(formula, mf, start, wts):
singular gradient matrix at initial parameter estimates
Calls: nls -> switch -> nlsModel
I might be missing something with regard to the different algorithms or
I just have to try again with some different start values. If anyone
finds some time to have a look at that, any advice would be very much
appreciated.
##Brice-Model
bricemod<-nls(log10(D) ~ log10(Do)*exp(((-4*pi*En*Na)*((ro/2)*(r-ro)^2
+(1/3)*(r-ro)^3))/(R*T)), start=list(Do=0.8, En=390,
ro=1.03),trace=TRUE)
Na<-6.0221415*10^23 # Avrogadro's number
T<-1010 # Temp in K
R<-8.3144 #Gas constant [J mol^-1 K^-1]
r<-data.matrix(PRr[c("r")])
r
La 1.160
Ce 1.143
Pr 1.109
Nd 1.093
Sm1.079
Eu 1.066
Gd 1.053
Tb 1.040
Dy 1.027
Ho 1.015
Er 1.004
Tm 0.994
Yb 0.985
Lu 0.977
D<-data.matrix(PRr[c("D")])
D
La 0.1806551
Ce 0.2703113
Pr 0.3757225
Nd 0.5271811
Sm 0.8665835
Eu 1.0812568
Gd 1.0612762
Tb 1.0726612
Dy 1.1679270
Ho 1.1910920
Er 1.1336938
Tm 1.1215107
Yb 0.9619603
Lu 0.8315467
Niklaus
--
Niklaus Hürlimann
Doctorant-(cand.PhD)
Université de Lausanne
Institut de Minéralogie et Géochimie
L'Anthropole
CH-1015 Lausanne
Suisse
E-mail: niklaus.hurlim...@unil.ch
Tel:+41(0)21 692 4452
[[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.
Re: [R] Fitting a half-ellipse curve
Hello Michael,
Thanks very much for the hint to my problem this helps really a lot.
Cheers
Niklaus
--
Niklaus Hürlimann
Doctorant-PhD
Université de Lausanne
Institut de Minéralogie et Géochimie
L'Anthropole
CH-1015 Lausanne
Suisse
E-mail: niklaus.hurlim...@unil.ch
Tel:+41(0)21 692 4452
--- Begin Message ---
Hello Niklaus,
I'm not sure if the following is the sort of thing you are looking for (?)
You can fit an ellipse to your data using a deterministic least
squares method. The following is a function that I use to do this...
fit.ellipse <- function (x, y = NULL)
{
# Least squares fitting of an ellipse to point data
# using the algorithm described in:
# Radim Halir & Jan Flusser. 1998.
# Numerically stable direct least squares fitting of ellipses.
# Proceedings of the 6th International Conference in Central Europe
# on Computer Graphics and Visualization. WSCG '98, p. 125-132
#
# Adapted from the original Matlab code by Michael Bedward
#
#
# Arguments:
# x, y - the x and y coordinates of the data points
#
# Returns a list with the following elements:
#
# coef - coefficients of the ellipse as described by the general
#quadratic: ax^2 + bxy + cy^2 + dx + ey + f = 0
#
# center - center x and y
#
# major - major semi-axis length
#
# minor - minor semi-axis length
#
EPS <- 1.0e-8
dat <- xy.coords(x, y)
D1 <- cbind(dat$x * dat$x, dat$x * dat$y, dat$y * dat$y)
D2 <- cbind(dat$x, dat$y, 1)
S1 <- t(D1) %*% D1
S2 <- t(D1) %*% D2
S3 <- t(D2) %*% D2
T <- -solve(S3) %*% t(S2)
M <- S1 + S2 %*% T
M <- rbind(M[3,] / 2, -M[2,], M[1,] / 2)
evec <- eigen(M)$vec
cond <- 4 * evec[1,] * evec[3,] - evec[2,]^2
a1 <- evec[, which(cond > 0)]
f <- c(a1, T %*% a1)
names(f) <- letters[1:6]
# calculate the center and lengths of the semi-axes
b2 <- f[2]^2 / 4
center <- c((f[3] * f[4] / 2 - b2 * f[5]), (f[1] * f[5] / 2 - f[2] *
f[4] / 4)) / (b2 - f[1] * f[3])
names(center) <- c("x", "y")
num <- 2 * (f[1] * f[5]^2 / 4 + f[3] * f[4]^2 / 4 + f[6] * b2 -
f[2]*f[4]*f[5]/4 - f[1]*f[3]*f[6])
den1 <- (b2 - f[1]*f[3])
den2 <- sqrt((f[1] - f[3])^2 + 4*b2)
den3 <- f[1] + f[3]
semi.axes <- sqrt(c( num / (den1 * (den2 - den3)), num / (den1 *
(-den2 - den3)) ))
# calculate the angle of rotation
term <- (f[1] - f[3]) / f[2]
angle <- atan(1 / term) / 2
list(coef=f, center = center, major = max(semi.axes), minor =
min(semi.axes), angle = unname(angle))
}
There are quite a few functions / packages to draw ellipses in R, but
the following is will work directly with the output of the above
function...
get.ellipse <- function ( fit, n=360 )
{
# Calculate points on an ellipse described by
# the fit argument as returned by fit.ellipse
#
# n is the number of points to render
tt <- seq(0, 2*pi, length=n)
sa <- sin(fit$angle)
ca <- cos(fit$angle)
ct <- cos(tt)
st <- sin(tt)
x <- fit$center[1] + fit$maj * ct * ca - fit$min * st * sa
y <- fit$center[2] + fit$maj * ct * sa + fit$min * st * ca
cbind(x=x, y=y)
}
So if your data were in a matrix or data.frame X...
efit <- fit.ellipse( X )
e <- get.ellipse( efit )
plot(X)
lines(e, col="red")
Hope this helps,
Michael
On 29 September 2010 23:45, Niklaus Hurlimann wrote:
> Dear mailing list,
>
> I have following array:
>
> X2 Y2
> [1,] 422.7900 6.0
> [2,] 469.8007 10.5
> [3,] 483.9428 11.0
> [4,] 532.4917 25.5
> [5,] 596.1942 33.5
> [6,] 630.8496 40.5
> [7,] 733.2996 45.0
> [8,] 946.4779 32.0
> [9,] 996.8068 35.5
> [10,] 1074.3310 23.0
>
> I do afterwards the following:
>
> plot.new()
>
> plot.window(xlim=c(min(X1)-50,max(X1)+50),
> ylim=c(min(Y1)-25,max(Y1)+25))
>
> axis(2, cex.axis=1.2)
> axis(1, cex.axis=1.2)
>
> points(X2, Y2, type="p", pch=0, cex=1.2, col="black")
>
> title(main="Dyke Geometry Along Strike", cex.main=1.2, font.main=4)
> title(xlab="distance [m]", cex.lab=1.2)
> title(ylab="half-thickness [m]", cex.lab=1.2)
>
> box()
>
>
> I would like curve fitting where I proceeded with a non
> linear-regression with the according function below:
>
> nls(formula = Y2 ~ abs(b*abs(1-X2^2/a^2)^(1/2)), start = list( a=282,
> b=40))
>
> The formula should give the y-positive part of an ellipse in my opinion
> or I might be completely wrong.
>
> In the end I would like to fit a curve of half an ellipse through the
> data. I might could do this as well with a 2nd order polynomial fit. I
> am grateful for any suggestions and comments to my problem.
>
>
> Cheers
>
>
>
> --
> Niklaus Hürlimann
> Doctorant-PhD
>
> Université de Lausanne
> Instit
[R] Fitting a half-ellipse curve
Dear mailing list, I have following array: X2 Y2 [1,] 422.7900 6.0 [2,] 469.8007 10.5 [3,] 483.9428 11.0 [4,] 532.4917 25.5 [5,] 596.1942 33.5 [6,] 630.8496 40.5 [7,] 733.2996 45.0 [8,] 946.4779 32.0 [9,] 996.8068 35.5 [10,] 1074.3310 23.0 I do afterwards the following: plot.new() plot.window(xlim=c(min(X1)-50,max(X1)+50), ylim=c(min(Y1)-25,max(Y1)+25)) axis(2, cex.axis=1.2) axis(1, cex.axis=1.2) points(X2, Y2, type="p", pch=0, cex=1.2, col="black") title(main="Dyke Geometry Along Strike", cex.main=1.2, font.main=4) title(xlab="distance [m]", cex.lab=1.2) title(ylab="half-thickness [m]", cex.lab=1.2) box() I would like curve fitting where I proceeded with a non linear-regression with the according function below: nls(formula = Y2 ~ abs(b*abs(1-X2^2/a^2)^(1/2)), start = list( a=282, b=40)) The formula should give the y-positive part of an ellipse in my opinion or I might be completely wrong. In the end I would like to fit a curve of half an ellipse through the data. I might could do this as well with a 2nd order polynomial fit. I am grateful for any suggestions and comments to my problem. Cheers -- Niklaus Hürlimann Doctorant-PhD Université de Lausanne Institut de Minéralogie et Géochimie L'Anthropole CH-1015 Lausanne Suisse E-mail: niklaus.hurlim...@unil.ch Tel:+41(0)21 692 4452 [[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.
[R] plotting density in same plot in loop iteration
Hi R-mailing list
I would have the following set-up below with a simplified data-frame.
Through a loop which includes certain criteria for the densities I would
like to plot the different density-distributions in the same plot. Of
course I hope I don't do any mistakes with all the indexes of the
dataframe.
All I would like to have is the different densities in the same plot
with a general background density d.
In a second step there might be also the need to combine all the
different densities within the loop to one single distribution.
I would be grateful for any help on this.
Cheers
> a
[1] 1 4 7 11 16 21 23
> b
[1] 15 17 13 23 14 11 12
> c
[1] 0 15 0 0 0 0 0
> d
[1] 0 0 7 0 16 0 0
> e
[1] 0 0 0 11 0 21 0
> f
[1] 0 0 123 124 122 121 0
q<- data.frame(a,b,c,d,e,f)
a b c d e f
1 1 15 0 0 0 0
2 4 17 15 0 0 0
3 7 13 0 7 0 123
4 11 23 0 0 11 124
5 16 14 0 16 0 122
6 21 11 0 0 21 121
7 23 12 0 0 0 0
indx<-which(q[,"a"]==q[,"d"])
indy<-which(q[,"a"]==q[,"e"])
s<-[c(indx,indy),]
std<-sd(S[,"f"])
t<-round(max(q[,"a"]),0)
u<-round(min(q[,"a"]),0)
colnames(v)<-c("a")
D<-q[,"a"]
d<-density(D)
pdf('L')
for (h in length(s[["a"]])) {
for (i in length(S[["a"]])) {
inda<-which(q[h,"a"]==q[h,"d"])
indb<-which(q[i,"a"]==q[i,"e"])
R<-(s[c(inda),"f"]+s[c(indb),"f"])/2
indh<-which(abs(q[,"b"]-R)<=(2*std))
indi<-which(abs(q[,"c"]-R)<=(2*std))
A<-q[indh,c("a","b")]
colnames(A)<- c("a","g")
B<-q[indy,c("a","c")]
colnames(B)<- c("a","g")
C<-rbind(A,B)
c<-density(C[,"a"])
plot(d, col="blue",xlim=c(u,t),ylim=c(0,0.5))
polygon(d, col="blue",border="blue")
lines(c, col="red")
}
}
dev.off()
--
Niklaus Hürlimann
Université de Lausanne
Institut de Minéralogie et Géochimie
L'Anthropole
CH-1015 Lausanne
Suisse
E-mail: niklaus.hurlim...@unil.ch
Tel:+41(0)21 692 4452
[[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.
[R] simple loop iteration
Hi R mailing list,
probably a very basic problem here, I try to do the following:
> Q<-c(1,2,3)
> P<-c(4,5,6)
> A<- data.frame(Q,P)
> A
Q P
1 1 4
2 2 5
3 3 6
this is my simplified data.frame (matrix) now I try to create following
loop for subtraction of element within the data.frame:
> for(i in length(A[,"P"]-1){
delta[i]<- A[i,"P"]-A[i+1,"P"]
}
All I get is a vector of the correct length but with no readings.
Thanks for any help on this.
--
Niklaus Hürlimann
Université de Lausanne
Institut de Minéralogie et Géochimie
L'Anthropole
CH-1015 Lausanne
Suisse
E-mail: niklaus.hurlim...@unil.ch
Tel:+41(0)21 692 4452
[[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.
