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

