R help es - Dec 2010 - errores al encontrar parámetros de ajuste iniciales en R

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]]

Daniel,

según la expresión de la  ley de Planck, encontrada en la Wikipedia,

http://es.wikipedia.org/wiki/Ley_de_Planck


tu variable x debe estar relacionada con el periodo de la onda y T es  
la temperatura,
Si es así, me parece que el problema es simplemente que tu modelo no  
está bien escrito y debería ser más bien algo así como:


y ~ (b1/x^3) / (exp(b2/x/T) - 1 )

está función es continua y estrictamente decreciente (además es  
convexa),
por lo tanto es inyectiva, lo cual asegura  que al menos  
asintóticamente el estimador de mínimos cuadrados de tus parámetros  
converge y no tiene sesgo.
En fin no debería haber errores, a no ser que tengas muy pocos puntos  
y que estén demasiado cercanos (periodos de la onda similares).

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ó:
> (b1/x^5) * (1 / (exp(b2/x*T)

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@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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es

	[[alternative HTML version deleted]]

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:

[image: \lambda=\frac{c}{\nu}]

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)
<http://en.wikipedia.org/wiki/Planck''s_law>
El 5 de diciembre de 2010 15:08, Olivier Nuñez
<onunez@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@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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
>
>
	[[alternative HTML version deleted]]

Hola Daniel,

Acabo de ver que este hilo ha continuado y quería dar un par de detalles de
los resultados que obtengo a partir de los datos experimentales que has
aportado y de los ajustes del modelo que ha dado Olivier.

1. El modelo que he ajustado y los resultados son estos:
> fun.ajust<-nls( Densidad ~ (b1/(Longitud*r)^3) / (exp(b2/(Longitud*r*T))
-
1 ), data=datos,
+                 trace=TRUE, start=list(b1 = 0.01, b2 =0.02)  )
175174548787 :  0.01 0.02
166560945244 :  0.4610861 0.9420121
130902950887 :   6.364498 13.889852
32872478327 :   6.840776 27.923805
7805195855 :   7.408324 54.668227
1743852966 :    8.63834 105.06266
339150528 :   11.58223 196.54116
47775729 :   19.66186 355.41115
6544606 :   44.83457 608.58879
4927633 :  104.1179 844.4815
2561367 :  106.9241 738.4609
2066585 :  120.5978 792.7247
2064439 :  117.7953 786.3291
2064278 :  118.5403 788.0546
2064267 :  118.3512 787.5998
2064266 :  118.4019 787.7205
2064266 :  118.3885 787.6885
2064266 :  118.3920 787.6970
2064266 :  118.3911 787.6948
>
> coef(fun.ajust)
      b1       b2
118.3911 787.6948
> summary(fun.ajust)
Formula: Densidad ~ (b1/(Longitud * r)^3)/(exp(b2/(Longitud * r * T)) -
    1)

Parameters:
   Estimate Std. Error t value Pr(>|t|)
b1   118.39      12.08   9.804 3.41e-12 ***
b2   787.69      26.69  29.509  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 227.2 on 40 degrees of freedom

Number of iterations to convergence: 18
Achieved convergence tolerance: 3.514e-06

########### Discusión de resultados ##################
*A)* En mi caso los valores que ajustan el modelo son:
b1= 118.39
b2=787.69

*B)* Como ves he iniciado el modelo a partir de valores de b1=0.01 y
b2=0.02.
Cuando iniciamos este hilo de conversación utilicé estos mismos datos
iniciales, porque aparecían en un ejemplo de la función nls.
Pero ves que nls va ajustando el modelo a partir de esos datos hasta llagar
a los 118 y 787 para b1 y b2 respectivamente.

*C) *Para calcular algún límite de confianza sobre estos valores de b1 y de
b2, puedes utilizar bien la función nlsBoot o nlsJack de la librería
nlstools.

Este es el código:

boot.dat <- nlsBoot(fun.ajust)
plot(boot.dat)
plot(boot.dat, type = "boxplot", ask=FALSE)
summary.nlsBoot(boot.dat)

jack.dat <- nlsJack(fun.ajust)
plot(jack.dat)
summary(jack.dat)


 Los resultados que obtienes son:
> summary.nlsBoot(boot.dat)
------
Bootstrap estimates
      b1       b2
119.3822 789.6808

------
Bootstrap confidence intervals
       2.5%    97.5%
b1  99.3526 141.6616
b2 740.1372 835.6497

> summary(jack.dat)
------
Jackknife estimates
      b1       b2
118.9262 788.8243

------
Jackknife confidence intervals
        Low       Up
b1 100.7941 137.0584
b2 751.9638 825.6847



*D)* Con estos datos del modelo, puedes pintar tanto los resultados
experimentales como los de la función "ajustada" con estos parámetros
b1 y
b2.

Este es el código:

