[R-es] errores al encontrar parámetros de ajuste iniciales en R
Olivier Nuñez
onunez en iberstat.es
Dom Dic 5 16:58:26 CET 2010
Daniel,
el problema puede ser numérico ya que los valores de la longitud de
onda son próximos a cero y se elevan a la potencia (-3) en la función
de regresión.
Una posible solución consiste en cambiar la escala de las ondas.
Obviamente, la interpretación de los valores de tus parámetros cambia
Lo siguiente funciona (con la temperatura T=1):
> datos=read.csv2("datos.csv",header = TRUE)
> r=10^5 #escala
> fun.ajust<-nls( Densidad ~ (b1/(Longitud*r)^3) / (exp(b2/
(Longitud*r)) - 1 ), data=datos,
+ trace=TRUE, start=c(b1 = .01, b2 =.02) )
36337144 : 0.01 0.02
36329892 : 0.01737898 0.03284804
36323918 : 0.03397161 0.05792367
36292979 : 0.05488755 0.08228909
36259496 : 0.1073412 0.1289229
36242715 : 0.2681239 0.2188145
36112967 : 0.6061090 0.3068132
35903433 : 1.7838682 0.4459792
34712724 : 5.4598786 0.5322451
28658965 : 11.7197031 0.4452781
24193198 : 25.1715444 0.4982637
14199427 : 48.276335 0.491594
5219701 : 83.5498948 0.4951479
2062835 : 118.2081692 0.4929353
2062131 : 118.6022237 0.4941806
2062117 : 118.3849321 0.4938575
2062116 : 118.4424209 0.4939433
2062116 : 118.4272371 0.4939205
2062116 : 118.4312645 0.4939266
2062116 : 118.430198 0.493925
Adjunto el fichero datos.csv que utilicé sobre la base de los datos
que me mandaste.
Un saludo. Olivier
--
____________________________________
Olivier G. Nuñez
Email: onunez en iberstat.es
Tel : +34 663 03 69 09
Web: http://www.iberstat.es
____________________________________
------------ próxima parte ------------
A non-text attachment was scrubbed...
Name: datos.csv
Type: application/octet-stream
Size: 698 bytes
Desc: no disponible
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20101205/ba79f70f/attachment.obj>
------------ próxima parte ------------
El 05/12/2010, a las 16:14, Daniel Arismendi escribió:
>
> Son 42 datos que e obtenido.
>
> Longitud de onda Densidad de energía
> 6.26480e-07 68.44
> 6.93670e-07 97.05
> 7.08640e-07 114.23
> 7.75970e-07 191.52
> 7.98400e-07 211.55
> 8.80990e-07 397.66
> 8.81060e-07 420.57
> 9.93700e-07 683.99
> 9.86330e-07 709.77
> 1.13643e-06 1030.45
> 1.14396e-06 1056.22
> 1.21909e-06 1242.33
> 1.46641e-06 1657.47
> 1.67548e-06 1760.47
> 1.81717e-06 1769.00
> 2.00341e-06 1720.23
> 2.13748e-06 1677.21
> 2.23427e-06 1634.21
> 2.36072e-06 1539.65
> 2.50934e-06 1379.21
> 2.77700e-06 1138.53
> 2.77691e-06 1109.89
> 2.86625e-06 1066.89
> 2.87362e-06 1041.11
> 3.00016e-06 975.19
> 3.10437e-06 917.86
> 3.29785e-06 800.36
> 3.40211e-06 763.08
> 3.49898e-06 742.99
> 3.60323e-06 702.85
> 3.69260e-06 671.31
> 3.96816e-06 573.82
> 4.36287e-06 433.31
> 4.55657e-06 384.54
> 4.75032e-06 352.95
> 4.89933e-06 321.38
> 5.06326e-06 289.80
> 5.22718e-06 255.36
> 5.39859e-06 232.37
> 5.55510e-06 215.12
> 5.69665e-06 180.69
> 6.09174e-06 160.46
>
> El 5 de diciembre de 2010 16:10, Olivier Nuñez <onunez en iberstat.es>
> escribió:
> Daniel,
>
> para que te podamos ayudar,
> creo que lo más sencillo sería que mandes una muestra (o al menos
> una submuestra) de tus datos.
>
> Un saludo. Olivier
>
> -- ____________________________________
>
> Olivier G. Nuñez
> Email: onunez en iberstat.es
> Tel : +34 663 03 69 09
> Web: http://www.iberstat.es
>
> ____________________________________
>
>
>
>
> El 05/12/2010, a las 15:26, Daniel Arismendi escribió:
>
> Buenas Olivier.
>
> Gracias por tu rápida respuesta
>
> Todo lo que has dicho me parece excelente pero el detalle esta en
> que mi data experimental esta tomada con respecto a la longitud de
> onda y la densidad de energía.
>
> Lo que tu afirmas, tomando en cuenta la frecuencia, es una forma de
> escribir la ecuación para la densidad de energía en función de
> dicha frecuencia que básicamente proviene de decir que:
>
>
>
> donde: lamda es la longitud de onda
> c es la velocidad de la luz
> v es la frecuencia.
>
> que sucede si en este caso trabajo con la frecuencia podría estar
> introduciendo errores sistemáticos al problema y pues eso me hace
> pensar que el valor de h para la constante de Planck no sera el mas
> óptimo.
>
>
> puedes corroborarlo viendo este enlace sobre el tema:
> http://en.wikipedia.org/wiki/Planck's_law
>
>
> Lo que estoy tratando entonces de hacer es encontrar los parámetros
> iniciales para la función de ajuste pero por las dos maneras que e
> planteado me encuentro con esos 2 errores que coloque en mi primer
> mensaje y pues e recurrido a ustedes por falta de bibliografía
> acerca de estos en Internet.(que de seguro debe haber pero no
> encuentro una solución clara al problema que se me presenta)
>
> El 5 de diciembre de 2010 15:08, Olivier Nuñez <onunez en iberstat.es>
> escribió:
> Daniel
>
> al contestar rápidamente cometí un error.
> Es más bien la función
>
> (b1, b2) -> { (b1/x_i^3) / (exp(b2/x_i/T) - 1 ) }_{ i = 1..n }
>
> que ha de ser inyectiva para que todo vaya bien (donde los x_i son
> los periodos de onda considerados en tu experimento).
> Asegurarse de que dicha función es inyectiva no es trivial pero con
> un poco de algebra se hace.
> Un saludo. Olivier
>
>
> -- ____________________________________
>
> Olivier G. Nuñez
> Email: onunez en iberstat.es
> Tel : +34 663 03 69 09
> Web: http://www.iberstat.es
>
> ____________________________________
>
>
>
>
> El 05/12/2010, a las 13:39, Daniel Arismendi escribió:
>
> Buenas dias comunidad antes que nada gracias a las ultimas respuestas
> obtenidas por parte de ustedes.
>
> Mientras me leia las paginas que me dejaron para guiarme me dedique
> a buscar
> el libro correspondiente a nlrwr
>
> Nonlinear Regression with R de Christian Ritz ? Jens Carl Streibig
> el cual
> les puedo enviar por correo al que lo desee para que lo lean y lo
> tengan
> dentro de su repertorio bibliografico pues es al que hacen
> referencia en
> estos casos de regresiones no lineales la gente de R
>
> Estuve resolviendo lo del problema del parametro de ajuste para la
> constante
> de Planck mediante minimos cuadrados y pues me tope con estos
> problemas.
>
> En primer lugar como no soy tan agil encontrando los parametros
> para mi
> funcion de ajuste en este caso
>
> fun.ajust<-nls( y ~ (b1/x^5) * (1 / (exp(b2/x*T) - 1) ), data=datos,
> trace=TRUE,
> start=c(b1 = 0.01, b2 = 0.02) )
>
> me encuentro siempre con estos errores :
>
> 1*- Error in nls(... :
> step factor 0.000488281 reduced below 'minFactor' of 0.000976562
> 2*- Error in nls(... :
> singular gradient
> 3*- Error in nlsModel(formula, mf, start, wts) :
> singular gradient matrix at initial parameter estimates
>
> Estos solamente sin usar algun algoritmo para encontrar los minimos
> cuadrados pues cuando use los algoritmos siguientes encontraba otros
> errores:
>
> algorithm: specification of estimation algorithm:
> ? "default": the Gauss-Newton algorithm (the default)
> ? "plinear": for models with conditionally linear parameters
> - "port": for models with constraints on the parameters
> (can be used with the arguments lower/upper)
>
> Dado que tirar piedras no me funciona pues los errores siguen
> saliendo y la
> experiencia es algo que te dice como podrian ser los parametros
> iniciales de
> ajuste me dedique a usar 2 maneras a traves de las cuales pudiera
> encontrar
> mis parametros de ajuste.
>
>
> 1*-Hago uso de la libreria nlstools
> library(nlstools)
> modelo<-V2~(b1/V1^5) * (1 / (exp(b2/V1*T) - 1)
> preview(formula = modelo, data = datos, start = start =
> list(b1=0.01,b2=0.02))
> Error: no se pudo encontrar la función "preview"
>
> Referencia a la funcion preview:
> Details
>
> The function preview helps defining the parameter starting values
> prior
> fitting the model. It provides a superimposed plot of observed
> (circles) and
> predicted (crosses) values of the dependent variable versus one of the
> independent variables with the model evaluated at the starting
> values of the
> parameters. The function overview returns the parameters estimates,
> their
> standard errors as well as their asymptotic confidence intervals
> and the
> correlation matrix (alternately, the function confint provides better
> confidence interval estimates whenever it converges).
> plotfitdisplays a
> superimposed plot of the dependent variable versus one the independent
> variables together with the fitted model
>
>
> 2*- parecida a preview pero haciendo uso de una funcion predefinida
> y la
> funcion curve superpuesta sobre los puntos de los datos
> experimentales:
>
> datos<-read.table("datos.txt")
>
> attach(datos)
>
> plot(V2~V1,data=datos,xlab="longitud de onda",ylab="densidad de
> energia")
> T<-1595
> modelo1<- function(V1,b1,b2){(b1/V1^5) * (1 / (exp(b2/V1*T) - 1) )}
> plot(V2~V1,data=datos,ylim=c(0,1.85e3),xlab="longitud de
> onda",ylab="densidad de energia")
> curve(modelo(V1,b1=0.01,b2=0.02),add=TRUE)
> Error en curve(modelo(V1, b1 = 0.01, b2 = 0.02), add = TRUE) :
> 'expr' must be a function, call or an expression containing 'x'
>
>
> Referencia a esta segunda forma la encontre en el libro de nlrwr.
>
> Agradeceria si alguien conoce de estos errores como podria
> solucionarlos
> pues es poca la bibliografia respecto a estos en la web y no logro
> terminar
> de hacer la funcion de ajuste.
>
> Saludos a toda la comunidad.
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
>
>
>
Más información sobre la lista de distribución R-help-es