##### Pinto teorico y real
plot(Densidad~Longitud, data=datos, xlab="Longitud de Onda",
ylab="Densidad
Energía", col="blue")
x.teo<-seq(0,10e-6,0.05e-6)
b1.exp<-as.vector(coef(fun.ajust)[1])
b2.exp<-as.vector(coef(fun.ajust)[2])
#b1.exp<-42.19
#b2.exp<-1388
y.teo<-(b1.exp/(x.teo*r)^3) / (exp(b2.exp/(x.teo*r*T))-1)
points( x.teo, y.teo, col="red")

En el código tengo comentados los valores que aparecen en tu último correo.
Verás que con los datos calculados, el ajuste *no* es muy bueno para toda la
zona del máximo de la curva y toda la porción decreciente.

Con los datos de tu correo, creo que debe de haber un error, el ajuste con
ellos se aleja mucho de la curva experimental.

*E)* Y con estos valores, el cálculo de h sería el siguiente:

#########################
### Calculo de h
k.dat<-1.3806504e-23
c.dat<-3e8

#b2 = hc/k
#h= k b2/c

#b1=8*pi*h*c
# h = b1 / (8*pi*c)

h.exp<- k.dat * b2.exp / c.dat
h.exp.2<- b1.exp / (8 * 3.14 * c.dat)
print(h.exp)
print(h.exp.2)

Las variables k.dat y c.dat contienen los valores de la constante de
Boltzman y la velocidad de la luz, en el Sistema Internacional.
Entiendo que en tus cálculos tú también has utilizado el mismo sistema de
referencia.

Los resultados que obtengo mediante b1 (h.exp.2) y b2 (h.exp) son los
siguientes:
> print(h.exp)
[1] 3.625104e-29
> print(h.exp.2)
[1] 1.571007e-08
>
Un detalle que aunque no supone un cambio sustancial sobre los resultados,
pero sí para tenerlo en cuenta es que realmente el significado de b1. b1 = 8
* pi * h.

El primer factor de la energía de cuerpo negro es: 8 pi h f^3 / c^3, si la
frecuencia (f) se sustituye por c/l (l es lambda), se tiene: 8 pi h/l^3.

La sustitución que hemos hecho es de b1= 8*pi*h


*F) * He intentado reproducir el ajuste sin hacer la transformación que ha
sugerido Olivier, pero he obtenido el mismo tipo de errores que tú.
>
>
names(datos)<-c(''Longitud'',''Densidad'')
> T<-1595
>
> nls.control(maxiter = 150, tol = 1e-15, minFactor = 1/2024,
+             printEval = FALSE, warnOnly = FALSE)
$maxiter
[1] 150

$tol
[1] 1e-15

$minFactor
[1] 0.0004940711

$printEval
[1] FALSE

$warnOnly
[1] FALSE
>
> fun.ajust<-nls( Densidad ~ (b1/(Longitud)^3) / (exp(b2/(Longitud*T)) - 1
), data=datos,
+                 trace=TRUE, start=list(b1 = 1, b2 = 7) )
Error en nlsModel(formula, mf, start, wts) :
  singular gradient matrix at initial parameter
estimates
>
En la función "nls.control()" he reducido el paso, e incrementado el
número
posible de iteraciones para poder aumentar la precisión del cálculo... pero
como ves obtengo error.
Igualmente, ves que he comenzado desde otros puntos b1=1 y b2=7.
He probado con b1=100 y b2=700, que son valores próximo al otro cálculo y el
error es el mismo.



Saludos,
Carlos Ortega
www.qualityexcellence.es
http://www.datanalytics.com/blog/

2010/12/5 Daniel Arismendi <darismendi@cmc.org.ve>
> Efectivamente con la forma en que tu presentas esto cuadra perfectamente
> con
> el detalle de que  es (b1/(Longitud*r)^5) y lo otro T no puede ser 1 pues
> nunca encontraría el espectro de radiación de cuerpo negro, pero no importa
> yo le cambie (exp(b2/(Longitud*r*T) donde T es 1595.
>
> Pues cuando intento obtener h se que b1=8*pi*h*c
> en este caso el valor de b1 que obtengo es de 42.19
> y para b2 es 1388.23 donde b2 = hc/k
>
> Siendo k la constante de Boltzmann c la velocidad de la luz y h la
> constante
> de Planck que necesito hallar.
>
> Si la obtengo de b1  tengo un valor de h= 6e-9
> Si la obtengo de b2 h= 6.388868e-29
>
> Si valor a obtener debería ser próximo a  h= 6.2606896e-34 cosa que en mis
> resultados anteriores es parecido pero los ordenes de magnitud son muy
> diferentes.
>
> Ahora mi pregunta es: hay alguna forma de que R trabaje los valores de mis
> longitud de onda con una precisión mucho mayor? para que pueda trabajar con
> la regresión lineal de forma normal.
>
> Pues al trabajar con los valores multiplicados por 10^5 no es lo mejor y
> pues como decías tu la interpretación de los valores de mis parámetros
> cambia pero se debe a que estamos introduciendo errores sistemáticas
> asociados a un calculo indirecto para encontrar la regresión no lineal.
>
> Otra cosa Olivier como hiciste para darle los parámetros iniciales a la
> regresión no lineal? pues veo que pusiste b1=.01 y b2=.02 los encontraste
> de
> forma empírica o llegaste a un valor óptimo asociándolos a las dos maneras
> de encontrar los parámetros iniciales como lo propuse en mi mensaje
> inicial.?
>
>
> De igual forma gracias nuevamente por tu respuesta
> Saludos Olivier.
>
> El 5 de diciembre de 2010 16:58, Olivier Nuñez <onunez@iberstat.es
> >escribió:
>
> > 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@iberstat.es
> > Tel : +34 663 03 69 09
> > Web: http://www.iberstat.es
> >
> > ____________________________________
> >
> >
> >
> >
> >
> > 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@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@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@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@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@r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-help-es
> >>
> >>
> >>
> >>
> >>
> >
> >
>
>        [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>
>
	[[alternative HTML version deleted]]

Carlos me pareció excelente la forma en que lo desglosas todo se nota la
diferencia entre tu un usuario experimentado y yo una usuario que se inicia
en el mundo de R gracias por esos detalles.

Ahora a lo que voy y ya para ir finalizando este hilo.

De acuerdo a lo que empecé mi primer mensaje:
La idea era encontrar los parámetros iniciales mediante una forma que no
fuera la de colocar los parámetros iniciales que me ofrecía la data del
paquete NISTLS pues a pesar de que la curva que estos presentaban tiene
características muy similares a la que yo obtengo experimentalmente, mi idea
era haciendo ajustes directamente sobre la gráfica con estas dos maneras
pero me encuentro con estos errores que en un principio fueron los que
plantee al escribir el primer mensaje.


1*-Hago uso de la librería 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"
2*- parecida a preview pero haciendo uso de una función predefinida y la
función 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''

De estas dos formas que coloco alguno la ha utilizado o conoce el porque se
podrían estar presentando estos errores.
pues e cargado las librerías y todo lo demás pero nada que logro encontrar
los parámetros iniciales óptimos para el ajuste.

El propósito es  obtener los parámetros iniciales de tal forma que pudiera
presentar una gráfica en la que se observara como cambiando los parámetros
iniciales me acerco de mejor forma a la curva de ajuste óptima.

Gracias nuevamente y espero poder ayudarle en lo que necesiten ahora y en un
futuro cuando tenga mucha mas experiencia.

El 5 de diciembre de 2010 23:37, Carlos Ortega <coforfe@gmail.com>
escribió:
> Hola Daniel,
>
> Acabo de ver que este hilo ha continuado y quería dar un par de detalles de
> los resultados que obtengo a partir de los datos experimentales que has
> aportado y de los ajustes del modelo que ha dado Olivier.
>
> 1. El modelo que he ajustado y los resultados son estos:
>
> > fun.ajust<-nls( Densidad ~ (b1/(Longitud*r)^3) /
(exp(b2/(Longitud*r*T))
> - 1 ), data=datos,
> +                 trace=TRUE, start=list(b1 = 0.01, b2 =0.02)  )
> 175174548787 :  0.01 0.02
> 166560945244 :  0.4610861 0.9420121
> 130902950887 :   6.364498 13.889852
> 32872478327 :   6.840776 27.923805
> 7805195855 :   7.408324 54.668227
> 1743852966 :    8.63834 105.06266
> 339150528 :   11.58223 196.54116
> 47775729 :   19.66186 355.41115
> 6544606 :   44.83457 608.58879
> 4927633 :  104.1179 844.4815
> 2561367 :  106.9241 738.4609
> 2066585 :  120.5978 792.7247
> 2064439 :  117.7953 786.3291
> 2064278 :  118.5403 788.0546
> 2064267 :  118.3512 787.5998
> 2064266 :  118.4019 787.7205
> 2064266 :  118.3885 787.6885
> 2064266 :  118.3920 787.6970
> 2064266 :  118.3911 787.6948
> >
> > coef(fun.ajust)
>       b1       b2
> 118.3911 787.6948
> > summary(fun.ajust)
>
> Formula: Densidad ~ (b1/(Longitud * r)^3)/(exp(b2/(Longitud * r * T)) -
>     1)
>
> Parameters:
>    Estimate Std. Error t value Pr(>|t|)
> b1   118.39      12.08   9.804 3.41e-12 ***
> b2   787.69      26.69  29.509  < 2e-16 ***
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Residual standard error: 227.2 on 40 degrees of freedom
>
> Number of iterations to convergence: 18
> Achieved convergence tolerance: 3.514e-06
>
> ########### Discusión de resultados ##################
> *A)* En mi caso los valores que ajustan el modelo son:
> b1= 118.39
> b2=787.69
>
> *B)* Como ves he iniciado el modelo a partir de valores de b1=0.01 y
> b2=0.02.
> Cuando iniciamos este hilo de conversación utilicé estos mismos datos
> iniciales, porque aparecían en un ejemplo de la función nls.
> Pero ves que nls va ajustando el modelo a partir de esos datos hasta llagar
> a los 118 y 787 para b1 y b2 respectivamente.
>
> *C) *Para calcular algún límite de confianza sobre estos valores de b1 y
> de b2, puedes utilizar bien la función nlsBoot o nlsJack de la librería
> nlstools.
>
> Este es el código:
>
> boot.dat <- nlsBoot(fun.ajust)
> plot(boot.dat)
> plot(boot.dat, type = "boxplot", ask=FALSE)
> summary.nlsBoot(boot.dat)
>
> jack.dat <- nlsJack(fun.ajust)
> plot(jack.dat)
> summary(jack.dat)
>
>
>  Los resultados que obtienes son:
>
> > summary.nlsBoot(boot.dat)
>
> ------
> Bootstrap estimates
>       b1       b2
> 119.3822 789.6808
>
> ------
> Bootstrap confidence intervals
>        2.5%    97.5%
> b1  99.3526 141.6616
> b2 740.1372 835.6497
>
>
> > summary(jack.dat)
>
> ------
> Jackknife estimates
>       b1       b2
> 118.9262 788.8243
>
> ------
> Jackknife confidence intervals
>         Low       Up
> b1 100.7941 137.0584
> b2 751.9638 825.6847
>
>
>
> *D)* Con estos datos del modelo, puedes pintar tanto los resultados
> experimentales como los de la función "ajustada" con estos
parámetros b1 y
> b2.
>
> Este es el código:
>
> ##### Pinto teorico y real
> plot(Densidad~Longitud, data=datos, xlab="Longitud de Onda",
ylab="Densidad
> Energía", col="blue")
> x.teo<-seq(0,10e-6,0.05e-6)
> b1.exp<-as.vector(coef(fun.ajust)[1])
> b2.exp<-as.vector(coef(fun.ajust)[2])
> #b1.exp<-42.19
> #b2.exp<-1388
> y.teo<-(b1.exp/(x.teo*r)^3) / (exp(b2.exp/(x.teo*r*T))-1)
> points( x.teo, y.teo, col="red")
>
> En el código tengo comentados los valores que aparecen en tu último correo.
> Verás que con los datos calculados, el ajuste *no* es muy bueno para toda
> la zona del máximo de la curva y toda la porción decreciente.
>
> Con los datos de tu correo, creo que debe de haber un error, el ajuste con
> ellos se aleja mucho de la curva experimental.
>
> *E)* Y con estos valores, el cálculo de h sería el siguiente:
>
> #########################
> ### Calculo de h
> k.dat<-1.3806504e-23
> c.dat<-3e8
>
> #b2 = hc/k
> #h= k b2/c
>
> #b1=8*pi*h*c
> # h = b1 / (8*pi*c)
>
> h.exp<- k.dat * b2.exp / c.dat
> h.exp.2<- b1.exp / (8 * 3.14 * c.dat)
> print(h.exp)
> print(h.exp.2)
>
> Las variables k.dat y c.dat contienen los valores de la constante de
> Boltzman y la velocidad de la luz, en el Sistema Internacional.
> Entiendo que en tus cálculos tú también has utilizado el mismo sistema de
> referencia.
>
> Los resultados que obtengo mediante b1 (h.exp.2) y b2 (h.exp) son los
> siguientes:
>
> > print(h.exp)
> [1] 3.625104e-29
> > print(h.exp.2)
> [1] 1.571007e-08
> >
>
> Un detalle que aunque no supone un cambio sustancial sobre los resultados,
> pero sí para tenerlo en cuenta es que realmente el significado de b1. b1 =
8
> * pi * h.
>
> El primer factor de la energía de cuerpo negro es: 8 pi h f^3 / c^3, si la
> frecuencia (f) se sustituye por c/l (l es lambda), se tiene: 8 pi h/l^3.
>
> La sustitución que hemos hecho es de b1= 8*pi*h
>
>
> *F) * He intentado reproducir el ajuste sin hacer la transformación que ha
> sugerido Olivier, pero he obtenido el mismo tipo de errores que tú.
>
> >
> >
names(datos)<-c(''Longitud'',''Densidad'')
> > T<-1595
> >
> > nls.control(maxiter = 150, tol = 1e-15, minFactor = 1/2024,
> +             printEval = FALSE, warnOnly = FALSE)
> $maxiter
> [1] 150
>
> $tol
> [1] 1e-15
>
> $minFactor
> [1] 0.0004940711
>
>  $printEval
> [1] FALSE
>
> $warnOnly
> [1] FALSE
>
> >
> > fun.ajust<-nls( Densidad ~ (b1/(Longitud)^3) /
(exp(b2/(Longitud*T)) - 1
> ), data=datos,
> +                 trace=TRUE, start=list(b1 = 1, b2 = 7) )
> Error en nlsModel(formula, mf, start, wts) :
>   singular gradient matrix at initial parameter estimates
> >
>
> En la función "nls.control()" he reducido el paso, e incrementado
el número
> posible de iteraciones para poder aumentar la precisión del cálculo... pero
> como ves obtengo error.
> Igualmente, ves que he comenzado desde otros puntos b1=1 y b2=7.
> He probado con b1=100 y b2=700, que son valores próximo al otro cálculo y
> el error es el mismo.
>
>
>
> Saludos,
> Carlos Ortega
> www.qualityexcellence.es
> http://www.datanalytics.com/blog/
>
> 2010/12/5 Daniel Arismendi <darismendi@cmc.org.ve>
>
> Efectivamente con la forma en que tu presentas esto cuadra perfectamente
>> con
>> el detalle de que  es (b1/(Longitud*r)^5) y lo otro T no puede ser 1
pues
>> nunca encontraría el espectro de radiación de cuerpo negro, pero no
>> importa
>> yo le cambie (exp(b2/(Longitud*r*T) donde T es 1595.
>>
>> Pues cuando intento obtener h se que b1=8*pi*h*c
>> en este caso el valor de b1 que obtengo es de 42.19
>> y para b2 es 1388.23 donde b2 = hc/k
>>
>> Siendo k la constante de Boltzmann c la velocidad de la luz y h la
>> constante
>> de Planck que necesito hallar.
>>
>> Si la obtengo de b1  tengo un valor de h= 6e-9
>> Si la obtengo de b2 h= 6.388868e-29
>>
>> Si valor a obtener debería ser próximo a  h= 6.2606896e-34 cosa que en
mis
>> resultados anteriores es parecido pero los ordenes de magnitud son muy
>> diferentes.
>>
>> Ahora mi pregunta es: hay alguna forma de que R trabaje los valores de
mis
>> longitud de onda con una precisión mucho mayor? para que pueda trabajar
>> con
>> la regresión lineal de forma normal.
>>
>> Pues al trabajar con los valores multiplicados por 10^5 no es lo mejor
y
>> pues como decías tu la interpretación de los valores de mis parámetros
>> cambia pero se debe a que estamos introduciendo errores sistemáticas
>> asociados a un calculo indirecto para encontrar la regresión no lineal.
>>
>> Otra cosa Olivier como hiciste para darle los parámetros iniciales a la
>> regresión no lineal? pues veo que pusiste b1=.01 y b2=.02 los
encontraste
>> de
>> forma empírica o llegaste a un valor óptimo asociándolos a las dos
maneras
>> de encontrar los parámetros iniciales como lo propuse en mi mensaje
>> inicial.?
>>
>>
>> De igual forma gracias nuevamente por tu respuesta
>> Saludos Olivier.
>>
>> El 5 de diciembre de 2010 16:58, Olivier Nuñez <onunez@iberstat.es
>> >escribió:
>>
>> > 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@iberstat.es
>> > Tel : +34 663 03 69 09
>> > Web: http://www.iberstat.es
>> >
>> > ____________________________________
>> >
>> >
>> >
>> >
>> >
>> > 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@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@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<http://en.wikipedia.org/wiki/Planck%27s_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@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@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@r-project.org
>> >> https://stat.ethz.ch/mailman/listinfo/r-help-es
>> >>
>> >>
>> >>
>> >>
>> >>
>> >
>> >
>>
>>        [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-help-es mailing list
>> R-help-es@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>
>>
>
	[[alternative HTML version deleted]]

Hola Daniel,

Ayer lo había visto en la Wikipedia, pero era ya muy tarde para mirarlo más
despacio. Pero hoy en tu última entrada aparece de nuevo. La dependencia de
la Energía es con lambda^5 y no a la 3 como había hecho.

He repetido los cálculos y obtengo vuestros mismos resultados:

> fun.ajust<-nls( Densidad ~ (b1/(Longitud*r)^5) / (exp(b2/(Longitud*r*T))
-
1 ), data=datos,
+                 trace=TRUE, start=list(b1 = 0.01, b2 =0.02)  )
6.28589e+15 :  0.01 0.02
1.007657e+13 :  0.000417183 0.020834104
7.45834e+12 :  0.002579168 0.149620578
6.064184e+12 :  0.02519974 1.61045063
1.490792e+12 :  0.02496692 3.19389164
370967482202 :  0.0252321 6.3732852
91831608255 :   0.02575734 12.68644838
22501368337 :   0.02683438 25.13917990
5401393817 :   0.02912530 49.41740581
1243373917 :   0.03430239 95.90404916
264551276 :    0.0475513 183.3972677
55448124 :    0.09146908 352.13942960
30830335 :    0.3499503 755.8767479
30815787 :     1.123102 1043.435579
30230174 :     3.915183 1366.382307
18962968 :    13.45007 1420.04000
3137772 :    27.65976 1363.23174
77909.39 :    42.01506 1397.57997
29043.31 :    42.15555 1387.67983
29016.14 :    42.18324 1388.09104
>
> coef(fun.ajust)
        b1         b2
  42.18324 1388.09104
> summary(fun.ajust)
Formula: Densidad ~ (b1/(Longitud * r)^5)/(exp(b2/(Longitud * r * T)) -
    1)

Parameters:
    Estimate Std. Error t value Pr(>|t|)
b1   42.1832     0.5446   77.45   <2e-16 ***
b2 1388.0910     3.6865  376.54   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.93 on 40 degrees of freedom

Number of iterations to convergence: 19
Achieved convergence tolerance: 5.801e-06


*B)* Intervalos de confianza de b1 y b2
Además de las funciones nlsBoot() y nlsJack(), hay otra función de la
librería nlrwr que permite el cálculo directo de los intervalos de confianza
para los parámetros del ajuste. La función se llama confint2().

Como verás los resultados de los tres métodos de cálculo son equivalentes.
> confint2(fun.ajust)
        2.5 %     97.5 %
b1   41.08252   43.28396
b2 1380.64043 1395.54165

> summary.nlsBoot(boot.dat)
------
Bootstrap estimates
        b1         b2
  42.18289 1388.22464

------
Bootstrap confidence intervals
         2.5%      97.5%
b1   41.16342   43.23346
b2 1381.38576 1395.24134
> summary(jack.dat)
------
Jackknife estimates
       b1        b2
  42.1864 1388.1491

------
Jackknife confidence intervals
          Low         Up
b1   40.86558   43.50722
b2 1380.77507 1395.52312


*C)* Gráficos numérico y experimentales - Uso de Curve.
Con estos nuevos resultados, el ajuste numérico y experimental son muy, muy
buenos.
Al final verás el ejemplo de cómo he utilizado curve.

##### Pinto teorico y real
plot(Densidad~Longitud, data=datos, xlab="Longitud de Onda",
ylab="Densidad
Energía",
     col="blue", pch="+", lty=2)
x.teo<-seq(0,10e-6,0.05e-6)
b1.exp<-as.vector(coef(fun.ajust)[1])
b2.exp<-as.vector(coef(fun.ajust)[2])
#b1.exp<-110
#b2.exp<-720
y.teo<-(b1.exp/(x.teo*r)^5) / (exp(b2.exp/(x.teo*r*T))-1)
#points( x.teo, y.teo, col="red")

curve( b1.exp/(x*r)^5 / (exp( b2.exp / (x*r*T) )-1), 0, 1e-5, n=500, add=T,
lwd=2, col="violet")

*D)* Resultados de la constante de Planck:
Ahora sí, el valor de h a partir del resultado de b2 se aproxima bastante al
valor de la constante.
>
> #########################
> ### Calculo de h
> k.dat<-1.3806504e-23
> c.dat<-3e8
>
> #b2 = hc/k
> #h= k b2/c
>
> #b1=8*pi*h*c
> # h = b1 / (8*pi*c)
>
> h.exp<- k.dat * b2.exp / c.dat
> #h.exp.2<- b1.exp / (8 * 3.14 * c.dat)
> h.exp.2<- b1.exp / (8 * 3.14 )
> print(h.exp)
[1] 6.388228e-29
> print(h.exp.2)
[1] 1.679269
>



*E)* Uso de preview.
Ayer también intenté usar preview para obtener unos parámetros iniciales de
arranque para el modelo, pero también tuve problemas parecidos a los tuyos.
Voy a seguir intentándolo un poco más.


Saludos,
Carlos Ortega
www.qualityexcellence.es
http://www.datanalytics.com/blog/




2010/12/6 Daniel Arismendi <darismendi@cmc.org.ve>
> Carlos me pareció excelente la forma en que lo desglosas todo se nota la
> diferencia entre tu un usuario experimentado y yo una usuario que se inicia
> en el mundo de R gracias por esos detalles.
>
> Ahora a lo que voy y ya para ir finalizando este hilo.
>
> De acuerdo a lo que empecé mi primer mensaje:
> La idea era encontrar los parámetros iniciales mediante una forma que no
> fuera la de colocar los parámetros iniciales que me ofrecía la data del
> paquete NISTLS pues a pesar de que la curva que estos presentaban tiene
> características muy similares a la que yo obtengo experimentalmente, mi
idea
> era haciendo ajustes directamente sobre la gráfica con estas dos maneras
> pero me encuentro con estos errores que en un principio fueron los que
> plantee al escribir el primer mensaje.
>
>
> 1*-Hago uso de la librería 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"
>
> 2*- parecida a preview pero haciendo uso de una función predefinida y la
> función 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''
>
> De estas dos formas que coloco alguno la ha utilizado o conoce el porque se
> podrían estar presentando estos errores.
> pues e cargado las librerías y todo lo demás pero nada que logro encontrar
> los parámetros iniciales óptimos para el ajuste.
>
> El propósito es  obtener los parámetros iniciales de tal forma que pudiera
> presentar una gráfica en la que se observara como cambiando los parámetros
> iniciales me acerco de mejor forma a la curva de ajuste óptima.
>
> Gracias nuevamente y espero poder ayudarle en lo que necesiten ahora y en
> un futuro cuando tenga mucha mas experiencia.
>
> El 5 de diciembre de 2010 23:37, Carlos Ortega
<coforfe@gmail.com>escribió:
>
>  Hola Daniel,
>>
>> Acabo de ver que este hilo ha continuado y quería dar un par de
detalles
>> de los resultados que obtengo a partir de los datos experimentales que
has
>> aportado y de los ajustes del modelo que ha dado Olivier.
>>
>> 1. El modelo que he ajustado y los resultados son estos:
>>
>> > fun.ajust<-nls( Densidad ~ (b1/(Longitud*r)^3) /
(exp(b2/(Longitud*r*T))
>> - 1 ), data=datos,
>> +                 trace=TRUE, start=list(b1 = 0.01, b2 =0.02)  )
>> 175174548787 :  0.01 0.02
>> 166560945244 :  0.4610861 0.9420121
>> 130902950887 :   6.364498 13.889852
>> 32872478327 :   6.840776 27.923805
>> 7805195855 :   7.408324 54.668227
>> 1743852966 :    8.63834 105.06266
>> 339150528 :   11.58223 196.54116
>> 47775729 :   19.66186 355.41115
>> 6544606 :   44.83457 608.58879
>> 4927633 :  104.1179 844.4815
>> 2561367 :  106.9241 738.4609
>> 2066585 :  120.5978 792.7247
>> 2064439 :  117.7953 786.3291
>> 2064278 :  118.5403 788.0546
>> 2064267 :  118.3512 787.5998
>> 2064266 :  118.4019 787.7205
>> 2064266 :  118.3885 787.6885
>> 2064266 :  118.3920 787.6970
>> 2064266 :  118.3911 787.6948
>> >
>> > coef(fun.ajust)
>>       b1       b2
>> 118.3911 787.6948
>> > summary(fun.ajust)
>>
>> Formula: Densidad ~ (b1/(Longitud * r)^3)/(exp(b2/(Longitud * r * T)) -
>>     1)
>>
>> Parameters:
>>    Estimate Std. Error t value Pr(>|t|)
>> b1   118.39      12.08   9.804 3.41e-12 ***
>> b2   787.69      26.69  29.509  < 2e-16 ***
>> ---
>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>> Residual standard error: 227.2 on 40 degrees of freedom
>>
>> Number of iterations to convergence: 18
>> Achieved convergence tolerance: 3.514e-06
>>
>> ########### Discusión de resultados ##################
>> *A)* En mi caso los valores que ajustan el modelo son:
>> b1= 118.39
>> b2=787.69
>>
>> *B)* Como ves he iniciado el modelo a partir de valores de b1=0.01 y
>> b2=0.02.
>> Cuando iniciamos este hilo de conversación utilicé estos mismos datos
>> iniciales, porque aparecían en un ejemplo de la función nls.
>> Pero ves que nls va ajustando el modelo a partir de esos datos hasta
>> llagar a los 118 y 787 para b1 y b2 respectivamente.
>>
>> *C) *Para calcular algún límite de confianza sobre estos valores de b1
y
>> de b2, puedes utilizar bien la función nlsBoot o nlsJack de la librería
>> nlstools.
>>
>> Este es el código:
>>
>> boot.dat <- nlsBoot(fun.ajust)
>> plot(boot.dat)
>> plot(boot.dat, type = "boxplot", ask=FALSE)
>> summary.nlsBoot(boot.dat)
>>
>> jack.dat <- nlsJack(fun.ajust)
>> plot(jack.dat)
>> summary(jack.dat)
>>
>>
>>  Los resultados que obtienes son:
>>
>> > summary.nlsBoot(boot.dat)
>>
>> ------
>> Bootstrap estimates
>>       b1       b2
>> 119.3822 789.6808
>>
>> ------
>> Bootstrap confidence intervals
>>        2.5%    97.5%
>> b1  99.3526 141.6616
>> b2 740.1372 835.6497
>>
>>
>> > summary(jack.dat)
>>
>> ------
>> Jackknife estimates
>>       b1       b2
>> 118.9262 788.8243
>>
>> ------
>> Jackknife confidence intervals
>>         Low       Up
>> b1 100.7941 137.0584
>> b2 751.9638 825.6847
>>
>>
>>
>> *D)* Con estos datos del modelo, puedes pintar tanto los resultados
>> experimentales como los de la función "ajustada" con estos
parámetros b1 y
>> b2.
>>
>> Este es el código:
>>
>> ##### Pinto teorico y real
>> plot(Densidad~Longitud, data=datos, xlab="Longitud de Onda",
>> ylab="Densidad Energía", col="blue")
>> x.teo<-seq(0,10e-6,0.05e-6)
>> b1.exp<-as.vector(coef(fun.ajust)[1])
>> b2.exp<-as.vector(coef(fun.ajust)[2])
>> #b1.exp<-42.19
>> #b2.exp<-1388
>> y.teo<-(b1.exp/(x.teo*r)^3) / (exp(b2.exp/(x.teo*r*T))-1)
>> points( x.teo, y.teo, col="red")
>>
>> En el código tengo comentados los valores que aparecen en tu último
>> correo.
>> Verás que con los datos calculados, el ajuste *no* es muy bueno para
toda
>> la zona del máximo de la curva y toda la porción decreciente.
>>
>> Con los datos de tu correo, creo que debe de haber un error, el ajuste
con
>> ellos se aleja mucho de la curva experimental.
>>
>> *E)* Y con estos valores, el cálculo de h sería el siguiente:
>>
>> #########################
>> ### Calculo de h
>> k.dat<-1.3806504e-23
>> c.dat<-3e8
>>
>> #b2 = hc/k
>> #h= k b2/c
>>
>> #b1=8*pi*h*c
>> # h = b1 / (8*pi*c)
>>
>> h.exp<- k.dat * b2.exp / c.dat
>> h.exp.2<- b1.exp / (8 * 3.14 * c.dat)
>> print(h.exp)
>> print(h.exp.2)
>>
>> Las variables k.dat y c.dat contienen los valores de la constante de
>> Boltzman y la velocidad de la luz, en el Sistema Internacional.
>> Entiendo que en tus cálculos tú también has utilizado el mismo sistema
de
>> referencia.
>>
>> Los resultados que obtengo mediante b1 (h.exp.2) y b2 (h.exp) son los
>> siguientes:
>>
>> > print(h.exp)
>> [1] 3.625104e-29
>> > print(h.exp.2)
>> [1] 1.571007e-08
>> >
>>
>> Un detalle que aunque no supone un cambio sustancial sobre los
resultados,
>> pero sí para tenerlo en cuenta es que realmente el significado de b1.
b1 = 8
>> * pi * h.
>>
>> El primer factor de la energía de cuerpo negro es: 8 pi h f^3 / c^3, si
la
>> frecuencia (f) se sustituye por c/l (l es lambda), se tiene: 8 pi
h/l^3.
>>
>> La sustitución que hemos hecho es de b1= 8*pi*h
>>
>>
>> *F) * He intentado reproducir el ajuste sin hacer la transformación que
>> ha sugerido Olivier, pero he obtenido el mismo tipo de errores que tú.
>>
>> >
>> >
names(datos)<-c(''Longitud'',''Densidad'')
>> > T<-1595
>> >
>> > nls.control(maxiter = 150, tol = 1e-15, minFactor = 1/2024,
>> +             printEval = FALSE, warnOnly = FALSE)
>> $maxiter
>> [1] 150
>>
>> $tol
>> [1] 1e-15
>>
>> $minFactor
>> [1] 0.0004940711
>>
>>  $printEval
>> [1] FALSE
>>
>> $warnOnly
>> [1] FALSE
>>
>> >
>> > fun.ajust<-nls( Densidad ~ (b1/(Longitud)^3) /
(exp(b2/(Longitud*T)) - 1
>> ), data=datos,
>> +                 trace=TRUE, start=list(b1 = 1, b2 = 7) )
>> Error en nlsModel(formula, mf, start, wts) :
>>   singular gradient matrix at initial parameter estimates
>> >
>>
>> En la función "nls.control()" he reducido el paso, e
incrementado el
>> número posible de iteraciones para poder aumentar la precisión del
>> cálculo... pero como ves obtengo error.
>> Igualmente, ves que he comenzado desde otros puntos b1=1 y b2=7.
>> He probado con b1=100 y b2=700, que son valores próximo al otro cálculo
y
>> el error es el mismo.
>>
>>
>>
>> Saludos,
>> Carlos Ortega
>> www.qualityexcellence.es
>> http://www.datanalytics.com/blog/
>>
>> 2010/12/5 Daniel Arismendi <darismendi@cmc.org.ve>
>>
>> Efectivamente con la forma en que tu presentas esto cuadra
perfectamente
>>> con
>>> el detalle de que  es (b1/(Longitud*r)^5) y lo otro T no puede ser
1 pues
>>> nunca encontraría el espectro de radiación de cuerpo negro, pero no
>>> importa
>>> yo le cambie (exp(b2/(Longitud*r*T) donde T es 1595.
>>>
>>> Pues cuando intento obtener h se que b1=8*pi*h*c
>>> en este caso el valor de b1 que obtengo es de 42.19
>>> y para b2 es 1388.23 donde b2 = hc/k
>>>
>>> Siendo k la constante de Boltzmann c la velocidad de la luz y h la
>>> constante
>>> de Planck que necesito hallar.
>>>
>>> Si la obtengo de b1  tengo un valor de h= 6e-9
>>> Si la obtengo de b2 h= 6.388868e-29
>>>
>>> Si valor a obtener debería ser próximo a  h= 6.2606896e-34 cosa que
en
>>> mis
>>> resultados anteriores es parecido pero los ordenes de magnitud son
muy
>>> diferentes.
>>>
>>> Ahora mi pregunta es: hay alguna forma de que R trabaje los valores
de
>>> mis
>>> longitud de onda con una precisión mucho mayor? para que pueda
trabajar
>>> con
>>> la regresión lineal de forma normal.
>>>
>>> Pues al trabajar con los valores multiplicados por 10^5 no es lo
mejor y
>>> pues como decías tu la interpretación de los valores de mis
parámetros
>>> cambia pero se debe a que estamos introduciendo errores
sistemáticas
>>> asociados a un calculo indirecto para encontrar la regresión no
lineal.
>>>
>>> Otra cosa Olivier como hiciste para darle los parámetros iniciales
a la
>>> regresión no lineal? pues veo que pusiste b1=.01 y b2=.02 los
encontraste
>>> de
>>> forma empírica o llegaste a un valor óptimo asociándolos a las dos
>>> maneras
>>> de encontrar los parámetros iniciales como lo propuse en mi mensaje
>>> inicial.?
>>>
>>>
>>> De igual forma gracias nuevamente por tu respuesta
>>> Saludos Olivier.
>>>
>>> El 5 de diciembre de 2010 16:58, Olivier Nuñez
<onunez@iberstat.es
>>> >escribió:
>>>
>>> > 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@iberstat.es
>>> > Tel : +34 663 03 69 09
>>> > Web: http://www.iberstat.es
>>> >
>>> > ____________________________________
>>> >
>>> >
>>> >
>>> >
>>> >
>>> > 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@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@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<http://en.wikipedia.org/wiki/Planck%27s_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@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@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@r-project.org
>>> >> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>> >>
>>> >>
>>> >>
>>> >>
>>> >>
>>> >
>>> >
>>>
>>>        [[alternative HTML version deleted]]
>>>
>>>
>>> _______________________________________________
>>> R-help-es mailing list
>>> R-help-es@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>>
>>>
>>
>
	[[alternative HTML version deleted]]